adjointExtra_d.f90

Go to the documentation of this file.

!        generated by tapenade     (inria, ecuador team)
!  tapenade 3.16 (develop) - 22 aug 2023 15:51
!
module adjointextra_d
  implicit none
 
contains
!  differentiation of volume_block in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *vol
!   with respect to varying inputs: *x
!   rw status of diff variables: *x:in *vol:out
!   plus diff mem management of: x:in vol:in
  subroutine volume_block_d()
! this is copy of metric.f90. it was necessary to copy this file
! since there is debugging stuff in the original that is not
! necessary for ad.
    use constants
    use blockpointers
    use cgnsgrid
    use communication
    use inputtimespectral
    implicit none
!
!      local parameter.
!
    real(kind=realtype), parameter :: thresvolume=1.e-2_realtype
    real(kind=realtype), parameter :: halocellratio=1e-10_realtype
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, n, m, l, ii
    integer(kind=inttype) :: mm
    real(kind=realtype) :: fact, mult
    real(kind=realtype) :: xp, yp, zp, vp1, vp2, vp3, vp4, vp5, vp6
    real(kind=realtype) :: xpd, ypd, zpd, vp1d, vp2d, vp3d, vp4d, vp5d, &
&   vp6d
    real(kind=realtype) :: xxp, yyp, zzp
    real(kind=realtype), dimension(3) :: v1, v2
    intrinsic abs
! compute the volumes. the hexahedron is split into 6 pyramids
! whose volumes are computed. the volume is positive for a
! right handed block.
! initialize the volumes to zero. the reasons is that the second
! level halo's must be initialized to zero and for convenience
! all the volumes are set to zero.
    vold = 0.0_8
    vol = zero
    if (associated(vold)) vold = 0.0_8
    vp1d = 0.0_8
    vp2d = 0.0_8
    vp3d = 0.0_8
    vp4d = 0.0_8
    vp5d = 0.0_8
    vp6d = 0.0_8
    do k=1,ke
      n = k - 1
      do j=1,je
        m = j - 1
        do i=1,ie
          l = i - 1
! compute the coordinates of the center of gravity.
          xpd = eighth*(xd(i, j, k, 1)+xd(i, m, k, 1)+xd(i, m, n, 1)+xd(&
&           i, j, n, 1)+xd(l, j, k, 1)+xd(l, m, k, 1)+xd(l, m, n, 1)+xd(&
&           l, j, n, 1))
          xp = eighth*(x(i, j, k, 1)+x(i, m, k, 1)+x(i, m, n, 1)+x(i, j&
&           , n, 1)+x(l, j, k, 1)+x(l, m, k, 1)+x(l, m, n, 1)+x(l, j, n&
&           , 1))
          ypd = eighth*(xd(i, j, k, 2)+xd(i, m, k, 2)+xd(i, m, n, 2)+xd(&
&           i, j, n, 2)+xd(l, j, k, 2)+xd(l, m, k, 2)+xd(l, m, n, 2)+xd(&
&           l, j, n, 2))
          yp = eighth*(x(i, j, k, 2)+x(i, m, k, 2)+x(i, m, n, 2)+x(i, j&
&           , n, 2)+x(l, j, k, 2)+x(l, m, k, 2)+x(l, m, n, 2)+x(l, j, n&
&           , 2))
          zpd = eighth*(xd(i, j, k, 3)+xd(i, m, k, 3)+xd(i, m, n, 3)+xd(&
&           i, j, n, 3)+xd(l, j, k, 3)+xd(l, m, k, 3)+xd(l, m, n, 3)+xd(&
&           l, j, n, 3))
          zp = eighth*(x(i, j, k, 3)+x(i, m, k, 3)+x(i, m, n, 3)+x(i, j&
&           , n, 3)+x(l, j, k, 3)+x(l, m, k, 3)+x(l, m, n, 3)+x(l, j, n&
&           , 3))
! compute the volumes of the 6 sub pyramids. the
! arguments of volpym must be such that for a (regular)
! right handed hexahedron all volumes are positive.
          call volpym_d(x(i, j, k, 1), xd(i, j, k, 1), x(i, j, k, 2), xd&
&                 (i, j, k, 2), x(i, j, k, 3), xd(i, j, k, 3), x(i, j, n&
&                 , 1), xd(i, j, n, 1), x(i, j, n, 2), xd(i, j, n, 2), x&
&                 (i, j, n, 3), xd(i, j, n, 3), x(i, m, n, 1), xd(i, m, &
&                 n, 1), x(i, m, n, 2), xd(i, m, n, 2), x(i, m, n, 3), &
&                 xd(i, m, n, 3), x(i, m, k, 1), xd(i, m, k, 1), x(i, m&
&                 , k, 2), xd(i, m, k, 2), x(i, m, k, 3), xd(i, m, k, 3)&
&                 , vp1, vp1d)
          call volpym_d(x(l, j, k, 1), xd(l, j, k, 1), x(l, j, k, 2), xd&
&                 (l, j, k, 2), x(l, j, k, 3), xd(l, j, k, 3), x(l, m, k&
&                 , 1), xd(l, m, k, 1), x(l, m, k, 2), xd(l, m, k, 2), x&
&                 (l, m, k, 3), xd(l, m, k, 3), x(l, m, n, 1), xd(l, m, &
&                 n, 1), x(l, m, n, 2), xd(l, m, n, 2), x(l, m, n, 3), &
&                 xd(l, m, n, 3), x(l, j, n, 1), xd(l, j, n, 1), x(l, j&
&                 , n, 2), xd(l, j, n, 2), x(l, j, n, 3), xd(l, j, n, 3)&
&                 , vp2, vp2d)
          call volpym_d(x(i, j, k, 1), xd(i, j, k, 1), x(i, j, k, 2), xd&
&                 (i, j, k, 2), x(i, j, k, 3), xd(i, j, k, 3), x(l, j, k&
&                 , 1), xd(l, j, k, 1), x(l, j, k, 2), xd(l, j, k, 2), x&
&                 (l, j, k, 3), xd(l, j, k, 3), x(l, j, n, 1), xd(l, j, &
&                 n, 1), x(l, j, n, 2), xd(l, j, n, 2), x(l, j, n, 3), &
&                 xd(l, j, n, 3), x(i, j, n, 1), xd(i, j, n, 1), x(i, j&
&                 , n, 2), xd(i, j, n, 2), x(i, j, n, 3), xd(i, j, n, 3)&
&                 , vp3, vp3d)
          call volpym_d(x(i, m, k, 1), xd(i, m, k, 1), x(i, m, k, 2), xd&
&                 (i, m, k, 2), x(i, m, k, 3), xd(i, m, k, 3), x(i, m, n&
&                 , 1), xd(i, m, n, 1), x(i, m, n, 2), xd(i, m, n, 2), x&
&                 (i, m, n, 3), xd(i, m, n, 3), x(l, m, n, 1), xd(l, m, &
&                 n, 1), x(l, m, n, 2), xd(l, m, n, 2), x(l, m, n, 3), &
&                 xd(l, m, n, 3), x(l, m, k, 1), xd(l, m, k, 1), x(l, m&
&                 , k, 2), xd(l, m, k, 2), x(l, m, k, 3), xd(l, m, k, 3)&
&                 , vp4, vp4d)
          call volpym_d(x(i, j, k, 1), xd(i, j, k, 1), x(i, j, k, 2), xd&
&                 (i, j, k, 2), x(i, j, k, 3), xd(i, j, k, 3), x(i, m, k&
&                 , 1), xd(i, m, k, 1), x(i, m, k, 2), xd(i, m, k, 2), x&
&                 (i, m, k, 3), xd(i, m, k, 3), x(l, m, k, 1), xd(l, m, &
&                 k, 1), x(l, m, k, 2), xd(l, m, k, 2), x(l, m, k, 3), &
&                 xd(l, m, k, 3), x(l, j, k, 1), xd(l, j, k, 1), x(l, j&
&                 , k, 2), xd(l, j, k, 2), x(l, j, k, 3), xd(l, j, k, 3)&
&                 , vp5, vp5d)
          call volpym_d(x(i, j, n, 1), xd(i, j, n, 1), x(i, j, n, 2), xd&
&                 (i, j, n, 2), x(i, j, n, 3), xd(i, j, n, 3), x(l, j, n&
&                 , 1), xd(l, j, n, 1), x(l, j, n, 2), xd(l, j, n, 2), x&
&                 (l, j, n, 3), xd(l, j, n, 3), x(l, m, n, 1), xd(l, m, &
&                 n, 1), x(l, m, n, 2), xd(l, m, n, 2), x(l, m, n, 3), &
&                 xd(l, m, n, 3), x(i, m, n, 1), xd(i, m, n, 1), x(i, m&
&                 , n, 2), xd(i, m, n, 2), x(i, m, n, 3), xd(i, m, n, 3)&
&                 , vp6, vp6d)
! set the volume to 1/6 of the sum of the volumes of the
! pyramid. remember that volpym computes 6 times the
! volume.
          vold(i, j, k) = sixth*(vp1d+vp2d+vp3d+vp4d+vp5d+vp6d)
          vol(i, j, k) = sixth*(vp1+vp2+vp3+vp4+vp5+vp6)
          if (vol(i, j, k) .ge. 0.) then
            vol(i, j, k) = vol(i, j, k)
          else
            vold(i, j, k) = -vold(i, j, k)
            vol(i, j, k) = -vol(i, j, k)
          end if
        end do
      end do
    end do
