adjointExtra_b.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_b
  implicit none
 
contains
!  differentiation of volume_block in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *x *vol
!   with respect to varying inputs: *x *vol
!   rw status of diff variables: *x:incr *vol:in-out
!   plus diff mem management of: x:in vol:in
  subroutine volume_block_b()
! 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
    real(kind=realtype) :: tempd
    real(kind=realtype) :: tmp
    real(kind=realtype) :: tmpd
    real(kind=realtype) :: tmp0
    real(kind=realtype) :: tmpd0
    real(kind=realtype) :: tmp1
    real(kind=realtype) :: tmpd1
    integer :: branch
! 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
      call pushinteger4(n)
      n = k - 1
      do j=1,je
        call pushinteger4(m)
        m = j - 1
        do i=1,ie
          l = i - 1
! compute the coordinates of the center of gravity.
          call pushreal8(xp)
          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))
          call pushreal8(yp)
          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))
          call pushreal8(zp)
          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
            call pushcontrol1b(0)
            vol(i, j, k) = vol(i, j, k)
          else
            vol(i, j, k) = -vol(i, j, k)
            call pushcontrol1b(1)
          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
          vol(1, j, k) = vol(2, j, k)
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
        end if
        if (vol(ie, j, k)/vol(il, j, k) .lt. halocellratio) then
          tmp = vol(il, j, k)
          vol(ie, j, k) = tmp
          call pushcontrol1b(1)
        else
          call pushcontrol1b(0)
        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
          vol(i, 1, k) = vol(i, 2, k)
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
        end if
        if (vol(i, je, k)/vol(i, jl, k) .lt. halocellratio) then
          tmp0 = vol(i, jl, k)
          vol(i, je, k) = tmp0
          call pushcontrol1b(1)
        else
          call pushcontrol1b(0)
        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
          vol(i, j, 1) = vol(i, j, 2)
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
        end if
        if (vol(i, j, ke)/vol(i, j, kl) .lt. halocellratio) then
          tmp1 = vol(i, j, kl)
          vol(i, j, ke) = tmp1
          call pushcontrol1b(1)
        else
          call pushcontrol1b(0)
        end if
      end do
    end do
    do j=je,1,-1
      do i=ie,1,-1
        call popcontrol1b(branch)
        if (branch .ne. 0) then
          tmpd1 = vold(i, j, ke)
          vold(i, j, ke) = 0.0_8
          vold(i, j, kl) = vold(i, j, kl) + tmpd1
        end if
        call popcontrol1b(branch)
        if (branch .eq. 0) then
          vold(i, j, 2) = vold(i, j, 2) + vold(i, j, 1)
          vold(i, j, 1) = 0.0_8
        end if
      end do
    end do
    do k=kl,2,-1
      do i=ie,1,-1
        call popcontrol1b(branch)
        if (branch .ne. 0) then
          tmpd0 = vold(i, je, k)
          vold(i, je, k) = 0.0_8
          vold(i, jl, k) = vold(i, jl, k) + tmpd0
        end if
        call popcontrol1b(branch)
        if (branch .eq. 0) then
          vold(i, 2, k) = vold(i, 2, k) + vold(i, 1, k)
          vold(i, 1, k) = 0.0_8
        end if
      end do
    end do
    do k=kl,2,-1
      do j=jl,2,-1
        call popcontrol1b(branch)
        if (branch .ne. 0) then
          tmpd = vold(ie, j, k)
          vold(ie, j, k) = 0.0_8
          vold(il, j, k) = vold(il, j, k) + tmpd
        end if
        call popcontrol1b(branch)
        if (branch .eq. 0) then
          vold(2, j, k) = vold(2, j, k) + vold(1, j, k)
          vold(1, j, k) = 0.0_8
        end if
      end do
    end do
    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=ke,1,-1
      do j=je,1,-1
        do i=ie,1,-1
          call popcontrol1b(branch)
          if (branch .ne. 0) vold(i, j, k) = -vold(i, j, k)
          tempd = sixth*vold(i, j, k)
          vold(i, j, k) = 0.0_8
          vp1d = vp1d + tempd
          vp2d = vp2d + tempd
          vp3d = vp3d + tempd
          vp4d = vp4d + tempd
          vp5d = vp5d + tempd
          vp6d = vp6d + tempd
          l = i - 1
          zpd = 0.0_8
          ypd = 0.0_8
          xpd = 0.0_8
          call volpym_b(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)
          vp6d = 0.0_8
          call volpym_b(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)
          vp5d = 0.0_8
          call volpym_b(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)
          vp4d = 0.0_8
          call volpym_b(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)
          vp3d = 0.0_8
          call volpym_b(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)
          vp2d = 0.0_8
          call volpym_b(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)
          vp1d = 0.0_8
          call popreal8(zp)
          tempd = eighth*zpd
          xd(i, j, k, 3) = xd(i, j, k, 3) + tempd
          xd(i, m, k, 3) = xd(i, m, k, 3) + tempd
          xd(i, m, n, 3) = xd(i, m, n, 3) + tempd
          xd(i, j, n, 3) = xd(i, j, n, 3) + tempd
          xd(l, j, k, 3) = xd(l, j, k, 3) + tempd
          xd(l, m, k, 3) = xd(l, m, k, 3) + tempd
          xd(l, m, n, 3) = xd(l, m, n, 3) + tempd
          xd(l, j, n, 3) = xd(l, j, n, 3) + tempd
          call popreal8(yp)
          tempd = eighth*ypd
          xd(i, j, k, 2) = xd(i, j, k, 2) + tempd
          xd(i, m, k, 2) = xd(i, m, k, 2) + tempd
          xd(i, m, n, 2) = xd(i, m, n, 2) + tempd
          xd(i, j, n, 2) = xd(i, j, n, 2) + tempd
          xd(l, j, k, 2) = xd(l, j, k, 2) + tempd
          xd(l, m, k, 2) = xd(l, m, k, 2) + tempd
          xd(l, m, n, 2) = xd(l, m, n, 2) + tempd
          xd(l, j, n, 2) = xd(l, j, n, 2) + tempd
          call popreal8(xp)
          tempd = eighth*xpd
          xd(i, j, k, 1) = xd(i, j, k, 1) + tempd
          xd(i, m, k, 1) = xd(i, m, k, 1) + tempd
          xd(i, m, n, 1) = xd(i, m, n, 1) + tempd
          xd(i, j, n, 1) = xd(i, j, n, 1) + tempd
          xd(l, j, k, 1) = xd(l, j, k, 1) + tempd
          xd(l, m, k, 1) = xd(l, m, k, 1) + tempd
          xd(l, m, n, 1) = xd(l, m, n, 1) + tempd
          xd(l, j, n, 1) = xd(l, j, n, 1) + tempd
        end do
        call popinteger4(m)
      end do
      call popinteger4(n)
    end do
    vold = 0.0_8
 
  contains