! some additional safety stuff for halo volumes.
    do k=2,kl
      do j=2,jl
        if (vol(1, j, k)/vol(2, j, k) .lt. halocellratio) then
          vold(1, j, k) = vold(2, j, k)
          vol(1, j, k) = vol(2, j, k)
        end if
        if (vol(ie, j, k)/vol(il, j, k) .lt. halocellratio) then
          vold(ie, j, k) = vold(il, j, k)
          vol(ie, j, k) = vol(il, j, k)
        end if
      end do
    end do
    do k=2,kl
      do i=1,ie
        if (vol(i, 1, k)/vol(i, 2, k) .lt. halocellratio) then
          vold(i, 1, k) = vold(i, 2, k)
          vol(i, 1, k) = vol(i, 2, k)
        end if
        if (vol(i, je, k)/vol(i, jl, k) .lt. halocellratio) then
          vold(i, je, k) = vold(i, jl, k)
          vol(i, je, k) = vol(i, jl, k)
        end if
      end do
    end do
    do j=1,je
      do i=1,ie
        if (vol(i, j, 1)/vol(i, j, 2) .lt. halocellratio) then
          vold(i, j, 1) = vold(i, j, 2)
          vol(i, j, 1) = vol(i, j, 2)
        end if
        if (vol(i, j, ke)/vol(i, j, kl) .lt. halocellratio) then
          vold(i, j, ke) = vold(i, j, kl)
          vol(i, j, ke) = vol(i, j, kl)
        end if
      end do
    end do
 
  contains
!  differentiation of volpym in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: volume
!   with respect to varying inputs: xp yp zp xa xb xc xd ya yb
!                yc yd za zb zc zd
    subroutine volpym_d(xa, xad, ya, yad, za, zad, xb, xbd, yb, ybd, zb&
&     , zbd, xc, xcd, yc, ycd, zc, zcd, xd, xdd, yd, ydd, zd, zdd, &
&     volume, volumed)
!
!         volpym computes 6 times the volume of a pyramid. node p,
!         whose coordinates are set in the subroutine metric itself,
!         is the top node and a-b-c-d is the quadrilateral surface.
!         it is assumed that the cross product vca * vdb points in
!         the direction of the top node. here vca is the diagonal
!         running from node c to node a and vdb the diagonal from
!         node d to node b.
!
      use precision
      implicit none
!
!        function type.
!
      real(kind=realtype) :: volume
      real(kind=realtype) :: volumed
!
!        function arguments.
!
      real(kind=realtype), intent(in) :: xa, ya, za, xb, yb, zb
      real(kind=realtype), intent(in) :: xad, yad, zad, xbd, ybd, zbd
      real(kind=realtype), intent(in) :: xc, yc, zc, xd, yd, zd
      real(kind=realtype), intent(in) :: xcd, ycd, zcd, xdd, ydd, zdd
      real(kind=realtype) :: temp
      real(kind=realtype) :: temp0
      real(kind=realtype) :: temp1
      real(kind=realtype) :: temp2
      real(kind=realtype) :: temp3
      real(kind=realtype) :: temp4
      temp = (ya-yc)*(zb-zd) - (za-zc)*(yb-yd)
      temp0 = xp - fourth*(xa+xb+xc+xd)
      temp1 = (za-zc)*(xb-xd) - (xa-xc)*(zb-zd)
      temp2 = yp - fourth*(ya+yb+yc+yd)
      temp3 = (xa-xc)*(yb-yd) - (ya-yc)*(xb-xd)
      temp4 = zp - fourth*(za+zb+zc+zd)
      volumed = temp*(xpd-fourth*(xad+xbd+xcd+xdd)) + temp0*((zb-zd)*(&
&       yad-ycd)+(ya-yc)*(zbd-zdd)-(yb-yd)*(zad-zcd)-(za-zc)*(ybd-ydd)) &
&       + temp1*(ypd-fourth*(yad+ybd+ycd+ydd)) + temp2*((xb-xd)*(zad-zcd&
&       )+(za-zc)*(xbd-xdd)-(zb-zd)*(xad-xcd)-(xa-xc)*(zbd-zdd)) + temp3&
&       *(zpd-fourth*(zad+zbd+zcd+zdd)) + temp4*((yb-yd)*(xad-xcd)+(xa-&
&       xc)*(ybd-ydd)-(xb-xd)*(yad-ycd)-(ya-yc)*(xbd-xdd))
      volume = temp0*temp + temp2*temp1 + temp4*temp3
    end subroutine volpym_d
 
    subroutine volpym(xa, ya, za, xb, yb, zb, xc, yc, zc, xd, yd, zd, &
&     volume)
!
!         volpym computes 6 times the volume of a pyramid. node p,
!         whose coordinates are set in the subroutine metric itself,
!         is the top node and a-b-c-d is the quadrilateral surface.
!         it is assumed that the cross product vca * vdb points in
!         the direction of the top node. here vca is the diagonal
!         running from node c to node a and vdb the diagonal from
!         node d to node b.
!
      use precision
      implicit none
!
!        function type.
!
      real(kind=realtype) :: volume
!
!        function arguments.
!
      real(kind=realtype), intent(in) :: xa, ya, za, xb, yb, zb
      real(kind=realtype), intent(in) :: xc, yc, zc, xd, yd, zd
      volume = (xp-fourth*(xa+xb+xc+xd))*((ya-yc)*(zb-zd)-(za-zc)*(yb-yd&
&       )) + (yp-fourth*(ya+yb+yc+yd))*((za-zc)*(xb-xd)-(xa-xc)*(zb-zd))&
&       + (zp-fourth*(za+zb+zc+zd))*((xa-xc)*(yb-yd)-(ya-yc)*(xb-xd))
    end subroutine volpym
 
  end subroutine volume_block_d
 
  subroutine volume_block()
! this is copy of metric.f90. it was necessary to copy this file
! since there is debugging stuff in the original that is not
! necessary for ad.
    use constants
    use blockpointers
    use cgnsgrid
    use communication
    use inputtimespectral
    implicit none
!
!      local parameter.
!
    real(kind=realtype), parameter :: thresvolume=1.e-2_realtype
    real(kind=realtype), parameter :: halocellratio=1e-10_realtype
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, n, m, l, ii
    integer(kind=inttype) :: mm
    real(kind=realtype) :: fact, mult
    real(kind=realtype) :: xp, yp, zp, vp1, vp2, vp3, vp4, vp5, vp6
    real(kind=realtype) :: xxp, yyp, zzp
    real(kind=realtype), dimension(3) :: v1, v2
    intrinsic abs
! compute the volumes. the hexahedron is split into 6 pyramids
! whose volumes are computed. the volume is positive for a
! right handed block.
! initialize the volumes to zero. the reasons is that the second
! level halo's must be initialized to zero and for convenience
! all the volumes are set to zero.
    vol = zero
    do k=1,ke
      n = k - 1
      do j=1,je
        m = j - 1
        do i=1,ie
          l = i - 1
! compute the coordinates of the center of gravity.
          xp = eighth*(x(i, j, k, 1)+x(i, m, k, 1)+x(i, m, n, 1)+x(i, j&
&           , n, 1)+x(l, j, k, 1)+x(l, m, k, 1)+x(l, m, n, 1)+x(l, j, n&
&           , 1))
          yp = eighth*(x(i, j, k, 2)+x(i, m, k, 2)+x(i, m, n, 2)+x(i, j&
&           , n, 2)+x(l, j, k, 2)+x(l, m, k, 2)+x(l, m, n, 2)+x(l, j, n&
&           , 2))
          zp = eighth*(x(i, j, k, 3)+x(i, m, k, 3)+x(i, m, n, 3)+x(i, j&
&           , n, 3)+x(l, j, k, 3)+x(l, m, k, 3)+x(l, m, n, 3)+x(l, j, n&
&           , 3))
! compute the volumes of the 6 sub pyramids. the
! arguments of volpym must be such that for a (regular)
! right handed hexahedron all volumes are positive.
          call volpym(x(i, j, k, 1), x(i, j, k, 2), x(i, j, k, 3), x(i, &
&               j, n, 1), x(i, j, n, 2), x(i, j, n, 3), x(i, m, n, 1), x&
&               (i, m, n, 2), x(i, m, n, 3), x(i, m, k, 1), x(i, m, k, 2&
&               ), x(i, m, k, 3), vp1)
          call volpym(x(l, j, k, 1), x(l, j, k, 2), x(l, j, k, 3), x(l, &
&               m, k, 1), x(l, m, k, 2), x(l, m, k, 3), x(l, m, n, 1), x&
&               (l, m, n, 2), x(l, m, n, 3), x(l, j, n, 1), x(l, j, n, 2&
&               ), x(l, j, n, 3), vp2)
          call volpym(x(i, j, k, 1), x(i, j, k, 2), x(i, j, k, 3), x(l, &
&               j, k, 1), x(l, j, k, 2), x(l, j, k, 3), x(l, j, n, 1), x&
&               (l, j, n, 2), x(l, j, n, 3), x(i, j, n, 1), x(i, j, n, 2&
&               ), x(i, j, n, 3), vp3)
          call volpym(x(i, m, k, 1), x(i, m, k, 2), x(i, m, k, 3), x(i, &
&               m, n, 1), x(i, m, n, 2), x(i, m, n, 3), x(l, m, n, 1), x&
&               (l, m, n, 2), x(l, m, n, 3), x(l, m, k, 1), x(l, m, k, 2&
&               ), x(l, m, k, 3), vp4)
          call volpym(x(i, j, k, 1), x(i, j, k, 2), x(i, j, k, 3), x(i, &
&               m, k, 1), x(i, m, k, 2), x(i, m, k, 3), x(l, m, k, 1), x&
&               (l, m, k, 2), x(l, m, k, 3), x(l, j, k, 1), x(l, j, k, 2&
&               ), x(l, j, k, 3), vp5)
          call volpym(x(i, j, n, 1), x(i, j, n, 2), x(i, j, n, 3), x(l, &
&               j, n, 1), x(l, j, n, 2), x(l, j, n, 3), x(l, m, n, 1), x&
&               (l, m, n, 2), x(l, m, n, 3), x(i, m, n, 1), x(i, m, n, 2&
&               ), x(i, m, n, 3), vp6)
! set the volume to 1/6 of the sum of the volumes of the
! pyramid. remember that volpym computes 6 times the
! volume.
          vol(i, j, k) = sixth*(vp1+vp2+vp3+vp4+vp5+vp6)
          if (vol(i, j, k) .ge. 0.) then
            vol(i, j, k) = vol(i, j, k)
          else
            vol(i, j, k) = -vol(i, j, k)
          end if
        end do
      end do
    end do