!  differentiation of volpym in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: xp yp zp xa xb xc xd ya yb
!                yc yd za zb zc zd volume
!   with respect to varying inputs: xp yp zp xa xb xc xd ya yb
!                yc yd za zb zc zd
    subroutine volpym_b(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) :: xad, yad, zad, xbd, ybd, zbd
      real(kind=realtype), intent(in) :: xc, yc, zc, xd, yd, zd
      real(kind=realtype) :: xcd, ycd, zcd, xdd, ydd, zdd
      real(kind=realtype) :: tempd
      real(kind=realtype) :: tempd0
      real(kind=realtype) :: tempd1
      real(kind=realtype) :: tempd2
      real(kind=realtype) :: tempd3
      real(kind=realtype) :: tempd4
      real(kind=realtype) :: tempd5
      real(kind=realtype) :: tempd6
      real(kind=realtype) :: tempd7
      tempd = ((ya-yc)*(zb-zd)-(za-zc)*(yb-yd))*volumed
      tempd1 = (xp-fourth*(xa+xb+xc+xd))*volumed
      tempd2 = ((za-zc)*(xb-xd)-(xa-xc)*(zb-zd))*volumed
      tempd4 = (yp-fourth*(ya+yb+yc+yd))*volumed
      tempd5 = ((xa-xc)*(yb-yd)-(ya-yc)*(xb-xd))*volumed
      tempd7 = (zp-fourth*(za+zb+zc+zd))*volumed
      zpd = zpd + tempd5
      tempd6 = -(fourth*tempd5)
      zad = zad + tempd6 + (xb-xd)*tempd4 - (yb-yd)*tempd1
      zbd = zbd + tempd6 + (ya-yc)*tempd1 - (xa-xc)*tempd4
      zcd = zcd + tempd6 + (yb-yd)*tempd1 - (xb-xd)*tempd4
      zdd = zdd + tempd6 + (xa-xc)*tempd4 - (ya-yc)*tempd1
      ypd = ypd + tempd2
      tempd3 = -(fourth*tempd2)
      ybd = ybd + (xa-xc)*tempd7 + tempd3 - (za-zc)*tempd1
      ydd = ydd + tempd3 - (xa-xc)*tempd7 + (za-zc)*tempd1
      yad = yad + tempd3 - (xb-xd)*tempd7 + (zb-zd)*tempd1
      ycd = ycd + (xb-xd)*tempd7 + tempd3 - (zb-zd)*tempd1
      xpd = xpd + tempd
      tempd0 = -(fourth*tempd)
      xad = xad + (yb-yd)*tempd7 + tempd0 - (zb-zd)*tempd4
      xcd = xcd + (zb-zd)*tempd4 - (yb-yd)*tempd7 + tempd0
      xbd = xbd + (za-zc)*tempd4 - (ya-yc)*tempd7 + tempd0
      xdd = xdd + (ya-yc)*tempd7 + tempd0 - (za-zc)*tempd4
    end subroutine volpym_b
 
    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_b
 
  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 reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *x *si *sj *sk
!   with respect to varying inputs: *x *si *sj *sk
!   rw status of diff variables: *x:incr *si:in-out *sj:in-out
!                *sk:in-out
!   plus diff mem management of: x:in si:in sj:in sk:in
  subroutine metric_block_b()
    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
    real(kind=realtype) :: tempd
! 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
    call pushreal8array(v1, 3)
    call pushreal8array(v2, 3)
!$fwd-of 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
    call pushreal8array(v1, 3)
    call pushreal8array(v2, 3)
!$fwd-of 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
    v1d = 0.0_8
    v2d = 0.0_8
!$bwd-of ii-loop 
    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.
      tempd = fact*skd(i, j, k, 3)
      skd(i, j, k, 3) = 0.0_8
      v1d(1) = v1d(1) + v2(2)*tempd
      v2d(2) = v2d(2) + v1(1)*tempd
      v1d(2) = v1d(2) - v2(1)*tempd
      v2d(1) = v2d(1) - v1(2)*tempd
      tempd = fact*skd(i, j, k, 2)
      skd(i, j, k, 2) = 0.0_8
      v1d(3) = v1d(3) + v2(1)*tempd
      v2d(1) = v2d(1) + v1(3)*tempd
      v1d(1) = v1d(1) - v2(3)*tempd
      v2d(3) = v2d(3) - v1(1)*tempd
      tempd = fact*skd(i, j, k, 1)
      skd(i, j, k, 1) = 0.0_8
      v1d(2) = v1d(2) + v2(3)*tempd
      v2d(3) = v2d(3) + v1(2)*tempd
      v1d(3) = v1d(3) - v2(2)*tempd
      v2d(2) = v2d(2) - v1(3)*tempd
      xd(l, j, k, 3) = xd(l, j, k, 3) + v2d(3)
      xd(i, m, k, 3) = xd(i, m, k, 3) - v2d(3)
      v2d(3) = 0.0_8
      xd(l, j, k, 2) = xd(l, j, k, 2) + v2d(2)
      xd(i, m, k, 2) = xd(i, m, k, 2) - v2d(2)
      v2d(2) = 0.0_8
      xd(l, j, k, 1) = xd(l, j, k, 1) + v2d(1)
      xd(i, m, k, 1) = xd(i, m, k, 1) - v2d(1)
      v2d(1) = 0.0_8
      xd(i, j, k, 3) = xd(i, j, k, 3) + v1d(3)
      xd(l, m, k, 3) = xd(l, m, k, 3) - v1d(3)
      v1d(3) = 0.0_8
      xd(i, j, k, 2) = xd(i, j, k, 2) + v1d(2)
      xd(l, m, k, 2) = xd(l, m, k, 2) - v1d(2)
      v1d(2) = 0.0_8
      xd(i, j, k, 1) = xd(i, j, k, 1) + v1d(1)
      xd(l, m, k, 1) = xd(l, m, k, 1) - v1d(1)
      v1d(1) = 0.0_8
    end do
    call popreal8array(v2, 3)
    call popreal8array(v1, 3)