! some additional safety stuff for halo volumes.
    do k=2,kl
      do j=2,jl
        if (vol(1, j, k)/vol(2, j, k) .lt. halocellratio) vol(1, j, k)&
&          = vol(2, j, k)
        if (vol(ie, j, k)/vol(il, j, k) .lt. halocellratio) vol(ie, j, k&
&         ) = vol(il, j, k)
      end do
    end do
    do k=2,kl
      do i=1,ie
        if (vol(i, 1, k)/vol(i, 2, k) .lt. halocellratio) vol(i, 1, k)&
&          = vol(i, 2, k)
        if (vol(i, je, k)/vol(i, jl, k) .lt. halocellratio) vol(i, je, k&
&         ) = vol(i, jl, k)
      end do
    end do
    do j=1,je
      do i=1,ie
        if (vol(i, j, 1)/vol(i, j, 2) .lt. halocellratio) vol(i, j, 1)&
&          = vol(i, j, 2)
        if (vol(i, j, ke)/vol(i, j, kl) .lt. halocellratio) vol(i, j, ke&
&         ) = vol(i, j, kl)
      end do
    end do
 
  contains
    subroutine volpym(xa, ya, za, xb, yb, zb, xc, yc, zc, xd, yd, zd, &
&     volume)
!
!         volpym computes 6 times the volume of a pyramid. node p,
!         whose coordinates are set in the subroutine metric itself,
!         is the top node and a-b-c-d is the quadrilateral surface.
!         it is assumed that the cross product vca * vdb points in
!         the direction of the top node. here vca is the diagonal
!         running from node c to node a and vdb the diagonal from
!         node d to node b.
!
      use precision
      implicit none
!
!        function type.
!
      real(kind=realtype) :: volume
!
!        function arguments.
!
      real(kind=realtype), intent(in) :: xa, ya, za, xb, yb, zb
      real(kind=realtype), intent(in) :: xc, yc, zc, xd, yd, zd
      volume = (xp-fourth*(xa+xb+xc+xd))*((ya-yc)*(zb-zd)-(za-zc)*(yb-yd&
&       )) + (yp-fourth*(ya+yb+yc+yd))*((za-zc)*(xb-xd)-(xa-xc)*(zb-zd))&
&       + (zp-fourth*(za+zb+zc+zd))*((xa-xc)*(yb-yd)-(ya-yc)*(xb-xd))
    end subroutine volpym
 
  end subroutine volume_block
 
!  differentiation of metric_block in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *si *sj *sk
!   with respect to varying inputs: *x
!   rw status of diff variables: *x:in *si:out *sj:out *sk:out
!   plus diff mem management of: x:in si:in sj:in sk:in
  subroutine metric_block_d()
    use constants
    use blockpointers
    implicit none
! local variables.
    integer(kind=inttype) :: i, j, k, n, m, l, ii
    real(kind=realtype) :: fact
    real(kind=realtype) :: xxp, yyp, zzp
    real(kind=realtype), dimension(3) :: v1, v2
    real(kind=realtype), dimension(3) :: v1d, v2d
    intrinsic mod
! set the factor in the surface normals computation. for a
! left handed block this factor is negative, such that the
! normals still point in the direction of increasing index.
! the formulae used later on assume a right handed block
! and fact is used to correct this for a left handed block,
! as well as the scaling factor of 0.5
    if (righthanded) then
      fact = half
    else
      fact = -half
    end if
    if (associated(sid)) sid = 0.0_8
    v1d = 0.0_8
    v2d = 0.0_8
!
!  computation of the face normals in i-, j- and k-direction.
!  formula's are valid for a right handed block; for a left
!  handed block the correct orientation is obtained via fact.
!  the normals point in the direction of increasing index.
!  the absolute value of fact is 0.5, because the cross
!  product of the two diagonals is twice the normal vector.
!  note that also the normals of the first level halo cells
!  are computed. these are needed for the viscous fluxes.
!
! projected areas of cell faces in the i direction.
    do ii=0,ke*je*(ie+1)-1
! 0:ie
      i = mod(ii, ie + 1) + 0
!1:je
      j = mod(ii/(ie+1), je) + 1
!1:ke
      k = ii/((ie+1)*je) + 1
      n = k - 1
      m = j - 1
! determine the two diagonal vectors of the face.
      v1d(1) = xd(i, j, n, 1) - xd(i, m, k, 1)
      v1(1) = x(i, j, n, 1) - x(i, m, k, 1)
      v1d(2) = xd(i, j, n, 2) - xd(i, m, k, 2)
      v1(2) = x(i, j, n, 2) - x(i, m, k, 2)
      v1d(3) = xd(i, j, n, 3) - xd(i, m, k, 3)
      v1(3) = x(i, j, n, 3) - x(i, m, k, 3)
      v2d(1) = xd(i, j, k, 1) - xd(i, m, n, 1)
      v2(1) = x(i, j, k, 1) - x(i, m, n, 1)
      v2d(2) = xd(i, j, k, 2) - xd(i, m, n, 2)
      v2(2) = x(i, j, k, 2) - x(i, m, n, 2)
      v2d(3) = xd(i, j, k, 3) - xd(i, m, n, 3)
      v2(3) = x(i, j, k, 3) - x(i, m, n, 3)
! the face normal, which is the cross product of the two
! diagonal vectors times fact; remember that fact is
! either -0.5 or 0.5.
      sid(i, j, k, 1) = fact*(v2(3)*v1d(2)+v1(2)*v2d(3)-v2(2)*v1d(3)-v1(&
&       3)*v2d(2))
      si(i, j, k, 1) = fact*(v1(2)*v2(3)-v1(3)*v2(2))
      sid(i, j, k, 2) = fact*(v2(1)*v1d(3)+v1(3)*v2d(1)-v2(3)*v1d(1)-v1(&
&       1)*v2d(3))
      si(i, j, k, 2) = fact*(v1(3)*v2(1)-v1(1)*v2(3))
      sid(i, j, k, 3) = fact*(v2(2)*v1d(1)+v1(1)*v2d(2)-v2(1)*v1d(2)-v1(&
&       2)*v2d(1))
      si(i, j, k, 3) = fact*(v1(1)*v2(2)-v1(2)*v2(1))
    end do
    if (associated(sjd)) sjd = 0.0_8
! projected areas of cell faces in the j direction
    do ii=0,ke*(je+1)*ie-1
! 1:ie
      i = mod(ii, ie) + 1
!0:je
      j = mod(ii/ie, je + 1) + 0
!1:ke
      k = ii/(ie*(je+1)) + 1
      n = k - 1
      l = i - 1
! determine the two diagonal vectors of the face.
      v1d(1) = xd(i, j, n, 1) - xd(l, j, k, 1)
      v1(1) = x(i, j, n, 1) - x(l, j, k, 1)
      v1d(2) = xd(i, j, n, 2) - xd(l, j, k, 2)
      v1(2) = x(i, j, n, 2) - x(l, j, k, 2)
      v1d(3) = xd(i, j, n, 3) - xd(l, j, k, 3)
      v1(3) = x(i, j, n, 3) - x(l, j, k, 3)
      v2d(1) = xd(l, j, n, 1) - xd(i, j, k, 1)
      v2(1) = x(l, j, n, 1) - x(i, j, k, 1)
      v2d(2) = xd(l, j, n, 2) - xd(i, j, k, 2)
      v2(2) = x(l, j, n, 2) - x(i, j, k, 2)
      v2d(3) = xd(l, j, n, 3) - xd(i, j, k, 3)
      v2(3) = x(l, j, n, 3) - x(i, j, k, 3)
! the face normal, which is the cross product of the two
! diagonal vectors times fact; remember that fact is
! either -0.5 or 0.5.
      sjd(i, j, k, 1) = fact*(v2(3)*v1d(2)+v1(2)*v2d(3)-v2(2)*v1d(3)-v1(&
&       3)*v2d(2))
      sj(i, j, k, 1) = fact*(v1(2)*v2(3)-v1(3)*v2(2))
      sjd(i, j, k, 2) = fact*(v2(1)*v1d(3)+v1(3)*v2d(1)-v2(3)*v1d(1)-v1(&
&       1)*v2d(3))
      sj(i, j, k, 2) = fact*(v1(3)*v2(1)-v1(1)*v2(3))
      sjd(i, j, k, 3) = fact*(v2(2)*v1d(1)+v1(1)*v2d(2)-v2(1)*v1d(2)-v1(&
&       2)*v2d(1))
      sj(i, j, k, 3) = fact*(v1(1)*v2(2)-v1(2)*v2(1))
    end do
    if (associated(skd)) skd = 0.0_8
! projected areas of cell faces in the k direction.
    do ii=0,(ke+1)*je*ie-1
! 1:ie
      i = mod(ii, ie) + 1
!1:je
      j = mod(ii/ie, je) + 1
!0:ke
      k = ii/(ie*je) + 0
      m = j - 1
      l = i - 1
! determine the two diagonal vectors of the face.
      v1d(1) = xd(i, j, k, 1) - xd(l, m, k, 1)
      v1(1) = x(i, j, k, 1) - x(l, m, k, 1)
      v1d(2) = xd(i, j, k, 2) - xd(l, m, k, 2)
      v1(2) = x(i, j, k, 2) - x(l, m, k, 2)
      v1d(3) = xd(i, j, k, 3) - xd(l, m, k, 3)
      v1(3) = x(i, j, k, 3) - x(l, m, k, 3)
      v2d(1) = xd(l, j, k, 1) - xd(i, m, k, 1)
      v2(1) = x(l, j, k, 1) - x(i, m, k, 1)
      v2d(2) = xd(l, j, k, 2) - xd(i, m, k, 2)
      v2(2) = x(l, j, k, 2) - x(i, m, k, 2)
      v2d(3) = xd(l, j, k, 3) - xd(i, m, k, 3)
      v2(3) = x(l, j, k, 3) - x(i, m, k, 3)
! the face normal, which is the cross product of the two
! diagonal vectors times fact; remember that fact is
! either -0.5 or 0.5.
      skd(i, j, k, 1) = fact*(v2(3)*v1d(2)+v1(2)*v2d(3)-v2(2)*v1d(3)-v1(&
&       3)*v2d(2))
      sk(i, j, k, 1) = fact*(v1(2)*v2(3)-v1(3)*v2(2))
      skd(i, j, k, 2) = fact*(v2(1)*v1d(3)+v1(3)*v2d(1)-v2(3)*v1d(1)-v1(&
&       1)*v2d(3))
      sk(i, j, k, 2) = fact*(v1(3)*v2(1)-v1(1)*v2(3))
      skd(i, j, k, 3) = fact*(v2(2)*v1d(1)+v1(1)*v2d(2)-v2(1)*v1d(2)-v1(&
&       2)*v2d(1))
      sk(i, j, k, 3) = fact*(v1(1)*v2(2)-v1(2)*v2(1))
    end do
  end subroutine metric_block_d
 
  subroutine metric_block()
    use constants
    use blockpointers
    implicit none
! local variables.
    integer(kind=inttype) :: i, j, k, n, m, l, ii
    real(kind=realtype) :: fact
    real(kind=realtype) :: xxp, yyp, zzp
    real(kind=realtype), dimension(3) :: v1, v2
    intrinsic mod
! set the factor in the surface normals computation. for a
! left handed block this factor is negative, such that the
! normals still point in the direction of increasing index.
! the formulae used later on assume a right handed block
! and fact is used to correct this for a left handed block,
! as well as the scaling factor of 0.5
    if (righthanded) then
      fact = half
    else
      fact = -half
    end if
!$ad ii-loop
!
!  computation of the face normals in i-, j- and k-direction.
!  formula's are valid for a right handed block; for a left
!  handed block the correct orientation is obtained via fact.
!  the normals point in the direction of increasing index.
!  the absolute value of fact is 0.5, because the cross
!  product of the two diagonals is twice the normal vector.
!  note that also the normals of the first level halo cells
!  are computed. these are needed for the viscous fluxes.
!
! projected areas of cell faces in the i direction.
    do ii=0,ke*je*(ie+1)-1
! 0:ie
      i = mod(ii, ie + 1) + 0
!1:je
      j = mod(ii/(ie+1), je) + 1
!1:ke
      k = ii/((ie+1)*je) + 1
      n = k - 1
      m = j - 1
! determine the two diagonal vectors of the face.
      v1(1) = x(i, j, n, 1) - x(i, m, k, 1)
      v1(2) = x(i, j, n, 2) - x(i, m, k, 2)
      v1(3) = x(i, j, n, 3) - x(i, m, k, 3)
      v2(1) = x(i, j, k, 1) - x(i, m, n, 1)
      v2(2) = x(i, j, k, 2) - x(i, m, n, 2)
      v2(3) = x(i, j, k, 3) - x(i, m, n, 3)
! the face normal, which is the cross product of the two
! diagonal vectors times fact; remember that fact is
! either -0.5 or 0.5.
      si(i, j, k, 1) = fact*(v1(2)*v2(3)-v1(3)*v2(2))
      si(i, j, k, 2) = fact*(v1(3)*v2(1)-v1(1)*v2(3))
      si(i, j, k, 3) = fact*(v1(1)*v2(2)-v1(2)*v2(1))
    end do
!$ad ii-loop
! projected areas of cell faces in the j direction
    do ii=0,ke*(je+1)*ie-1
! 1:ie
      i = mod(ii, ie) + 1
!0:je
      j = mod(ii/ie, je + 1) + 0
!1:ke
      k = ii/(ie*(je+1)) + 1
      n = k - 1
      l = i - 1
! determine the two diagonal vectors of the face.
      v1(1) = x(i, j, n, 1) - x(l, j, k, 1)
      v1(2) = x(i, j, n, 2) - x(l, j, k, 2)
      v1(3) = x(i, j, n, 3) - x(l, j, k, 3)
      v2(1) = x(l, j, n, 1) - x(i, j, k, 1)
      v2(2) = x(l, j, n, 2) - x(i, j, k, 2)
      v2(3) = x(l, j, n, 3) - x(i, j, k, 3)
! the face normal, which is the cross product of the two
! diagonal vectors times fact; remember that fact is
! either -0.5 or 0.5.
      sj(i, j, k, 1) = fact*(v1(2)*v2(3)-v1(3)*v2(2))
      sj(i, j, k, 2) = fact*(v1(3)*v2(1)-v1(1)*v2(3))
      sj(i, j, k, 3) = fact*(v1(1)*v2(2)-v1(2)*v2(1))
    end do
!$ad ii-loop
! projected areas of cell faces in the k direction.
    do ii=0,(ke+1)*je*ie-1
! 1:ie
      i = mod(ii, ie) + 1
!1:je
      j = mod(ii/ie, je) + 1
!0:ke
      k = ii/(ie*je) + 0
      m = j - 1
      l = i - 1
! determine the two diagonal vectors of the face.
      v1(1) = x(i, j, k, 1) - x(l, m, k, 1)
      v1(2) = x(i, j, k, 2) - x(l, m, k, 2)
      v1(3) = x(i, j, k, 3) - x(l, m, k, 3)
      v2(1) = x(l, j, k, 1) - x(i, m, k, 1)
      v2(2) = x(l, j, k, 2) - x(i, m, k, 2)
      v2(3) = x(l, j, k, 3) - x(i, m, k, 3)
! the face normal, which is the cross product of the two
! diagonal vectors times fact; remember that fact is
! either -0.5 or 0.5.
      sk(i, j, k, 1) = fact*(v1(2)*v2(3)-v1(3)*v2(2))
      sk(i, j, k, 2) = fact*(v1(3)*v2(1)-v1(1)*v2(3))
      sk(i, j, k, 3) = fact*(v1(1)*v2(2)-v1(2)*v2(1))
    end do
  end subroutine metric_block
 
!  differentiation of boundarynormals in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *(*bcdata.norm)
!   with respect to varying inputs: *si *sj *sk
!   rw status of diff variables: *si:in *sj:in *sk:in *(*bcdata.norm):out
!   plus diff mem management of: si:in sj:in sk:in bcdata:in *bcdata.norm:in
  subroutine boundarynormals_d()
!  the unit normals on the boundary faces. these always point
!  out of the domain, so a multiplication by -1 is needed for
!  the imin, jmin and kmin boundaries.
!
    use constants
    use blockpointers
    use cgnsgrid
    use communication
    use inputtimespectral
    use diffsizes
!  hint: isize1ofdrfbcdata should be the size of dimension 1 of array *bcdata
    implicit none
! local variables.
    integer(kind=inttype) :: i, j, ii
    integer(kind=inttype) :: mm
    real(kind=realtype) :: fact, mult
    real(kind=realtype) :: factd
    real(kind=realtype) :: xxp, yyp, zzp
    real(kind=realtype) :: xxpd, yypd, zzpd
    intrinsic mod
    intrinsic sqrt
    real(kind=realtype) :: arg1
    real(kind=realtype) :: arg1d
    real(kind=realtype) :: temp
    integer :: ii1
    if (associated(bcdatad)) then
      do ii1=1,isize1ofdrfbcdata
        bcdatad(ii1)%norm = 0.0_8
      end do
    end if
    zzpd = 0.0_8
    yypd = 0.0_8
    xxpd = 0.0_8