!$bwd-of ii-loop 
    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.
      tempd = fact*sjd(i, j, k, 3)
      sjd(i, j, k, 3) = 0.0_8
      v1d(1) = v1d(1) + v2(2)*tempd
      v2d(2) = v2d(2) + v1(1)*tempd
      v1d(2) = v1d(2) - v2(1)*tempd
      v2d(1) = v2d(1) - v1(2)*tempd
      tempd = fact*sjd(i, j, k, 2)
      sjd(i, j, k, 2) = 0.0_8
      v1d(3) = v1d(3) + v2(1)*tempd
      v2d(1) = v2d(1) + v1(3)*tempd
      v1d(1) = v1d(1) - v2(3)*tempd
      v2d(3) = v2d(3) - v1(1)*tempd
      tempd = fact*sjd(i, j, k, 1)
      sjd(i, j, k, 1) = 0.0_8
      v1d(2) = v1d(2) + v2(3)*tempd
      v2d(3) = v2d(3) + v1(2)*tempd
      v1d(3) = v1d(3) - v2(2)*tempd
      v2d(2) = v2d(2) - v1(3)*tempd
      xd(l, j, n, 3) = xd(l, j, n, 3) + v2d(3)
      xd(i, j, k, 3) = xd(i, j, k, 3) - v2d(3)
      v2d(3) = 0.0_8
      xd(l, j, n, 2) = xd(l, j, n, 2) + v2d(2)
      xd(i, j, k, 2) = xd(i, j, k, 2) - v2d(2)
      v2d(2) = 0.0_8
      xd(l, j, n, 1) = xd(l, j, n, 1) + v2d(1)
      xd(i, j, k, 1) = xd(i, j, k, 1) - v2d(1)
      v2d(1) = 0.0_8
      xd(i, j, n, 3) = xd(i, j, n, 3) + v1d(3)
      xd(l, j, k, 3) = xd(l, j, k, 3) - v1d(3)
      v1d(3) = 0.0_8
      xd(i, j, n, 2) = xd(i, j, n, 2) + v1d(2)
      xd(l, j, k, 2) = xd(l, j, k, 2) - v1d(2)
      v1d(2) = 0.0_8
      xd(i, j, n, 1) = xd(i, j, n, 1) + v1d(1)
      xd(l, j, k, 1) = xd(l, j, k, 1) - v1d(1)
      v1d(1) = 0.0_8
    end do
    call popreal8array(v2, 3)
    call popreal8array(v1, 3)
!$bwd-of ii-loop 
    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.
      tempd = fact*sid(i, j, k, 3)
      sid(i, j, k, 3) = 0.0_8
      v1d(1) = v1d(1) + v2(2)*tempd
      v2d(2) = v2d(2) + v1(1)*tempd
      v1d(2) = v1d(2) - v2(1)*tempd
      v2d(1) = v2d(1) - v1(2)*tempd
      tempd = fact*sid(i, j, k, 2)
      sid(i, j, k, 2) = 0.0_8
      v1d(3) = v1d(3) + v2(1)*tempd
      v2d(1) = v2d(1) + v1(3)*tempd
      v1d(1) = v1d(1) - v2(3)*tempd
      v2d(3) = v2d(3) - v1(1)*tempd
      tempd = fact*sid(i, j, k, 1)
      sid(i, j, k, 1) = 0.0_8
      v1d(2) = v1d(2) + v2(3)*tempd
      v2d(3) = v2d(3) + v1(2)*tempd
      v1d(3) = v1d(3) - v2(2)*tempd
      v2d(2) = v2d(2) - v1(3)*tempd
      xd(i, j, k, 3) = xd(i, j, k, 3) + v2d(3)
      xd(i, m, n, 3) = xd(i, m, n, 3) - v2d(3)
      v2d(3) = 0.0_8
      xd(i, j, k, 2) = xd(i, j, k, 2) + v2d(2)
      xd(i, m, n, 2) = xd(i, m, n, 2) - v2d(2)
      v2d(2) = 0.0_8
      xd(i, j, k, 1) = xd(i, j, k, 1) + v2d(1)
      xd(i, m, n, 1) = xd(i, m, n, 1) - v2d(1)
      v2d(1) = 0.0_8
      xd(i, j, n, 3) = xd(i, j, n, 3) + v1d(3)
      xd(i, m, k, 3) = xd(i, m, k, 3) - v1d(3)
      v1d(3) = 0.0_8
      xd(i, j, n, 2) = xd(i, j, n, 2) + v1d(2)
      xd(i, m, k, 2) = xd(i, m, k, 2) - v1d(2)
      v1d(2) = 0.0_8
      xd(i, j, n, 1) = xd(i, j, n, 1) + v1d(1)
      xd(i, m, k, 1) = xd(i, m, k, 1) - v1d(1)
      v1d(1) = 0.0_8
    end do
  end subroutine metric_block_b
 
  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 reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *si *sj *sk *(*bcdata.norm)