!loop over the boundary subfaces of this block.
bocoloop:do mm=1,nbocos
! loop over the boundary faces of the subface.
      do ii=0,(bcdata(mm)%jcend-bcdata(mm)%jcbeg+1)*(bcdata(mm)%icend-&
&         bcdata(mm)%icbeg+1)-1
        i = mod(ii, bcdata(mm)%icend - bcdata(mm)%icbeg + 1) + bcdata(mm&
&         )%icbeg
        j = ii/(bcdata(mm)%icend-bcdata(mm)%icbeg+1) + bcdata(mm)%jcbeg
        select case  (bcfaceid(mm)) 
        case (imin) 
          mult = -one
          xxpd = sid(1, i, j, 1)
          xxp = si(1, i, j, 1)
          yypd = sid(1, i, j, 2)
          yyp = si(1, i, j, 2)
          zzpd = sid(1, i, j, 3)
          zzp = si(1, i, j, 3)
        case (imax) 
          mult = one
          xxpd = sid(il, i, j, 1)
          xxp = si(il, i, j, 1)
          yypd = sid(il, i, j, 2)
          yyp = si(il, i, j, 2)
          zzpd = sid(il, i, j, 3)
          zzp = si(il, i, j, 3)
        case (jmin) 
          mult = -one
          xxpd = sjd(i, 1, j, 1)
          xxp = sj(i, 1, j, 1)
          yypd = sjd(i, 1, j, 2)
          yyp = sj(i, 1, j, 2)
          zzpd = sjd(i, 1, j, 3)
          zzp = sj(i, 1, j, 3)
        case (jmax) 
          mult = one
          xxpd = sjd(i, jl, j, 1)
          xxp = sj(i, jl, j, 1)
          yypd = sjd(i, jl, j, 2)
          yyp = sj(i, jl, j, 2)
          zzpd = sjd(i, jl, j, 3)
          zzp = sj(i, jl, j, 3)
        case (kmin) 
          mult = -one
          xxpd = skd(i, j, 1, 1)
          xxp = sk(i, j, 1, 1)
          yypd = skd(i, j, 1, 2)
          yyp = sk(i, j, 1, 2)
          zzpd = skd(i, j, 1, 3)
          zzp = sk(i, j, 1, 3)
        case (kmax) 
          mult = one
          xxpd = skd(i, j, kl, 1)
          xxp = sk(i, j, kl, 1)
          yypd = skd(i, j, kl, 2)
          yyp = sk(i, j, kl, 2)
          zzpd = skd(i, j, kl, 3)
          zzp = sk(i, j, kl, 3)
        end select
! compute the inverse of the length of the normal vector
! and possibly correct for inward pointing.
        arg1d = 2*xxp*xxpd + 2*yyp*yypd + 2*zzp*zzpd
        arg1 = xxp*xxp + yyp*yyp + zzp*zzp
        temp = sqrt(arg1)
        if (arg1 .eq. 0.0_8) then
          factd = 0.0_8
        else
          factd = arg1d/(2.0*temp)
        end if
        fact = temp
        if (fact .gt. zero) then
          factd = -(mult*factd/fact**2)
          fact = mult/fact
        end if
! compute the unit normal.
        bcdatad(mm)%norm(i, j, 1) = xxp*factd + fact*xxpd
        bcdata(mm)%norm(i, j, 1) = fact*xxp
        bcdatad(mm)%norm(i, j, 2) = yyp*factd + fact*yypd
        bcdata(mm)%norm(i, j, 2) = fact*yyp
        bcdatad(mm)%norm(i, j, 3) = zzp*factd + fact*zzpd
        bcdata(mm)%norm(i, j, 3) = fact*zzp
      end do
    end do bocoloop
  end subroutine boundarynormals_d
 
  subroutine boundarynormals()
!  the unit normals on the boundary faces. these always point
!  out of the domain, so a multiplication by -1 is needed for
!  the imin, jmin and kmin boundaries.
!
    use constants
    use blockpointers
    use cgnsgrid
    use communication
    use inputtimespectral
    implicit none
! local variables.
    integer(kind=inttype) :: i, j, ii
    integer(kind=inttype) :: mm
    real(kind=realtype) :: fact, mult
    real(kind=realtype) :: xxp, yyp, zzp
    intrinsic mod
    intrinsic sqrt
    real(kind=realtype) :: arg1
!$ad ii-loop
!loop over the boundary subfaces of this block.
bocoloop:do mm=1,nbocos
!$ad ii-loop
! loop over the boundary faces of the subface.
      do ii=0,(bcdata(mm)%jcend-bcdata(mm)%jcbeg+1)*(bcdata(mm)%icend-&
&         bcdata(mm)%icbeg+1)-1
        i = mod(ii, bcdata(mm)%icend - bcdata(mm)%icbeg + 1) + bcdata(mm&
&         )%icbeg
        j = ii/(bcdata(mm)%icend-bcdata(mm)%icbeg+1) + bcdata(mm)%jcbeg
        select case  (bcfaceid(mm)) 
        case (imin) 
          mult = -one
          xxp = si(1, i, j, 1)
          yyp = si(1, i, j, 2)
          zzp = si(1, i, j, 3)
        case (imax) 
          mult = one
          xxp = si(il, i, j, 1)
          yyp = si(il, i, j, 2)
          zzp = si(il, i, j, 3)
        case (jmin) 
          mult = -one
          xxp = sj(i, 1, j, 1)
          yyp = sj(i, 1, j, 2)
          zzp = sj(i, 1, j, 3)
        case (jmax) 
          mult = one
          xxp = sj(i, jl, j, 1)
          yyp = sj(i, jl, j, 2)
          zzp = sj(i, jl, j, 3)
        case (kmin) 
          mult = -one
          xxp = sk(i, j, 1, 1)
          yyp = sk(i, j, 1, 2)
          zzp = sk(i, j, 1, 3)
        case (kmax) 
          mult = one
          xxp = sk(i, j, kl, 1)
          yyp = sk(i, j, kl, 2)
          zzp = sk(i, j, kl, 3)
        end select
! compute the inverse of the length of the normal vector
! and possibly correct for inward pointing.
        arg1 = xxp*xxp + yyp*yyp + zzp*zzp
        fact = sqrt(arg1)
        if (fact .gt. zero) fact = mult/fact
! compute the unit normal.
        bcdata(mm)%norm(i, j, 1) = fact*xxp
        bcdata(mm)%norm(i, j, 2) = fact*yyp
        bcdata(mm)%norm(i, j, 3) = fact*zzp
      end do
    end do bocoloop
  end subroutine boundarynormals
 
!  differentiation of xhalo_block in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *x
!   with respect to varying inputs: *x
!   rw status of diff variables: *x:in-out
!   plus diff mem management of: x:in
  subroutine xhalo_block_d()
!
!       xhalo determines the coordinates of the nodal halo's.
!       first it sets all halo coordinates by simple extrapolation,
!       then the symmetry planes are treated (also the unit normal of
!       symmetry planes are determined) and finally an exchange is
!       made for the internal halo's.
!
    use constants
    use blockpointers
    use communication
    use inputtimespectral
    implicit none
!
!      local variables.
!
    integer(kind=inttype) :: mm, i, j, k
    integer(kind=inttype) :: ibeg, iend, jbeg, jend, iimax, jjmax
    logical :: err
    real(kind=realtype) :: length, dot
    real(kind=realtype) :: dotd
    real(kind=realtype), dimension(3) :: v1, v2, norm
    real(kind=realtype), dimension(3) :: v1d
    intrinsic sqrt
    real(kind=realtype) :: arg1
! extrapolation in i-direction.
    do k=1,kl
      do j=1,jl
        xd(0, j, k, 1) = two*xd(1, j, k, 1) - xd(2, j, k, 1)
        x(0, j, k, 1) = two*x(1, j, k, 1) - x(2, j, k, 1)
        xd(0, j, k, 2) = two*xd(1, j, k, 2) - xd(2, j, k, 2)
        x(0, j, k, 2) = two*x(1, j, k, 2) - x(2, j, k, 2)
        xd(0, j, k, 3) = two*xd(1, j, k, 3) - xd(2, j, k, 3)
        x(0, j, k, 3) = two*x(1, j, k, 3) - x(2, j, k, 3)
        xd(ie, j, k, 1) = two*xd(il, j, k, 1) - xd(nx, j, k, 1)
        x(ie, j, k, 1) = two*x(il, j, k, 1) - x(nx, j, k, 1)
        xd(ie, j, k, 2) = two*xd(il, j, k, 2) - xd(nx, j, k, 2)
        x(ie, j, k, 2) = two*x(il, j, k, 2) - x(nx, j, k, 2)
        xd(ie, j, k, 3) = two*xd(il, j, k, 3) - xd(nx, j, k, 3)
        x(ie, j, k, 3) = two*x(il, j, k, 3) - x(nx, j, k, 3)
      end do
    end do
! extrapolation in j-direction.
    do k=1,kl
      do i=0,ie
        xd(i, 0, k, 1) = two*xd(i, 1, k, 1) - xd(i, 2, k, 1)
        x(i, 0, k, 1) = two*x(i, 1, k, 1) - x(i, 2, k, 1)
        xd(i, 0, k, 2) = two*xd(i, 1, k, 2) - xd(i, 2, k, 2)
        x(i, 0, k, 2) = two*x(i, 1, k, 2) - x(i, 2, k, 2)
        xd(i, 0, k, 3) = two*xd(i, 1, k, 3) - xd(i, 2, k, 3)
        x(i, 0, k, 3) = two*x(i, 1, k, 3) - x(i, 2, k, 3)
        xd(i, je, k, 1) = two*xd(i, jl, k, 1) - xd(i, ny, k, 1)
        x(i, je, k, 1) = two*x(i, jl, k, 1) - x(i, ny, k, 1)
        xd(i, je, k, 2) = two*xd(i, jl, k, 2) - xd(i, ny, k, 2)
        x(i, je, k, 2) = two*x(i, jl, k, 2) - x(i, ny, k, 2)
        xd(i, je, k, 3) = two*xd(i, jl, k, 3) - xd(i, ny, k, 3)
        x(i, je, k, 3) = two*x(i, jl, k, 3) - x(i, ny, k, 3)
      end do
    end do
! extrapolation in k-direction.
    do j=0,je
      do i=0,ie
        xd(i, j, 0, 1) = two*xd(i, j, 1, 1) - xd(i, j, 2, 1)
        x(i, j, 0, 1) = two*x(i, j, 1, 1) - x(i, j, 2, 1)
        xd(i, j, 0, 2) = two*xd(i, j, 1, 2) - xd(i, j, 2, 2)
        x(i, j, 0, 2) = two*x(i, j, 1, 2) - x(i, j, 2, 2)
        xd(i, j, 0, 3) = two*xd(i, j, 1, 3) - xd(i, j, 2, 3)
        x(i, j, 0, 3) = two*x(i, j, 1, 3) - x(i, j, 2, 3)
        xd(i, j, ke, 1) = two*xd(i, j, kl, 1) - xd(i, j, nz, 1)
        x(i, j, ke, 1) = two*x(i, j, kl, 1) - x(i, j, nz, 1)
        xd(i, j, ke, 2) = two*xd(i, j, kl, 2) - xd(i, j, nz, 2)
        x(i, j, ke, 2) = two*x(i, j, kl, 2) - x(i, j, nz, 2)
        xd(i, j, ke, 3) = two*xd(i, j, kl, 3) - xd(i, j, nz, 3)
        x(i, j, ke, 3) = two*x(i, j, kl, 3) - x(i, j, nz, 3)
      end do
    end do
    v1d = 0.0_8
!
!           mirror the halo coordinates adjacent to the symmetry
!           planes
!
! loop over boundary subfaces.
loopbocos:do mm=1,nbocos
! the actual correction of the coordinates only takes
! place for symmetry planes.
      if (bctype(mm) .eq. symm) then
! set some variables, depending on the block face on
! which the subface is located.
        norm(1) = bcdata(mm)%symnorm(1)
        norm(2) = bcdata(mm)%symnorm(2)
        norm(3) = bcdata(mm)%symnorm(3)
        arg1 = norm(1)**2 + norm(2)**2 + norm(3)**2
        length = sqrt(arg1)
! compute the unit normal of the subface.
        norm(1) = norm(1)/length
        norm(2) = norm(2)/length
        norm(3) = norm(3)/length
! see xhalo_block for comments for below:
        if (length .gt. eps) then
          select case  (bcfaceid(mm)) 
          case (imin) 
            ibeg = jnbeg(mm)
            iend = jnend(mm)
            iimax = jl
            jbeg = knbeg(mm)
            jend = knend(mm)
            jjmax = kl
            if (ibeg .eq. 1) ibeg = 0
            if (iend .eq. iimax) iend = iimax + 1
            if (jbeg .eq. 1) jbeg = 0
            if (jend .eq. jjmax) jend = jjmax + 1
            do j=jbeg,jend
              do i=ibeg,iend
                v1d(1) = xd(1, i, j, 1) - xd(2, i, j, 1)
                v1(1) = x(1, i, j, 1) - x(2, i, j, 1)
                v1d(2) = xd(1, i, j, 2) - xd(2, i, j, 2)
                v1(2) = x(1, i, j, 2) - x(2, i, j, 2)
                v1d(3) = xd(1, i, j, 3) - xd(2, i, j, 3)
                v1(3) = x(1, i, j, 3) - x(2, i, j, 3)
                dotd = two*(norm(1)*v1d(1)+norm(2)*v1d(2)+norm(3)*v1d(3)&
&                 )
                dot = two*(v1(1)*norm(1)+v1(2)*norm(2)+v1(3)*norm(3))
                xd(0, i, j, 1) = xd(2, i, j, 1) + norm(1)*dotd
                x(0, i, j, 1) = x(2, i, j, 1) + dot*norm(1)
                xd(0, i, j, 2) = xd(2, i, j, 2) + norm(2)*dotd
                x(0, i, j, 2) = x(2, i, j, 2) + dot*norm(2)
                xd(0, i, j, 3) = xd(2, i, j, 3) + norm(3)*dotd
                x(0, i, j, 3) = x(2, i, j, 3) + dot*norm(3)
              end do
            end do
          case (imax) 
            ibeg = jnbeg(mm)
            iend = jnend(mm)
            iimax = jl
            jbeg = knbeg(mm)
            jend = knend(mm)
            jjmax = kl
            if (ibeg .eq. 1) ibeg = 0
            if (iend .eq. iimax) iend = iimax + 1
            if (jbeg .eq. 1) jbeg = 0
            if (jend .eq. jjmax) jend = jjmax + 1
            do j=jbeg,jend
              do i=ibeg,iend
                v1d(1) = xd(il, i, j, 1) - xd(nx, i, j, 1)
                v1(1) = x(il, i, j, 1) - x(nx, i, j, 1)
                v1d(2) = xd(il, i, j, 2) - xd(nx, i, j, 2)
                v1(2) = x(il, i, j, 2) - x(nx, i, j, 2)
                v1d(3) = xd(il, i, j, 3) - xd(nx, i, j, 3)
                v1(3) = x(il, i, j, 3) - x(nx, i, j, 3)
                dotd = two*(norm(1)*v1d(1)+norm(2)*v1d(2)+norm(3)*v1d(3)&
&                 )
                dot = two*(v1(1)*norm(1)+v1(2)*norm(2)+v1(3)*norm(3))
                xd(ie, i, j, 1) = xd(nx, i, j, 1) + norm(1)*dotd
                x(ie, i, j, 1) = x(nx, i, j, 1) + dot*norm(1)
                xd(ie, i, j, 2) = xd(nx, i, j, 2) + norm(2)*dotd
                x(ie, i, j, 2) = x(nx, i, j, 2) + dot*norm(2)
                xd(ie, i, j, 3) = xd(nx, i, j, 3) + norm(3)*dotd
                x(ie, i, j, 3) = x(nx, i, j, 3) + dot*norm(3)
              end do
            end do
          case (jmin) 
            ibeg = inbeg(mm)
            iend = inend(mm)
            iimax = il
            jbeg = knbeg(mm)
            jend = knend(mm)
            jjmax = kl
            if (ibeg .eq. 1) ibeg = 0
            if (iend .eq. iimax) iend = iimax + 1
            if (jbeg .eq. 1) jbeg = 0
            if (jend .eq. jjmax) jend = jjmax + 1
            do j=jbeg,jend
              do i=ibeg,iend
                v1d(1) = xd(i, 1, j, 1) - xd(i, 2, j, 1)
                v1(1) = x(i, 1, j, 1) - x(i, 2, j, 1)
                v1d(2) = xd(i, 1, j, 2) - xd(i, 2, j, 2)
                v1(2) = x(i, 1, j, 2) - x(i, 2, j, 2)
                v1d(3) = xd(i, 1, j, 3) - xd(i, 2, j, 3)
                v1(3) = x(i, 1, j, 3) - x(i, 2, j, 3)
                dotd = two*(norm(1)*v1d(1)+norm(2)*v1d(2)+norm(3)*v1d(3)&
&                 )
                dot = two*(v1(1)*norm(1)+v1(2)*norm(2)+v1(3)*norm(3))
                xd(i, 0, j, 1) = xd(i, 2, j, 1) + norm(1)*dotd
                x(i, 0, j, 1) = x(i, 2, j, 1) + dot*norm(1)
                xd(i, 0, j, 2) = xd(i, 2, j, 2) + norm(2)*dotd
                x(i, 0, j, 2) = x(i, 2, j, 2) + dot*norm(2)
                xd(i, 0, j, 3) = xd(i, 2, j, 3) + norm(3)*dotd
                x(i, 0, j, 3) = x(i, 2, j, 3) + dot*norm(3)
              end do
            end do
          case (jmax) 
            ibeg = inbeg(mm)
            iend = inend(mm)
            iimax = il
            jbeg = knbeg(mm)
            jend = knend(mm)
            jjmax = kl
            if (ibeg .eq. 1) ibeg = 0
            if (iend .eq. iimax) iend = iimax + 1
            if (jbeg .eq. 1) jbeg = 0
            if (jend .eq. jjmax) jend = jjmax + 1
            do j=jbeg,jend
              do i=ibeg,iend
                v1d(1) = xd(i, jl, j, 1) - xd(i, ny, j, 1)
                v1(1) = x(i, jl, j, 1) - x(i, ny, j, 1)
                v1d(2) = xd(i, jl, j, 2) - xd(i, ny, j, 2)
                v1(2) = x(i, jl, j, 2) - x(i, ny, j, 2)
                v1d(3) = xd(i, jl, j, 3) - xd(i, ny, j, 3)
                v1(3) = x(i, jl, j, 3) - x(i, ny, j, 3)
                dotd = two*(norm(1)*v1d(1)+norm(2)*v1d(2)+norm(3)*v1d(3)&
&                 )
                dot = two*(v1(1)*norm(1)+v1(2)*norm(2)+v1(3)*norm(3))
                xd(i, je, j, 1) = xd(i, ny, j, 1) + norm(1)*dotd
                x(i, je, j, 1) = x(i, ny, j, 1) + dot*norm(1)
                xd(i, je, j, 2) = xd(i, ny, j, 2) + norm(2)*dotd
                x(i, je, j, 2) = x(i, ny, j, 2) + dot*norm(2)
                xd(i, je, j, 3) = xd(i, ny, j, 3) + norm(3)*dotd
                x(i, je, j, 3) = x(i, ny, j, 3) + dot*norm(3)
              end do
            end do
          case (kmin) 
            ibeg = inbeg(mm)
            iend = inend(mm)
            iimax = il
            jbeg = jnbeg(mm)
            jend = jnend(mm)
            jjmax = jl
            if (ibeg .eq. 1) ibeg = 0
            if (iend .eq. iimax) iend = iimax + 1
            if (jbeg .eq. 1) jbeg = 0
            if (jend .eq. jjmax) jend = jjmax + 1
            do j=jbeg,jend
              do i=ibeg,iend
                v1d(1) = xd(i, j, 1, 1) - xd(i, j, 2, 1)
                v1(1) = x(i, j, 1, 1) - x(i, j, 2, 1)
                v1d(2) = xd(i, j, 1, 2) - xd(i, j, 2, 2)
                v1(2) = x(i, j, 1, 2) - x(i, j, 2, 2)
                v1d(3) = xd(i, j, 1, 3) - xd(i, j, 2, 3)
                v1(3) = x(i, j, 1, 3) - x(i, j, 2, 3)
                dotd = two*(norm(1)*v1d(1)+norm(2)*v1d(2)+norm(3)*v1d(3)&
&                 )
                dot = two*(v1(1)*norm(1)+v1(2)*norm(2)+v1(3)*norm(3))
                xd(i, j, 0, 1) = xd(i, j, 2, 1) + norm(1)*dotd
                x(i, j, 0, 1) = x(i, j, 2, 1) + dot*norm(1)
                xd(i, j, 0, 2) = xd(i, j, 2, 2) + norm(2)*dotd
                x(i, j, 0, 2) = x(i, j, 2, 2) + dot*norm(2)
                xd(i, j, 0, 3) = xd(i, j, 2, 3) + norm(3)*dotd
                x(i, j, 0, 3) = x(i, j, 2, 3) + dot*norm(3)
              end do
            end do
          case (kmax) 
            ibeg = inbeg(mm)
            iend = inend(mm)
            iimax = il
            jbeg = jnbeg(mm)
            jend = jnend(mm)
            jjmax = jl
            if (ibeg .eq. 1) ibeg = 0
            if (iend .eq. iimax) iend = iimax + 1
            if (jbeg .eq. 1) jbeg = 0
            if (jend .eq. jjmax) jend = jjmax + 1
            do j=jbeg,jend
              do i=ibeg,iend
                v1d(1) = xd(i, j, kl, 1) - xd(i, j, nz, 1)
                v1(1) = x(i, j, kl, 1) - x(i, j, nz, 1)
                v1d(2) = xd(i, j, kl, 2) - xd(i, j, nz, 2)
                v1(2) = x(i, j, kl, 2) - x(i, j, nz, 2)
                v1d(3) = xd(i, j, kl, 3) - xd(i, j, nz, 3)
                v1(3) = x(i, j, kl, 3) - x(i, j, nz, 3)
                dotd = two*(norm(1)*v1d(1)+norm(2)*v1d(2)+norm(3)*v1d(3)&
&                 )
                dot = two*(v1(1)*norm(1)+v1(2)*norm(2)+v1(3)*norm(3))
                xd(i, j, ke, 1) = xd(i, j, nz, 1) + norm(1)*dotd
                x(i, j, ke, 1) = x(i, j, nz, 1) + dot*norm(1)
                xd(i, j, ke, 2) = xd(i, j, nz, 2) + norm(2)*dotd
                x(i, j, ke, 2) = x(i, j, nz, 2) + dot*norm(2)
                xd(i, j, ke, 3) = xd(i, j, nz, 3) + norm(3)*dotd
                x(i, j, ke, 3) = x(i, j, nz, 3) + dot*norm(3)
              end do
            end do
          end select
        end if
      end if
    end do loopbocos
  end subroutine xhalo_block_d
 
  subroutine xhalo_block()