!   with respect to varying inputs: *si *sj *sk *(*bcdata.norm)
!   rw status of diff variables: *si:incr *sj:incr *sk:incr *(*bcdata.norm):in-out
!   plus diff mem management of: si:in sj:in sk:in bcdata:in *bcdata.norm:in
  subroutine boundarynormals_b()
!  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) :: factd
    real(kind=realtype) :: xxp, yyp, zzp
    real(kind=realtype) :: xxpd, yypd, zzpd
    intrinsic mod
    intrinsic sqrt
    real(kind=realtype) :: tempd
    integer :: branch
    zzpd = 0.0_8
    yypd = 0.0_8
    xxpd = 0.0_8
!$bwd-of ii-loop 
    do mm=1,nbocos
!$bwd-of ii-loop 
      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)
          call pushcontrol3b(1)
        case (imax) 
          mult = one
          xxp = si(il, i, j, 1)
          yyp = si(il, i, j, 2)
          zzp = si(il, i, j, 3)
          call pushcontrol3b(2)
        case (jmin) 
          mult = -one
          xxp = sj(i, 1, j, 1)
          yyp = sj(i, 1, j, 2)
          zzp = sj(i, 1, j, 3)
          call pushcontrol3b(3)
        case (jmax) 
          mult = one
          xxp = sj(i, jl, j, 1)
          yyp = sj(i, jl, j, 2)
          zzp = sj(i, jl, j, 3)
          call pushcontrol3b(4)
        case (kmin) 
          mult = -one
          xxp = sk(i, j, 1, 1)
          yyp = sk(i, j, 1, 2)
          zzp = sk(i, j, 1, 3)
          call pushcontrol3b(5)
        case (kmax) 
          mult = one
          xxp = sk(i, j, kl, 1)
          yyp = sk(i, j, kl, 2)
          zzp = sk(i, j, kl, 3)
          call pushcontrol3b(6)
        case default
          call pushcontrol3b(0)
        end select
! compute the inverse of the length of the normal vector
! and possibly correct for inward pointing.
        fact = sqrt(xxp*xxp + yyp*yyp + zzp*zzp)
        if (fact .gt. zero) then
          call pushreal8(fact)
          fact = mult/fact
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
        end if
        factd = zzp*bcdatad(mm)%norm(i, j, 3)
        zzpd = zzpd + fact*bcdatad(mm)%norm(i, j, 3)
        bcdatad(mm)%norm(i, j, 3) = 0.0_8
        factd = factd + yyp*bcdatad(mm)%norm(i, j, 2)
        yypd = yypd + fact*bcdatad(mm)%norm(i, j, 2)
        bcdatad(mm)%norm(i, j, 2) = 0.0_8
        factd = factd + xxp*bcdatad(mm)%norm(i, j, 1)
        xxpd = xxpd + fact*bcdatad(mm)%norm(i, j, 1)
        bcdatad(mm)%norm(i, j, 1) = 0.0_8
        call popcontrol1b(branch)
        if (branch .eq. 0) then
          call popreal8(fact)
          factd = -(mult*factd/fact**2)
        end if
        if (xxp**2 + yyp**2 + zzp**2 .eq. 0.0_8) then
          tempd = 0.0_8
        else
          tempd = factd/(2.0*sqrt(xxp**2+yyp**2+zzp**2))
        end if
        xxpd = xxpd + 2*xxp*tempd
        yypd = yypd + 2*yyp*tempd
        zzpd = zzpd + 2*zzp*tempd
        call popcontrol3b(branch)
        if (branch .lt. 3) then
          if (branch .ne. 0) then
            if (branch .eq. 1) then
              sid(1, i, j, 3) = sid(1, i, j, 3) + zzpd
              sid(1, i, j, 2) = sid(1, i, j, 2) + yypd
              sid(1, i, j, 1) = sid(1, i, j, 1) + xxpd
              zzpd = 0.0_8
              yypd = 0.0_8
              xxpd = 0.0_8
            else
              sid(il, i, j, 3) = sid(il, i, j, 3) + zzpd
              sid(il, i, j, 2) = sid(il, i, j, 2) + yypd
              sid(il, i, j, 1) = sid(il, i, j, 1) + xxpd
              zzpd = 0.0_8
              yypd = 0.0_8
              xxpd = 0.0_8
            end if
          end if
        else if (branch .lt. 5) then
          if (branch .eq. 3) then
            sjd(i, 1, j, 3) = sjd(i, 1, j, 3) + zzpd
            sjd(i, 1, j, 2) = sjd(i, 1, j, 2) + yypd
            sjd(i, 1, j, 1) = sjd(i, 1, j, 1) + xxpd
            zzpd = 0.0_8
            yypd = 0.0_8
            xxpd = 0.0_8
          else
            sjd(i, jl, j, 3) = sjd(i, jl, j, 3) + zzpd
            sjd(i, jl, j, 2) = sjd(i, jl, j, 2) + yypd
            sjd(i, jl, j, 1) = sjd(i, jl, j, 1) + xxpd
            zzpd = 0.0_8
            yypd = 0.0_8
            xxpd = 0.0_8
          end if
        else if (branch .eq. 5) then
          skd(i, j, 1, 3) = skd(i, j, 1, 3) + zzpd
          skd(i, j, 1, 2) = skd(i, j, 1, 2) + yypd
          skd(i, j, 1, 1) = skd(i, j, 1, 1) + xxpd
          zzpd = 0.0_8
          yypd = 0.0_8
          xxpd = 0.0_8
        else
          skd(i, j, kl, 3) = skd(i, j, kl, 3) + zzpd
          skd(i, j, kl, 2) = skd(i, j, kl, 2) + yypd
          skd(i, j, kl, 1) = skd(i, j, kl, 1) + xxpd
          zzpd = 0.0_8
          yypd = 0.0_8
          xxpd = 0.0_8
        end if
      end do
    end do
  end subroutine boundarynormals_b
 
  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
!$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.
        fact = sqrt(xxp*xxp + yyp*yyp + zzp*zzp)
        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 reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     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_b()