!
!       xhalo determines the coordinates of the nodal halo's.
!       first it sets all halo coordinates by simple extrapolation,
!       then the symmetry planes are treated (also the unit normal of
!       symmetry planes are determined) and finally an exchange is
!       made for the internal halo's.
!
    use constants
    use blockpointers
    use communication
    use inputtimespectral
    implicit none
!
!      local variables.
!
    integer(kind=inttype) :: mm, i, j, k
    integer(kind=inttype) :: ibeg, iend, jbeg, jend, iimax, jjmax
    logical :: err
    real(kind=realtype) :: length, dot
    real(kind=realtype), dimension(3) :: v1, v2, norm
    intrinsic sqrt
    real(kind=realtype) :: arg1
! extrapolation in i-direction.
    do k=1,kl
      do j=1,jl
        x(0, j, k, 1) = two*x(1, j, k, 1) - x(2, j, k, 1)
        x(0, j, k, 2) = two*x(1, j, k, 2) - x(2, j, k, 2)
        x(0, j, k, 3) = two*x(1, j, k, 3) - x(2, j, k, 3)
        x(ie, j, k, 1) = two*x(il, j, k, 1) - x(nx, j, k, 1)
        x(ie, j, k, 2) = two*x(il, j, k, 2) - x(nx, j, k, 2)
        x(ie, j, k, 3) = two*x(il, j, k, 3) - x(nx, j, k, 3)
      end do
    end do
! extrapolation in j-direction.
    do k=1,kl
      do i=0,ie
        x(i, 0, k, 1) = two*x(i, 1, k, 1) - x(i, 2, k, 1)
        x(i, 0, k, 2) = two*x(i, 1, k, 2) - x(i, 2, k, 2)
        x(i, 0, k, 3) = two*x(i, 1, k, 3) - x(i, 2, k, 3)
        x(i, je, k, 1) = two*x(i, jl, k, 1) - x(i, ny, k, 1)
        x(i, je, k, 2) = two*x(i, jl, k, 2) - x(i, ny, k, 2)
        x(i, je, k, 3) = two*x(i, jl, k, 3) - x(i, ny, k, 3)
      end do
    end do
! extrapolation in k-direction.
    do j=0,je
      do i=0,ie
        x(i, j, 0, 1) = two*x(i, j, 1, 1) - x(i, j, 2, 1)
        x(i, j, 0, 2) = two*x(i, j, 1, 2) - x(i, j, 2, 2)
        x(i, j, 0, 3) = two*x(i, j, 1, 3) - x(i, j, 2, 3)
        x(i, j, ke, 1) = two*x(i, j, kl, 1) - x(i, j, nz, 1)
        x(i, j, ke, 2) = two*x(i, j, kl, 2) - x(i, j, nz, 2)
        x(i, j, ke, 3) = two*x(i, j, kl, 3) - x(i, j, nz, 3)
      end do
    end do
!
!           mirror the halo coordinates adjacent to the symmetry
!           planes
!
! loop over boundary subfaces.
loopbocos:do mm=1,nbocos
! the actual correction of the coordinates only takes
! place for symmetry planes.
      if (bctype(mm) .eq. symm) then
! set some variables, depending on the block face on
! which the subface is located.
        norm(1) = bcdata(mm)%symnorm(1)
        norm(2) = bcdata(mm)%symnorm(2)
        norm(3) = bcdata(mm)%symnorm(3)
        arg1 = norm(1)**2 + norm(2)**2 + norm(3)**2
        length = sqrt(arg1)
! compute the unit normal of the subface.
        norm(1) = norm(1)/length
        norm(2) = norm(2)/length
        norm(3) = norm(3)/length
! see xhalo_block for comments for below:
        if (length .gt. eps) then
          select case  (bcfaceid(mm)) 
          case (imin) 
            ibeg = jnbeg(mm)
            iend = jnend(mm)
            iimax = jl
            jbeg = knbeg(mm)
            jend = knend(mm)
            jjmax = kl
            if (ibeg .eq. 1) ibeg = 0
            if (iend .eq. iimax) iend = iimax + 1
            if (jbeg .eq. 1) jbeg = 0
            if (jend .eq. jjmax) jend = jjmax + 1
            do j=jbeg,jend
              do i=ibeg,iend
                v1(1) = x(1, i, j, 1) - x(2, i, j, 1)
                v1(2) = x(1, i, j, 2) - x(2, i, j, 2)
                v1(3) = x(1, i, j, 3) - x(2, i, j, 3)
                dot = two*(v1(1)*norm(1)+v1(2)*norm(2)+v1(3)*norm(3))
                x(0, i, j, 1) = x(2, i, j, 1) + dot*norm(1)
                x(0, i, j, 2) = x(2, i, j, 2) + dot*norm(2)
                x(0, i, j, 3) = x(2, i, j, 3) + dot*norm(3)
              end do
            end do
          case (imax) 
            ibeg = jnbeg(mm)
            iend = jnend(mm)
            iimax = jl
            jbeg = knbeg(mm)
            jend = knend(mm)
            jjmax = kl
            if (ibeg .eq. 1) ibeg = 0
            if (iend .eq. iimax) iend = iimax + 1
            if (jbeg .eq. 1) jbeg = 0
            if (jend .eq. jjmax) jend = jjmax + 1
            do j=jbeg,jend
              do i=ibeg,iend
                v1(1) = x(il, i, j, 1) - x(nx, i, j, 1)
                v1(2) = x(il, i, j, 2) - x(nx, i, j, 2)
                v1(3) = x(il, i, j, 3) - x(nx, i, j, 3)
                dot = two*(v1(1)*norm(1)+v1(2)*norm(2)+v1(3)*norm(3))
                x(ie, i, j, 1) = x(nx, i, j, 1) + dot*norm(1)
                x(ie, i, j, 2) = x(nx, i, j, 2) + dot*norm(2)
                x(ie, i, j, 3) = x(nx, i, j, 3) + dot*norm(3)
              end do
            end do
          case (jmin) 
            ibeg = inbeg(mm)
            iend = inend(mm)
            iimax = il
            jbeg = knbeg(mm)
            jend = knend(mm)
            jjmax = kl
            if (ibeg .eq. 1) ibeg = 0
            if (iend .eq. iimax) iend = iimax + 1
            if (jbeg .eq. 1) jbeg = 0
            if (jend .eq. jjmax) jend = jjmax + 1
            do j=jbeg,jend
              do i=ibeg,iend
                v1(1) = x(i, 1, j, 1) - x(i, 2, j, 1)
                v1(2) = x(i, 1, j, 2) - x(i, 2, j, 2)
                v1(3) = x(i, 1, j, 3) - x(i, 2, j, 3)
                dot = two*(v1(1)*norm(1)+v1(2)*norm(2)+v1(3)*norm(3))
                x(i, 0, j, 1) = x(i, 2, j, 1) + dot*norm(1)
                x(i, 0, j, 2) = x(i, 2, j, 2) + dot*norm(2)
                x(i, 0, j, 3) = x(i, 2, j, 3) + dot*norm(3)
              end do
            end do
          case (jmax) 
            ibeg = inbeg(mm)
            iend = inend(mm)
            iimax = il
            jbeg = knbeg(mm)
            jend = knend(mm)
            jjmax = kl
            if (ibeg .eq. 1) ibeg = 0
            if (iend .eq. iimax) iend = iimax + 1
            if (jbeg .eq. 1) jbeg = 0
            if (jend .eq. jjmax) jend = jjmax + 1
            do j=jbeg,jend
              do i=ibeg,iend
                v1(1) = x(i, jl, j, 1) - x(i, ny, j, 1)
                v1(2) = x(i, jl, j, 2) - x(i, ny, j, 2)
                v1(3) = x(i, jl, j, 3) - x(i, ny, j, 3)
                dot = two*(v1(1)*norm(1)+v1(2)*norm(2)+v1(3)*norm(3))
                x(i, je, j, 1) = x(i, ny, j, 1) + dot*norm(1)
                x(i, je, j, 2) = x(i, ny, j, 2) + dot*norm(2)
                x(i, je, j, 3) = x(i, ny, j, 3) + dot*norm(3)
              end do
            end do
          case (kmin) 
            ibeg = inbeg(mm)
            iend = inend(mm)
            iimax = il
            jbeg = jnbeg(mm)
            jend = jnend(mm)
            jjmax = jl
            if (ibeg .eq. 1) ibeg = 0
            if (iend .eq. iimax) iend = iimax + 1
            if (jbeg .eq. 1) jbeg = 0
            if (jend .eq. jjmax) jend = jjmax + 1
            do j=jbeg,jend
              do i=ibeg,iend
                v1(1) = x(i, j, 1, 1) - x(i, j, 2, 1)
                v1(2) = x(i, j, 1, 2) - x(i, j, 2, 2)
                v1(3) = x(i, j, 1, 3) - x(i, j, 2, 3)
                dot = two*(v1(1)*norm(1)+v1(2)*norm(2)+v1(3)*norm(3))
                x(i, j, 0, 1) = x(i, j, 2, 1) + dot*norm(1)
                x(i, j, 0, 2) = x(i, j, 2, 2) + dot*norm(2)
                x(i, j, 0, 3) = x(i, j, 2, 3) + dot*norm(3)
              end do
            end do
          case (kmax) 
            ibeg = inbeg(mm)
            iend = inend(mm)
            iimax = il
            jbeg = jnbeg(mm)
            jend = jnend(mm)
            jjmax = jl
            if (ibeg .eq. 1) ibeg = 0
            if (iend .eq. iimax) iend = iimax + 1
            if (jbeg .eq. 1) jbeg = 0
            if (jend .eq. jjmax) jend = jjmax + 1
            do j=jbeg,jend
              do i=ibeg,iend
                v1(1) = x(i, j, kl, 1) - x(i, j, nz, 1)
                v1(2) = x(i, j, kl, 2) - x(i, j, nz, 2)
                v1(3) = x(i, j, kl, 3) - x(i, j, nz, 3)
                dot = two*(v1(1)*norm(1)+v1(2)*norm(2)+v1(3)*norm(3))
                x(i, j, ke, 1) = x(i, j, nz, 1) + dot*norm(1)
                x(i, j, ke, 2) = x(i, j, nz, 2) + dot*norm(2)
                x(i, j, ke, 3) = x(i, j, nz, 3) + dot*norm(3)
              end do
            end do
          end select
        end if
      end if
    end do loopbocos
  end subroutine xhalo_block
 
!  differentiation of resscale in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *dw
!   with respect to varying inputs: *dw
!   rw status of diff variables: *dw:in-out
!   plus diff mem management of: dw:in
  subroutine resscale_d()
    use constants
    use blockpointers, only : il, jl, kl, nx, ny, nz, volref, dw, dwd
    use flowvarrefstate, only : nwf, nt1, nt2
    use inputiteration, only : turbresscale
    implicit none
! local variables
    integer(kind=inttype) :: i, j, k, ii, nturb
    real(kind=realtype) :: ovol
! divide through by the reference volume
    nturb = nt2 - nt1 + 1
    do k=2,kl
      do j=2,jl
        do i=2,il
          ovol = one/volref(i, j, k)
          dwd(i, j, k, 1:nwf) = ovol*dwd(i, j, k, 1:nwf)
          dw(i, j, k, 1:nwf) = dw(i, j, k, 1:nwf)*ovol
          dwd(i, j, k, nt1:nt2) = ovol*turbresscale(1:nturb)*dwd(i, j, k&
&           , nt1:nt2)
          dw(i, j, k, nt1:nt2) = dw(i, j, k, nt1:nt2)*ovol*turbresscale(&
&           1:nturb)
        end do
      end do
    end do
  end subroutine resscale_d
 
  subroutine resscale()
    use constants
    use blockpointers, only : il, jl, kl, nx, ny, nz, volref, dw
    use flowvarrefstate, only : nwf, nt1, nt2
    use inputiteration, only : turbresscale
    implicit none
! local variables
    integer(kind=inttype) :: i, j, k, ii, nturb
    real(kind=realtype) :: ovol
! divide through by the reference volume
    nturb = nt2 - nt1 + 1
    do k=2,kl
      do j=2,jl
        do i=2,il
          ovol = one/volref(i, j, k)
          dw(i, j, k, 1:nwf) = dw(i, j, k, 1:nwf)*ovol
          dw(i, j, k, nt1:nt2) = dw(i, j, k, nt1:nt2)*ovol*turbresscale(&
&           1:nturb)
        end do
      end do
    end do
  end subroutine resscale
 
!  differentiation of sumdwandfw in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *dw
!   with respect to varying inputs: *dw *fw
!   rw status of diff variables: *dw:in-out *fw:in
!   plus diff mem management of: dw:in fw:in
  subroutine sumdwandfw_d()
    use constants
    use blockpointers, only : il, jl, kl, dw, dwd, fw, fwd, iblank
    use flowvarrefstate, only : nwf
    implicit none
! local variables
    integer(kind=inttype) :: i, j, k, l
    intrinsic real
    intrinsic max
    real(realtype) :: x1
    real(realtype) :: max1
    do l=1,nwf
      do k=2,kl
        do j=2,jl
          do i=2,il
            x1 = real(iblank(i, j, k), realtype)
            if (x1 .lt. zero) then
              max1 = zero
            else
              max1 = x1
            end if
            dwd(i, j, k, l) = max1*(dwd(i, j, k, l)+fwd(i, j, k, l))
            dw(i, j, k, l) = (dw(i, j, k, l)+fw(i, j, k, l))*max1
          end do
        end do
      end do
    end do
  end subroutine sumdwandfw_d
 
  subroutine sumdwandfw()
    use constants
    use blockpointers, only : il, jl, kl, dw, fw, iblank
    use flowvarrefstate, only : nwf
    implicit none
! local variables
    integer(kind=inttype) :: i, j, k, l
    intrinsic real
    intrinsic max
    real(realtype) :: x1
    real(realtype) :: max1
    do l=1,nwf
      do k=2,kl
        do j=2,jl
          do i=2,il
            x1 = real(iblank(i, j, k), realtype)
            if (x1 .lt. zero) then
              max1 = zero
            else
              max1 = x1
            end if
            dw(i, j, k, l) = (dw(i, j, k, l)+fw(i, j, k, l))*max1
          end do
        end do
      end do
    end do
  end subroutine sumdwandfw
 
end module adjointextra_d