!
!       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) :: tmp
    real(kind=realtype) :: tmpd
    real(kind=realtype) :: tmp0
    real(kind=realtype) :: tmpd0
    real(kind=realtype) :: tmp1
    real(kind=realtype) :: tmpd1
    real(kind=realtype) :: tmp2
    real(kind=realtype) :: tmpd2
    real(kind=realtype) :: tmp3
    real(kind=realtype) :: tmpd3
    real(kind=realtype) :: tmp4
    real(kind=realtype) :: tmpd4
    real(kind=realtype) :: tmp5
    real(kind=realtype) :: tmpd5
    real(kind=realtype) :: tmp6
    real(kind=realtype) :: tmpd6
    real(kind=realtype) :: tmp7
    real(kind=realtype) :: tmpd7
    real(kind=realtype) :: tempd
    real(kind=realtype) :: tmp8
    real(kind=realtype) :: tmpd8
    real(kind=realtype) :: tmp9
    real(kind=realtype) :: tmpd9
    real(kind=realtype) :: tmp10
    real(kind=realtype) :: tmpd10
    real(kind=realtype) :: tmp11
    real(kind=realtype) :: tmpd11
    real(kind=realtype) :: tmp12
    real(kind=realtype) :: tmpd12
    real(kind=realtype) :: tmp13
    real(kind=realtype) :: tmpd13
    real(kind=realtype) :: tmp14
    real(kind=realtype) :: tmpd14
    real(kind=realtype) :: tmp15
    real(kind=realtype) :: tmpd15
    real(kind=realtype) :: tmp16
    real(kind=realtype) :: tmpd16
    integer :: ad_from
    integer :: ad_to
    integer :: ad_from0
    integer :: ad_to0
    integer :: ad_from1
    integer :: ad_to1
    integer :: ad_from2
    integer :: ad_to2
    integer :: ad_from3
    integer :: ad_to3
    integer :: ad_from4
    integer :: ad_to4
    integer :: ad_from5
    integer :: ad_to5
    integer :: ad_from6
    integer :: ad_to6
    integer :: ad_from7
    integer :: ad_to7
    integer :: ad_from8
    integer :: ad_to8
    integer :: ad_from9
    integer :: ad_to9
    integer :: ad_from10
    integer :: ad_to10
    integer :: branch
!
!           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.
        call pushreal8(norm(1))
        norm(1) = bcdata(mm)%symnorm(1)
        call pushreal8(norm(2))
        norm(2) = bcdata(mm)%symnorm(2)
        call pushreal8(norm(3))
        norm(3) = bcdata(mm)%symnorm(3)
        length = sqrt(norm(1)**2 + norm(2)**2 + norm(3)**2)
! compute the unit normal of the subface.
        call pushreal8(norm(1))
        norm(1) = norm(1)/length
        call pushreal8(norm(2))
        norm(2) = norm(2)/length
        call pushreal8(norm(3))
        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
            ad_from0 = jbeg
            do j=ad_from0,jend
              ad_from = ibeg
              i = iend + 1
              call pushinteger4(i - 1)
              call pushinteger4(ad_from)
            end do
            call pushinteger4(j - 1)
            call pushinteger4(ad_from0)
            call pushcontrol4b(7)
          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
            ad_from2 = jbeg
            do j=ad_from2,jend
              ad_from1 = ibeg
              i = iend + 1
              call pushinteger4(i - 1)
              call pushinteger4(ad_from1)
            end do
            call pushinteger4(j - 1)
            call pushinteger4(ad_from2)
            call pushcontrol4b(6)
          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
            ad_from4 = jbeg
            do j=ad_from4,jend
              ad_from3 = ibeg
              i = iend + 1
              call pushinteger4(i - 1)
              call pushinteger4(ad_from3)
            end do
            call pushinteger4(j - 1)
            call pushinteger4(ad_from4)
            call pushcontrol4b(5)
          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
            ad_from6 = jbeg
            do j=ad_from6,jend
              ad_from5 = ibeg
              i = iend + 1
              call pushinteger4(i - 1)
              call pushinteger4(ad_from5)
            end do
            call pushinteger4(j - 1)
            call pushinteger4(ad_from6)
            call pushcontrol4b(4)
          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
            ad_from8 = jbeg
            do j=ad_from8,jend
              ad_from7 = ibeg
              i = iend + 1
              call pushinteger4(i - 1)
              call pushinteger4(ad_from7)
            end do
            call pushinteger4(j - 1)
            call pushinteger4(ad_from8)
            call pushcontrol4b(3)
          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
            ad_from10 = jbeg
            do j=ad_from10,jend
              ad_from9 = ibeg
              i = iend + 1
              call pushinteger4(i - 1)
              call pushinteger4(ad_from9)
            end do
            call pushinteger4(j - 1)
            call pushinteger4(ad_from10)
            call pushcontrol4b(2)
          case default
            call pushcontrol4b(8)
          end select
        else
          call pushcontrol4b(1)
        end if
      else
        call pushcontrol4b(0)
      end if
    end do loopbocos
    v1d = 0.0_8
    do 100 mm=nbocos,1,-1
      call popcontrol4b(branch)
      if (branch .lt. 4) then
        if (branch .lt. 2) then
          if (branch .eq. 0) goto 100
        else if (branch .eq. 2) then
          call popinteger4(ad_from10)
          call popinteger4(ad_to10)
          do j=ad_to10,ad_from10,-1
            call popinteger4(ad_from9)
            call popinteger4(ad_to9)
            do i=ad_to9,ad_from9,-1
              tmpd16 = xd(i, j, ke, 3)
              xd(i, j, ke, 3) = 0.0_8
              xd(i, j, nz, 3) = xd(i, j, nz, 3) + tmpd16
              tmpd15 = xd(i, j, ke, 2)
              xd(i, j, ke, 2) = 0.0_8
              xd(i, j, nz, 2) = xd(i, j, nz, 2) + tmpd15
              tmpd14 = xd(i, j, ke, 1)
              dotd = norm(3)*tmpd16 + norm(2)*tmpd15 + norm(1)*tmpd14
              xd(i, j, ke, 1) = 0.0_8
              xd(i, j, nz, 1) = xd(i, j, nz, 1) + tmpd14
              tempd = two*dotd
              v1d(1) = v1d(1) + norm(1)*tempd
              v1d(2) = v1d(2) + norm(2)*tempd
              v1d(3) = v1d(3) + norm(3)*tempd
              xd(i, j, kl, 3) = xd(i, j, kl, 3) + v1d(3)
              xd(i, j, nz, 3) = xd(i, j, nz, 3) - v1d(3)
              v1d(3) = 0.0_8
              xd(i, j, kl, 2) = xd(i, j, kl, 2) + v1d(2)
              xd(i, j, nz, 2) = xd(i, j, nz, 2) - v1d(2)
              v1d(2) = 0.0_8
              xd(i, j, kl, 1) = xd(i, j, kl, 1) + v1d(1)
              xd(i, j, nz, 1) = xd(i, j, nz, 1) - v1d(1)
              v1d(1) = 0.0_8
            end do
          end do
        else
          call popinteger4(ad_from8)
          call popinteger4(ad_to8)
          do j=ad_to8,ad_from8,-1
            call popinteger4(ad_from7)
            call popinteger4(ad_to7)
            do i=ad_to7,ad_from7,-1
              xd(i, j, 2, 3) = xd(i, j, 2, 3) + xd(i, j, 0, 3)
              dotd = norm(3)*xd(i, j, 0, 3)
              xd(i, j, 0, 3) = 0.0_8
              xd(i, j, 2, 2) = xd(i, j, 2, 2) + xd(i, j, 0, 2)
              dotd = dotd + norm(2)*xd(i, j, 0, 2)
              xd(i, j, 0, 2) = 0.0_8
              xd(i, j, 2, 1) = xd(i, j, 2, 1) + xd(i, j, 0, 1)
              dotd = dotd + norm(1)*xd(i, j, 0, 1)
              xd(i, j, 0, 1) = 0.0_8
              tempd = two*dotd
              v1d(1) = v1d(1) + norm(1)*tempd
              v1d(2) = v1d(2) + norm(2)*tempd
              v1d(3) = v1d(3) + norm(3)*tempd
              xd(i, j, 1, 3) = xd(i, j, 1, 3) + v1d(3)
              xd(i, j, 2, 3) = xd(i, j, 2, 3) - v1d(3)
              v1d(3) = 0.0_8
              xd(i, j, 1, 2) = xd(i, j, 1, 2) + v1d(2)
              xd(i, j, 2, 2) = xd(i, j, 2, 2) - v1d(2)
              v1d(2) = 0.0_8
              xd(i, j, 1, 1) = xd(i, j, 1, 1) + v1d(1)
              xd(i, j, 2, 1) = xd(i, j, 2, 1) - v1d(1)
              v1d(1) = 0.0_8
            end do
          end do
        end if
      else if (branch .lt. 6) then
        if (branch .eq. 4) then
          call popinteger4(ad_from6)
          call popinteger4(ad_to6)
          do j=ad_to6,ad_from6,-1
            call popinteger4(ad_from5)
            call popinteger4(ad_to5)
            do i=ad_to5,ad_from5,-1
              tmpd13 = xd(i, je, j, 3)
              xd(i, je, j, 3) = 0.0_8
              xd(i, ny, j, 3) = xd(i, ny, j, 3) + tmpd13
              tmpd12 = xd(i, je, j, 2)
              xd(i, je, j, 2) = 0.0_8
              xd(i, ny, j, 2) = xd(i, ny, j, 2) + tmpd12
              tmpd11 = xd(i, je, j, 1)
              dotd = norm(3)*tmpd13 + norm(2)*tmpd12 + norm(1)*tmpd11
              xd(i, je, j, 1) = 0.0_8
              xd(i, ny, j, 1) = xd(i, ny, j, 1) + tmpd11
              tempd = two*dotd
              v1d(1) = v1d(1) + norm(1)*tempd
              v1d(2) = v1d(2) + norm(2)*tempd
              v1d(3) = v1d(3) + norm(3)*tempd
              xd(i, jl, j, 3) = xd(i, jl, j, 3) + v1d(3)
              xd(i, ny, j, 3) = xd(i, ny, j, 3) - v1d(3)
              v1d(3) = 0.0_8
              xd(i, jl, j, 2) = xd(i, jl, j, 2) + v1d(2)
              xd(i, ny, j, 2) = xd(i, ny, j, 2) - v1d(2)
              v1d(2) = 0.0_8
              xd(i, jl, j, 1) = xd(i, jl, j, 1) + v1d(1)
              xd(i, ny, j, 1) = xd(i, ny, j, 1) - v1d(1)
              v1d(1) = 0.0_8
            end do
          end do
        else
          call popinteger4(ad_from4)
          call popinteger4(ad_to4)
          do j=ad_to4,ad_from4,-1
            call popinteger4(ad_from3)
            call popinteger4(ad_to3)
            do i=ad_to3,ad_from3,-1
              xd(i, 2, j, 3) = xd(i, 2, j, 3) + xd(i, 0, j, 3)
              dotd = norm(3)*xd(i, 0, j, 3)
              xd(i, 0, j, 3) = 0.0_8
              xd(i, 2, j, 2) = xd(i, 2, j, 2) + xd(i, 0, j, 2)
              dotd = dotd + norm(2)*xd(i, 0, j, 2)
              xd(i, 0, j, 2) = 0.0_8
              xd(i, 2, j, 1) = xd(i, 2, j, 1) + xd(i, 0, j, 1)
              dotd = dotd + norm(1)*xd(i, 0, j, 1)
              xd(i, 0, j, 1) = 0.0_8
              tempd = two*dotd
              v1d(1) = v1d(1) + norm(1)*tempd
              v1d(2) = v1d(2) + norm(2)*tempd
              v1d(3) = v1d(3) + norm(3)*tempd
              xd(i, 1, j, 3) = xd(i, 1, j, 3) + v1d(3)
              xd(i, 2, j, 3) = xd(i, 2, j, 3) - v1d(3)
              v1d(3) = 0.0_8
              xd(i, 1, j, 2) = xd(i, 1, j, 2) + v1d(2)
              xd(i, 2, j, 2) = xd(i, 2, j, 2) - v1d(2)
              v1d(2) = 0.0_8
              xd(i, 1, j, 1) = xd(i, 1, j, 1) + v1d(1)
              xd(i, 2, j, 1) = xd(i, 2, j, 1) - v1d(1)
              v1d(1) = 0.0_8
            end do
          end do
        end if
      else if (branch .eq. 6) then
        call popinteger4(ad_from2)
        call popinteger4(ad_to2)
        do j=ad_to2,ad_from2,-1
          call popinteger4(ad_from1)
          call popinteger4(ad_to1)
          do i=ad_to1,ad_from1,-1
            tmpd10 = xd(ie, i, j, 3)
            xd(ie, i, j, 3) = 0.0_8
            xd(nx, i, j, 3) = xd(nx, i, j, 3) + tmpd10
            tmpd9 = xd(ie, i, j, 2)
            xd(ie, i, j, 2) = 0.0_8
            xd(nx, i, j, 2) = xd(nx, i, j, 2) + tmpd9
            tmpd8 = xd(ie, i, j, 1)
            dotd = norm(3)*tmpd10 + norm(2)*tmpd9 + norm(1)*tmpd8
            xd(ie, i, j, 1) = 0.0_8
            xd(nx, i, j, 1) = xd(nx, i, j, 1) + tmpd8
            tempd = two*dotd
            v1d(1) = v1d(1) + norm(1)*tempd
            v1d(2) = v1d(2) + norm(2)*tempd
            v1d(3) = v1d(3) + norm(3)*tempd
            xd(il, i, j, 3) = xd(il, i, j, 3) + v1d(3)
            xd(nx, i, j, 3) = xd(nx, i, j, 3) - v1d(3)
            v1d(3) = 0.0_8
            xd(il, i, j, 2) = xd(il, i, j, 2) + v1d(2)
            xd(nx, i, j, 2) = xd(nx, i, j, 2) - v1d(2)
            v1d(2) = 0.0_8
            xd(il, i, j, 1) = xd(il, i, j, 1) + v1d(1)
            xd(nx, i, j, 1) = xd(nx, i, j, 1) - v1d(1)
            v1d(1) = 0.0_8
          end do
        end do
      else if (branch .eq. 7) then
        call popinteger4(ad_from0)
        call popinteger4(ad_to0)
        do j=ad_to0,ad_from0,-1
          call popinteger4(ad_from)
          call popinteger4(ad_to)
          do i=ad_to,ad_from,-1
            xd(2, i, j, 3) = xd(2, i, j, 3) + xd(0, i, j, 3)
            dotd = norm(3)*xd(0, i, j, 3)
            xd(0, i, j, 3) = 0.0_8
            xd(2, i, j, 2) = xd(2, i, j, 2) + xd(0, i, j, 2)
            dotd = dotd + norm(2)*xd(0, i, j, 2)
            xd(0, i, j, 2) = 0.0_8
            xd(2, i, j, 1) = xd(2, i, j, 1) + xd(0, i, j, 1)
            dotd = dotd + norm(1)*xd(0, i, j, 1)
            xd(0, i, j, 1) = 0.0_8
            tempd = two*dotd
            v1d(1) = v1d(1) + norm(1)*tempd
            v1d(2) = v1d(2) + norm(2)*tempd
            v1d(3) = v1d(3) + norm(3)*tempd
            xd(1, i, j, 3) = xd(1, i, j, 3) + v1d(3)
            xd(2, i, j, 3) = xd(2, i, j, 3) - v1d(3)
            v1d(3) = 0.0_8
            xd(1, i, j, 2) = xd(1, i, j, 2) + v1d(2)
            xd(2, i, j, 2) = xd(2, i, j, 2) - v1d(2)
            v1d(2) = 0.0_8
            xd(1, i, j, 1) = xd(1, i, j, 1) + v1d(1)
            xd(2, i, j, 1) = xd(2, i, j, 1) - v1d(1)
            v1d(1) = 0.0_8
          end do
        end do
      end if
      call popreal8(norm(3))
      call popreal8(norm(2))
      call popreal8(norm(1))
      call popreal8(norm(3))
      call popreal8(norm(2))
      call popreal8(norm(1))
 100 continue
    do j=je,0,-1
      do i=ie,0,-1
        tmpd7 = xd(i, j, ke, 3)
        xd(i, j, ke, 3) = 0.0_8
        xd(i, j, kl, 3) = xd(i, j, kl, 3) + two*tmpd7
        xd(i, j, nz, 3) = xd(i, j, nz, 3) - tmpd7
        tmpd6 = xd(i, j, ke, 2)
        xd(i, j, ke, 2) = 0.0_8
        xd(i, j, kl, 2) = xd(i, j, kl, 2) + two*tmpd6
        xd(i, j, nz, 2) = xd(i, j, nz, 2) - tmpd6
        tmpd5 = xd(i, j, ke, 1)
        xd(i, j, ke, 1) = 0.0_8
        xd(i, j, kl, 1) = xd(i, j, kl, 1) + two*tmpd5
        xd(i, j, nz, 1) = xd(i, j, nz, 1) - tmpd5
        xd(i, j, 1, 3) = xd(i, j, 1, 3) + two*xd(i, j, 0, 3)
        xd(i, j, 2, 3) = xd(i, j, 2, 3) - xd(i, j, 0, 3)
        xd(i, j, 0, 3) = 0.0_8
        xd(i, j, 1, 2) = xd(i, j, 1, 2) + two*xd(i, j, 0, 2)
        xd(i, j, 2, 2) = xd(i, j, 2, 2) - xd(i, j, 0, 2)
        xd(i, j, 0, 2) = 0.0_8
        xd(i, j, 1, 1) = xd(i, j, 1, 1) + two*xd(i, j, 0, 1)
        xd(i, j, 2, 1) = xd(i, j, 2, 1) - xd(i, j, 0, 1)
        xd(i, j, 0, 1) = 0.0_8
      end do
    end do
    do k=kl,1,-1
      do i=ie,0,-1
        tmpd4 = xd(i, je, k, 3)
        xd(i, je, k, 3) = 0.0_8
        xd(i, jl, k, 3) = xd(i, jl, k, 3) + two*tmpd4
        xd(i, ny, k, 3) = xd(i, ny, k, 3) - tmpd4
        tmpd3 = xd(i, je, k, 2)
        xd(i, je, k, 2) = 0.0_8
        xd(i, jl, k, 2) = xd(i, jl, k, 2) + two*tmpd3
        xd(i, ny, k, 2) = xd(i, ny, k, 2) - tmpd3
        tmpd2 = xd(i, je, k, 1)
        xd(i, je, k, 1) = 0.0_8
        xd(i, jl, k, 1) = xd(i, jl, k, 1) + two*tmpd2
        xd(i, ny, k, 1) = xd(i, ny, k, 1) - tmpd2
        xd(i, 1, k, 3) = xd(i, 1, k, 3) + two*xd(i, 0, k, 3)
        xd(i, 2, k, 3) = xd(i, 2, k, 3) - xd(i, 0, k, 3)
        xd(i, 0, k, 3) = 0.0_8
        xd(i, 1, k, 2) = xd(i, 1, k, 2) + two*xd(i, 0, k, 2)
        xd(i, 2, k, 2) = xd(i, 2, k, 2) - xd(i, 0, k, 2)
        xd(i, 0, k, 2) = 0.0_8
        xd(i, 1, k, 1) = xd(i, 1, k, 1) + two*xd(i, 0, k, 1)
        xd(i, 2, k, 1) = xd(i, 2, k, 1) - xd(i, 0, k, 1)
        xd(i, 0, k, 1) = 0.0_8
      end do
    end do
    do k=kl,1,-1
      do j=jl,1,-1
        tmpd1 = xd(ie, j, k, 3)
        xd(ie, j, k, 3) = 0.0_8
        xd(il, j, k, 3) = xd(il, j, k, 3) + two*tmpd1
        xd(nx, j, k, 3) = xd(nx, j, k, 3) - tmpd1
        tmpd0 = xd(ie, j, k, 2)
        xd(ie, j, k, 2) = 0.0_8
        xd(il, j, k, 2) = xd(il, j, k, 2) + two*tmpd0
        xd(nx, j, k, 2) = xd(nx, j, k, 2) - tmpd0
        tmpd = xd(ie, j, k, 1)
        xd(ie, j, k, 1) = 0.0_8
        xd(il, j, k, 1) = xd(il, j, k, 1) + two*tmpd
        xd(nx, j, k, 1) = xd(nx, j, k, 1) - tmpd
        xd(1, j, k, 3) = xd(1, j, k, 3) + two*xd(0, j, k, 3)
        xd(2, j, k, 3) = xd(2, j, k, 3) - xd(0, j, k, 3)
        xd(0, j, k, 3) = 0.0_8
        xd(1, j, k, 2) = xd(1, j, k, 2) + two*xd(0, j, k, 2)
        xd(2, j, k, 2) = xd(2, j, k, 2) - xd(0, j, k, 2)
        xd(0, j, k, 2) = 0.0_8
        xd(1, j, k, 1) = xd(1, j, k, 1) + two*xd(0, j, k, 1)
        xd(2, j, k, 1) = xd(2, j, k, 1) - xd(0, j, k, 1)
        xd(0, j, k, 1) = 0.0_8
      end do
    end do
  end subroutine xhalo_block_b
 
  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
! 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)
        length = sqrt(norm(1)**2 + norm(2)**2 + norm(3)**2)
! 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 reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     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_b()
    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
    intrinsic mod
! divide through by the reference volume
    nturb = nt2 - nt1 + 1
!$bwd-of ii-loop 
    do ii=0,nx*ny*nz-1
      i = mod(ii, nx) + 2
      j = mod(ii/nx, ny) + 2
      k = ii/(nx*ny) + 2
      ovol = one/volref(i, j, k)
      dwd(i, j, k, nt1:nt2) = ovol*turbresscale(1:nturb)*dwd(i, j, k, &
&       nt1:nt2)
      dwd(i, j, k, 1:nwf) = ovol*dwd(i, j, k, 1:nwf)
    end do
  end subroutine resscale_b
 
  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
    intrinsic mod
! divide through by the reference volume
    nturb = nt2 - nt1 + 1
!$ad ii-loop
    do ii=0,nx*ny*nz-1
      i = mod(ii, nx) + 2
      j = mod(ii/nx, ny) + 2
      k = ii/(nx*ny) + 2
      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 subroutine resscale
 
!  differentiation of sumdwandfw in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *dw
!   with respect to varying inputs: *dw *fw
!   rw status of diff variables: *dw:in-out *fw:out
!   plus diff mem management of: dw:in fw:in
  subroutine sumdwandfw_b()
    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
    integer :: branch
    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
              call pushreal8(max1)
              max1 = zero
              call pushcontrol1b(0)
            else
              call pushreal8(max1)
              max1 = x1
              call pushcontrol1b(1)
            end if
          end do
        end do
      end do
    end do
    if (associated(fwd)) fwd = 0.0_8
    do l=nwf,1,-1
      do k=kl,2,-1
        do j=jl,2,-1
          do i=il,2,-1
            fwd(i, j, k, l) = fwd(i, j, k, l) + max1*dwd(i, j, k, l)
            dwd(i, j, k, l) = max1*dwd(i, j, k, l)
            call popcontrol1b(branch)
            if (branch .eq. 0) then
              call popreal8(max1)
            else
              call popreal8(max1)
            end if
          end do
        end do
      end do
    end do
  end subroutine sumdwandfw_b
 
  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_b