flowUtils_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 flowutils_b
  implicit none
 
contains
!  differentiation of computettot in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: rgas p u v w ttot rho
!   with respect to varying inputs: rgas p u v w ttot rho
  subroutine computettot_b(rho, rhod, u, ud, v, vd, w, wd, p, pd, ttot, &
&   ttotd)
!
!       computettot computes the total temperature for the given
!       pressures, densities and velocities.
!
    use constants
    use inputphysics, only : cpmodel
    use flowvarrefstate, only : rgas, rgasd, gammainf
    use utils_b, only : terminate
    implicit none
!
!      subroutine arguments.
!
    real(kind=realtype), intent(in) :: rho, p, u, v, w
    real(kind=realtype) :: rhod, pd, ud, vd, wd
    real(kind=realtype) :: ttot
    real(kind=realtype) :: ttotd
!
!      local variables.
!
    integer(kind=inttype) :: i
    real(kind=realtype) :: govgm1, t, kin
    real(kind=realtype) :: td, kind0
    real(kind=realtype) :: tempd
    real(kind=realtype) :: temp
    real(kind=realtype) :: tempd0
! determine the cp model used.
    select case  (cpmodel) 
    case (cpconstant) 
! constant cp and thus constant gamma. the well-known
! formula is valid.
      govgm1 = gammainf/(gammainf-one)
      t = p/(rho*rgas)
      kin = half*(u*u+v*v+w*w)
!===============================================================
      temp = rho*kin/(govgm1*p)
      td = (one+temp)*ttotd
      tempd0 = t*ttotd/(govgm1*p)
      rhod = rhod + kin*tempd0
      kind0 = rho*tempd0
      tempd = td/(rho*rgas)
      pd = pd + tempd - govgm1*temp*tempd0
      tempd0 = half*kind0
      ud = ud + 2*u*tempd0
      vd = vd + 2*v*tempd0
      wd = wd + 2*w*tempd0
      tempd0 = -(p*tempd/(rho*rgas))
      rhod = rhod + rgas*tempd0
      rgasd = rgasd + rho*tempd0
      ttotd = 0.0_8
    end select
  end subroutine computettot_b
 
  subroutine computettot(rho, u, v, w, p, ttot)
!
!       computettot computes the total temperature for the given
!       pressures, densities and velocities.
!
    use constants
    use inputphysics, only : cpmodel
    use flowvarrefstate, only : rgas, gammainf
    use utils_b, only : terminate
    implicit none
!
!      subroutine arguments.
!
    real(kind=realtype), intent(in) :: rho, p, u, v, w
    real(kind=realtype), intent(out) :: ttot
!
!      local variables.
!
    integer(kind=inttype) :: i
    real(kind=realtype) :: govgm1, t, kin
! determine the cp model used.
    select case  (cpmodel) 
    case (cpconstant) 
! constant cp and thus constant gamma. the well-known
! formula is valid.
      govgm1 = gammainf/(gammainf-one)
      t = p/(rho*rgas)
      kin = half*(u*u+v*v+w*w)
      ttot = t*(one+rho*kin/(govgm1*p))
!===============================================================
    case (cptempcurvefits) 
! cp is a function of the temperature. the formula used for
! constant cp is not valid anymore and a more complicated
! procedure must be followed.
      call terminate('computettot', &
&              'variable cp formulation not implemented yet')
    end select
  end subroutine computettot
 
  subroutine computegamma(t, gamma, mm)
!
!       computegamma computes the corresponding values of gamma for
!       the given dimensional temperatures.
!
    use constants
    use cpcurvefits
    use inputphysics, only : cpmodel, gammaconstant
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: mm
    real(kind=realtype), dimension(mm), intent(in) :: t
    real(kind=realtype), dimension(mm), intent(out) :: gamma
!
!      local variables.
!
    integer(kind=inttype) :: i, ii, nn, start
    real(kind=realtype) :: cp, t2
! determine the cp model used in the computation.
    select case  (cpmodel) 
    case (cpconstant) 
! constant cp and thus constant gamma. set the values.
      do i=1,mm
        gamma(i) = gammaconstant
      end do
    end select
!        ================================================================
 
  end subroutine computegamma
 
!  differentiation of computeptot in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: p u v w ptot rho
!   with respect to varying inputs: p u v w ptot rho
  subroutine computeptot_b(rho, rhod, u, ud, v, vd, w, wd, p, pd, ptot, &
&   ptotd)
!
!       computeptot computes the total pressure for the given
!       pressures, densities and velocities.
!
    use constants
    use cpcurvefits
    use flowvarrefstate, only : tref, trefd, rgas, rgasd, gammainf
    use inputphysics, only : cpmodel
    implicit none
    real(kind=realtype), intent(in) :: rho, p, u, v, w
    real(kind=realtype) :: rhod, pd, ud, vd, wd
    real(kind=realtype) :: ptot
    real(kind=realtype) :: ptotd
!
!      local parameters.
!
    real(kind=realtype), parameter :: dtstop=0.01_realtype
!
!      local variables.
!
    integer(kind=inttype) :: i, ii, mm, nn, nnt, start
    real(kind=realtype) :: govgm1, kin
    real(kind=realtype) :: kind0
    real(kind=realtype) :: t, t2, tt, dt, h, htot, cp, scale, alp
    real(kind=realtype) :: intcport, intcportt, intcporttt
    real(kind=realtype) :: temp
    real(kind=realtype) :: tempd
!
! determine the cp model used.
    select case  (cpmodel) 
    case (cpconstant) 
! constant cp and thus constant gamma. the well-known
! formula is valid.
      govgm1 = gammainf/(gammainf-one)
      kin = half*(u*u+v*v+w*w)
!===============================================================
      temp = rho*kin/(govgm1*p)
      if (one + temp .le. 0.0_8 .and. (govgm1 .eq. 0.0_8 .or. govgm1 &
&         .ne. int(govgm1))) then
        tempd = 0.0_8
      else
        tempd = (one+temp)**(govgm1-1)*ptotd
      end if
      pd = pd + (one+temp)**govgm1*ptotd - govgm1*temp*tempd
      rhod = rhod + kin*tempd
      kind0 = rho*tempd
      tempd = half*kind0
      ud = ud + 2*u*tempd
      vd = vd + 2*v*tempd
      wd = wd + 2*w*tempd
      ptotd = 0.0_8
    end select
  end subroutine computeptot_b
 
  subroutine computeptot(rho, u, v, w, p, ptot)
!
!       computeptot computes the total pressure for the given
!       pressures, densities and velocities.
!
    use constants
    use cpcurvefits
    use flowvarrefstate, only : tref, rgas, gammainf
    use inputphysics, only : cpmodel
    implicit none
    real(kind=realtype), intent(in) :: rho, p, u, v, w
    real(kind=realtype), intent(out) :: ptot
!
!      local parameters.
!
    real(kind=realtype), parameter :: dtstop=0.01_realtype
!
!      local variables.
!
    integer(kind=inttype) :: i, ii, mm, nn, nnt, start
    real(kind=realtype) :: govgm1, kin
    real(kind=realtype) :: t, t2, tt, dt, h, htot, cp, scale, alp
    real(kind=realtype) :: intcport, intcportt, intcporttt
!
! determine the cp model used.
    select case  (cpmodel) 
    case (cpconstant) 
! constant cp and thus constant gamma. the well-known
! formula is valid.
      govgm1 = gammainf/(gammainf-one)
      kin = half*(u*u+v*v+w*w)
      ptot = p*(one+rho*kin/(govgm1*p))**govgm1
!===============================================================
    end select
  end subroutine computeptot
 
!  differentiation of computespeedofsoundsquared in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *aa *p *w
!   with respect to varying inputs: *aa *p *w
!   rw status of diff variables: *aa:in-out *p:incr *w:incr
!   plus diff mem management of: aa:in p:in w:in
  subroutine computespeedofsoundsquared_b()
!
!       computespeedofsoundsquared does what it says.
!
    use constants
    use blockpointers, only : ie, je, ke, w, wd, p, pd, aa, aad, gamma
    use utils_b, only : getcorrectfork
    implicit none
!
!      local variables.
!
    real(kind=realtype), parameter :: twothird=two*third
    integer(kind=inttype) :: i, j, k, ii
    real(kind=realtype) :: pp
    real(kind=realtype) :: ppd
    logical :: correctfork
    intrinsic mod
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype) :: tempd
    real(kind=realtype) :: tempd0
! determine if we need to correct for k
    correctfork = getcorrectfork()
    if (correctfork) then
!$bwd-of ii-loop 
      do ii=0,ie*je*ke-1
        i = mod(ii, ie) + 1
        j = mod(ii/ie, je) + 1
        k = ii/(ie*je) + 1
        pp = p(i, j, k) - twothird*w(i, j, k, irho)*w(i, j, k, itu1)
        temp = w(i, j, k, irho)
        tempd = gamma(i, j, k)*aad(i, j, k)/temp
        aad(i, j, k) = 0.0_8
        ppd = tempd
        wd(i, j, k, irho) = wd(i, j, k, irho) - pp*tempd/temp - w(i, j, &
&         k, itu1)*twothird*ppd
        pd(i, j, k) = pd(i, j, k) + ppd
        wd(i, j, k, itu1) = wd(i, j, k, itu1) - w(i, j, k, irho)*&
&         twothird*ppd
      end do
    else
!$bwd-of ii-loop 
      do ii=0,ie*je*ke-1
        i = mod(ii, ie) + 1
        j = mod(ii/ie, je) + 1
        k = ii/(ie*je) + 1
        temp0 = w(i, j, k, irho)
        tempd0 = gamma(i, j, k)*aad(i, j, k)/temp0
        aad(i, j, k) = 0.0_8
        pd(i, j, k) = pd(i, j, k) + tempd0
        wd(i, j, k, irho) = wd(i, j, k, irho) - p(i, j, k)*tempd0/temp0
      end do
    end if
  end subroutine computespeedofsoundsquared_b
 
  subroutine computespeedofsoundsquared()
!
!       computespeedofsoundsquared does what it says.
!
    use constants
    use blockpointers, only : ie, je, ke, w, p, aa, gamma
    use utils_b, only : getcorrectfork
    implicit none
!
!      local variables.
!
    real(kind=realtype), parameter :: twothird=two*third
    integer(kind=inttype) :: i, j, k, ii
    real(kind=realtype) :: pp
    logical :: correctfork
    intrinsic mod
! determine if we need to correct for k
    correctfork = getcorrectfork()
    if (correctfork) then
!$ad ii-loop
      do ii=0,ie*je*ke-1
        i = mod(ii, ie) + 1
        j = mod(ii/ie, je) + 1
        k = ii/(ie*je) + 1
        pp = p(i, j, k) - twothird*w(i, j, k, irho)*w(i, j, k, itu1)
        aa(i, j, k) = gamma(i, j, k)*pp/w(i, j, k, irho)
      end do
    else
!$ad ii-loop
      do ii=0,ie*je*ke-1
        i = mod(ii, ie) + 1
        j = mod(ii/ie, je) + 1
        k = ii/(ie*je) + 1
        aa(i, j, k) = gamma(i, j, k)*p(i, j, k)/w(i, j, k, irho)
      end do
    end if
  end subroutine computespeedofsoundsquared
 
!  differentiation of computeetotblock in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *p *w
!   with respect to varying inputs: *p *w
!   rw status of diff variables: *p:incr *w:in-out
!   plus diff mem management of: p:in w:in
  subroutine computeetotblock_b(istart, iend, jstart, jend, kstart, kend&
&   , correctfork)
!
!       computeetot computes the total energy from the given density,
!       velocity and presssure. for a calorically and thermally
!       perfect gas the well-known expression is used; for only a
!       thermally perfect gas, cp is a function of temperature, curve
!       fits are used and a more complex expression is obtained.
!       it is assumed that the pointers in blockpointers already
!       point to the correct block.
!
    use constants
    use blockpointers, only : w, wd, p, pd
    use flowvarrefstate, only : rgas, rgasd, tref, trefd
    use inputphysics, only : cpmodel, gammaconstant
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: istart, iend, jstart, jend
    integer(kind=inttype), intent(in) :: kstart, kend
    logical, intent(in) :: correctfork
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii, isize, jsize, ksize
    real(kind=realtype) :: ovgm1, factk, scale
    intrinsic mod
    real(kind=realtype) :: tmp
    real(kind=realtype) :: tmpd
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype) :: temp1
    real(kind=realtype) :: temp2
    real(kind=realtype) :: tempd
    real(kind=realtype) :: tmp0
    real(kind=realtype) :: tmpd0
    real(kind=realtype) :: tempd0
!
! determine the cp model used in the computation.
    select case  (cpmodel) 
    case (cpconstant) 
! constant cp and thus constant gamma.
! abbreviate 1/(gamma -1) a bit easier.
      ovgm1 = one/(gammaconstant-one)
! loop over the given range of the block and compute the first
! step of the energy.
      isize = iend - istart + 1
      jsize = jend - jstart + 1
      ksize = kend - kstart + 1
      factk = ovgm1*(five*third-gammaconstant)
      if (correctfork) then
!$bwd-of ii-loop 
        do ii=0,isize*jsize*ksize-1
          i = mod(ii, isize) + istart
          j = mod(ii/isize, jsize) + jstart
          k = ii/(isize*jsize) + kstart
          tmpd = wd(i, j, k, irhoe)
          wd(i, j, k, irhoe) = 0.0_8
          temp = w(i, j, k, ivz)
          temp0 = w(i, j, k, ivy)
          temp1 = w(i, j, k, ivx)
          pd(i, j, k) = pd(i, j, k) + ovgm1*tmpd
          wd(i, j, k, irho) = wd(i, j, k, irho) + ((temp1**2+temp0**2+&
&           temp**2)*half-w(i, j, k, itu1)*factk)*tmpd
          tempd = w(i, j, k, irho)*half*tmpd
          wd(i, j, k, itu1) = wd(i, j, k, itu1) - w(i, j, k, irho)*factk&
&           *tmpd
          wd(i, j, k, ivx) = wd(i, j, k, ivx) + 2*temp1*tempd
          wd(i, j, k, ivy) = wd(i, j, k, ivy) + 2*temp0*tempd
          wd(i, j, k, ivz) = wd(i, j, k, ivz) + 2*temp*tempd
        end do
      else
!$bwd-of ii-loop 
        do ii=0,isize*jsize*ksize-1
          i = mod(ii, isize) + istart
          j = mod(ii/isize, jsize) + jstart
          k = ii/(isize*jsize) + kstart
          tmpd0 = wd(i, j, k, irhoe)
          wd(i, j, k, irhoe) = 0.0_8
          temp2 = w(i, j, k, ivz)
          temp1 = w(i, j, k, ivy)
          temp0 = w(i, j, k, ivx)
          pd(i, j, k) = pd(i, j, k) + ovgm1*tmpd0
          wd(i, j, k, irho) = wd(i, j, k, irho) + (temp0**2+temp1**2+&
&           temp2**2)*half*tmpd0
          tempd0 = w(i, j, k, irho)*half*tmpd0
          wd(i, j, k, ivx) = wd(i, j, k, ivx) + 2*temp0*tempd0
          wd(i, j, k, ivy) = wd(i, j, k, ivy) + 2*temp1*tempd0
          wd(i, j, k, ivz) = wd(i, j, k, ivz) + 2*temp2*tempd0
        end do
      end if
    end select
  end subroutine computeetotblock_b
 
  subroutine computeetotblock(istart, iend, jstart, jend, kstart, kend, &
&   correctfork)
!
!       computeetot computes the total energy from the given density,
!       velocity and presssure. for a calorically and thermally
!       perfect gas the well-known expression is used; for only a
!       thermally perfect gas, cp is a function of temperature, curve
!       fits are used and a more complex expression is obtained.
!       it is assumed that the pointers in blockpointers already
!       point to the correct block.
!
    use constants
    use blockpointers, only : w, p
    use flowvarrefstate, only : rgas, tref
    use inputphysics, only : cpmodel, gammaconstant
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: istart, iend, jstart, jend
    integer(kind=inttype), intent(in) :: kstart, kend
    logical, intent(in) :: correctfork
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii, isize, jsize, ksize
    real(kind=realtype) :: ovgm1, factk, scale
    intrinsic mod
!
! determine the cp model used in the computation.
    select case  (cpmodel) 
    case (cpconstant) 
! constant cp and thus constant gamma.
! abbreviate 1/(gamma -1) a bit easier.
      ovgm1 = one/(gammaconstant-one)
! loop over the given range of the block and compute the first
! step of the energy.
      isize = iend - istart + 1
      jsize = jend - jstart + 1
      ksize = kend - kstart + 1
      factk = ovgm1*(five*third-gammaconstant)
      if (correctfork) then
!$ad ii-loop
        do ii=0,isize*jsize*ksize-1
          i = mod(ii, isize) + istart
          j = mod(ii/isize, jsize) + jstart
          k = ii/(isize*jsize) + kstart
          w(i, j, k, irhoe) = ovgm1*p(i, j, k) + half*w(i, j, k, irho)*(&
&           w(i, j, k, ivx)**2+w(i, j, k, ivy)**2+w(i, j, k, ivz)**2) - &
&           factk*w(i, j, k, irho)*w(i, j, k, itu1)
        end do
      else
!$ad ii-loop
        do ii=0,isize*jsize*ksize-1
          i = mod(ii, isize) + istart
          j = mod(ii/isize, jsize) + jstart
          k = ii/(isize*jsize) + kstart
          w(i, j, k, irhoe) = ovgm1*p(i, j, k) + half*w(i, j, k, irho)*(&
&           w(i, j, k, ivx)**2+w(i, j, k, ivy)**2+w(i, j, k, ivz)**2)
        end do
      end if
    end select
  end subroutine computeetotblock
 
!  differentiation of etot in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: k p u v w etotal rho
!   with respect to varying inputs: k p u v w etotal rho
  subroutine etot_b(rho, rhod, u, ud, v, vd, w, wd, p, pd, k, kd, etotal&
&   , etotald, correctfork)
!
!       etotarray computes the total energy from the given density,
!       velocity and presssure.
!       first the internal energy per unit mass is computed and after
!       that the kinetic energy is added as well the conversion to
!       energy per unit volume.
!
    use constants
    implicit none
!
!      subroutine arguments.
!
    real(kind=realtype), intent(in) :: rho, p, k
    real(kind=realtype) :: rhod, pd, kd
    real(kind=realtype), intent(in) :: u, v, w
    real(kind=realtype) :: ud, vd, wd
    real(kind=realtype) :: etotal
    real(kind=realtype) :: etotald
    logical, intent(in) :: correctfork
!
!      local variables.
!
    integer(kind=inttype) :: i
    real(kind=realtype) :: tempd
! compute the internal energy for unit mass.
    call pushreal8(etotal)
    call eint(rho, p, k, etotal, correctfork)
    rhod = rhod + (etotal+half*(u**2+v**2+w**2))*etotald
    tempd = half*rho*etotald
    etotald = rho*etotald
    ud = ud + 2*u*tempd
    vd = vd + 2*v*tempd
    wd = wd + 2*w*tempd
    call popreal8(etotal)
    call eint_b(rho, rhod, p, pd, k, kd, etotal, etotald, correctfork)
  end subroutine etot_b
 
  subroutine etot(rho, u, v, w, p, k, etotal, correctfork)
!
!       etotarray computes the total energy from the given density,
!       velocity and presssure.
!       first the internal energy per unit mass is computed and after
!       that the kinetic energy is added as well the conversion to
!       energy per unit volume.
!
    use constants
    implicit none
!
!      subroutine arguments.
!
    real(kind=realtype), intent(in) :: rho, p, k
    real(kind=realtype), intent(in) :: u, v, w
    real(kind=realtype), intent(out) :: etotal
    logical, intent(in) :: correctfork
!
!      local variables.
!
    integer(kind=inttype) :: i
! compute the internal energy for unit mass.
    call eint(rho, p, k, etotal, correctfork)
    etotal = rho*(etotal+half*(u*u+v*v+w*w))
  end subroutine etot
 
!  differentiation of eint in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: k p einternal rho
!   with respect to varying inputs: k p einternal rho
!      ==================================================================
  subroutine eint_b(rho, rhod, p, pd, k, kd, einternal, einternald, &
&   correctfork)
!
!       eintarray computes the internal energy per unit mass from the
!       given density and pressure (and possibly turbulent energy)
!       for a calorically and thermally perfect gas the well-known
!       expression is used; for only a thermally perfect gas, cp is a
!       function of temperature, curve fits are used and a more
!       complex expression is obtained.
!
    use constants
    use cpcurvefits
    use flowvarrefstate, only : rgas, rgasd, tref, trefd
    use inputphysics, only : cpmodel, gammaconstant
    implicit none
!
!      subroutine arguments.
!
    real(kind=realtype), intent(in) :: rho, p, k
    real(kind=realtype) :: rhod, pd, kd
    real(kind=realtype) :: einternal
    real(kind=realtype) :: einternald
    logical, intent(in) :: correctfork
!
!      local parameter.
!
    real(kind=realtype), parameter :: twothird=two*third
!
!      local variables.
!
    integer(kind=inttype) :: i, nn, mm, ii, start
    real(kind=realtype) :: ovgm1, factk, pp, t, t2, scale
    real(kind=realtype) :: tempd
! determine the cp model used in the computation.
    select case  (cpmodel) 
    case (cpconstant) 
! abbreviate 1/(gamma -1) a bit easier.
      ovgm1 = one/(gammaconstant-one)
! loop over the number of elements of the array and compute
! the total energy.
! second step. correct the energy in case a turbulent kinetic
! energy is present.
      if (correctfork) then
        factk = ovgm1*(five*third-gammaconstant)
        kd = kd - factk*einternald
      end if
      tempd = ovgm1*einternald/rho
      pd = pd + tempd
      rhod = rhod - p*tempd/rho
      einternald = 0.0_8
    end select
  end subroutine eint_b
 
!      ==================================================================
  subroutine eint(rho, p, k, einternal, correctfork)
!
!       eintarray computes the internal energy per unit mass from the
!       given density and pressure (and possibly turbulent energy)
!       for a calorically and thermally perfect gas the well-known
!       expression is used; for only a thermally perfect gas, cp is a
!       function of temperature, curve fits are used and a more
!       complex expression is obtained.
!
    use constants
    use cpcurvefits
    use flowvarrefstate, only : rgas, tref
    use inputphysics, only : cpmodel, gammaconstant
    implicit none
!
!      subroutine arguments.
!
    real(kind=realtype), intent(in) :: rho, p, k
    real(kind=realtype), intent(out) :: einternal
    logical, intent(in) :: correctfork
!
!      local parameter.
!
    real(kind=realtype), parameter :: twothird=two*third
!
!      local variables.
!
    integer(kind=inttype) :: i, nn, mm, ii, start
    real(kind=realtype) :: ovgm1, factk, pp, t, t2, scale
! determine the cp model used in the computation.
    select case  (cpmodel) 
    case (cpconstant) 
! abbreviate 1/(gamma -1) a bit easier.
      ovgm1 = one/(gammaconstant-one)
! loop over the number of elements of the array and compute
! the total energy.
      einternal = ovgm1*p/rho
! second step. correct the energy in case a turbulent kinetic
! energy is present.
      if (correctfork) then
        factk = ovgm1*(five*third-gammaconstant)
        einternal = einternal - factk*k
      end if
    end select
  end subroutine eint
 
!  differentiation of computepressuresimple in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: pinfcorr *p *w
!   with respect to varying inputs: pinfcorr *p *w
!   rw status of diff variables: pinfcorr:incr *p:in-out *w:incr
!   plus diff mem management of: p:in w:in
  subroutine computepressuresimple_b(includehalos)
! compute the pressure on a block with the pointers already set. this
! routine is used by the forward mode ad code only.
    use constants
    use blockpointers
    use flowvarrefstate
    use inputphysics
    implicit none
! input parameter
    logical, intent(in) :: includehalos
! local variables
    integer(kind=inttype) :: i, j, k, ii
    real(kind=realtype) :: gm1, v2
    real(kind=realtype) :: v2d
    integer(kind=inttype) :: ibeg, iend, isize, jbeg, jend, jsize, kbeg&
&   , kend, ksize
    intrinsic mod
    intrinsic max
    real(kind=realtype) :: tempd
! compute the pressures
    gm1 = gammaconstant - one
    if (includehalos) then
      ibeg = 0
      jbeg = 0
      kbeg = 0
      iend = ib
      jend = jb
      kend = kb
    else
      ibeg = 2
      jbeg = 2
      kbeg = 2
      iend = il
      jend = jl
      kend = kl
    end if
    isize = iend - ibeg + 1
    jsize = jend - jbeg + 1
    ksize = kend - kbeg + 1
!$bwd-of ii-loop 
    do ii=0,isize*jsize*ksize-1
      i = mod(ii, isize) + ibeg
      j = mod(ii/isize, jsize) + jbeg
      k = ii/(isize*jsize) + kbeg
      v2 = w(i, j, k, ivx)**2 + w(i, j, k, ivy)**2 + w(i, j, k, ivz)**2
      p(i, j, k) = gm1*(w(i, j, k, irhoe)-half*w(i, j, k, irho)*v2)
      if (p(i, j, k) .lt. 1.e-4_realtype*pinfcorr) then
        pinfcorrd = pinfcorrd + 1.e-4_realtype*pd(i, j, k)
        pd(i, j, k) = 0.0_8
      end if
      tempd = gm1*pd(i, j, k)
      pd(i, j, k) = 0.0_8
      wd(i, j, k, irhoe) = wd(i, j, k, irhoe) + tempd
      wd(i, j, k, irho) = wd(i, j, k, irho) - v2*half*tempd
      v2d = -(w(i, j, k, irho)*half*tempd)
      wd(i, j, k, ivx) = wd(i, j, k, ivx) + 2*w(i, j, k, ivx)*v2d
      wd(i, j, k, ivy) = wd(i, j, k, ivy) + 2*w(i, j, k, ivy)*v2d
      wd(i, j, k, ivz) = wd(i, j, k, ivz) + 2*w(i, j, k, ivz)*v2d
    end do
  end subroutine computepressuresimple_b
 
  subroutine computepressuresimple(includehalos)
! compute the pressure on a block with the pointers already set. this
! routine is used by the forward mode ad code only.
    use constants
    use blockpointers
    use flowvarrefstate
    use inputphysics
    implicit none
! input parameter
    logical, intent(in) :: includehalos
! local variables
    integer(kind=inttype) :: i, j, k, ii
    real(kind=realtype) :: gm1, v2
    integer(kind=inttype) :: ibeg, iend, isize, jbeg, jend, jsize, kbeg&
&   , kend, ksize
    intrinsic mod
    intrinsic max
! compute the pressures
    gm1 = gammaconstant - one
    if (includehalos) then
      ibeg = 0
      jbeg = 0
      kbeg = 0
      iend = ib
      jend = jb
      kend = kb
    else
      ibeg = 2
      jbeg = 2
      kbeg = 2
      iend = il
      jend = jl
      kend = kl
    end if
    isize = iend - ibeg + 1
    jsize = jend - jbeg + 1
    ksize = kend - kbeg + 1
!$ad ii-loop
    do ii=0,isize*jsize*ksize-1
      i = mod(ii, isize) + ibeg
      j = mod(ii/isize, jsize) + jbeg
      k = ii/(isize*jsize) + kbeg
      v2 = w(i, j, k, ivx)**2 + w(i, j, k, ivy)**2 + w(i, j, k, ivz)**2
      p(i, j, k) = gm1*(w(i, j, k, irhoe)-half*w(i, j, k, irho)*v2)
      if (p(i, j, k) .lt. 1.e-4_realtype*pinfcorr) then
        p(i, j, k) = 1.e-4_realtype*pinfcorr
      else
        p(i, j, k) = p(i, j, k)
      end if
    end do
  end subroutine computepressuresimple
 
  subroutine computepressure(ibeg, iend, jbeg, jend, kbeg, kend, &
&   pointeroffset)
!
!       computepressure computes the pressure from the total energy,
!       density and velocities in the given cell range of the block to
!       which the pointers in blockpointers currently point.
!       it is possible to specify a possible pointer offset, because
!       this routine is also used when reading a restart file.
!
    use constants
    use inputphysics, only : cpmodel, gammaconstant
    use blockpointers, only : w, p
    use flowvarrefstate, only : kpresent, pinf, rgas, tref
    use cpcurvefits
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: ibeg, iend, jbeg, jend, kbeg, &
&   kend
    integer(kind=inttype), intent(in) :: pointeroffset
!
!      local parameters.
!
    real(kind=realtype), parameter :: dtstop=0.01_realtype
    real(kind=realtype), parameter :: twothird=two*third
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ip, jp, kp, nn, ii, start
    real(kind=realtype) :: gm1, factk, v2, scale, e0, e
    real(kind=realtype) :: trefinv, t, dt, t2, alp, cv
    intrinsic max
    intrinsic log
    intrinsic abs
    real(kind=realtype) :: abs0
! determine the cp model used in the computation.
    select case  (cpmodel) 
    case (cpconstant) 
! constant cp and thus constant gamma. the relation
! eint = cv*t can be used and consequently the standard
! relation between pressure and internal energy is valid.
! abbreviate some constants that occur in the pressure
! computation.
      gm1 = gammaconstant - one
      factk = five*third - gammaconstant
! loop over the cells. take the possible pointer
! offset into account and store the pressure in the
! position w(:,:,:, irhoe).
      do k=kbeg,kend
        kp = k + pointeroffset
        do j=jbeg,jend
          jp = j + pointeroffset
          do i=ibeg,iend
            ip = i + pointeroffset
            v2 = w(ip, jp, kp, ivx)**2 + w(ip, jp, kp, ivy)**2 + w(ip, &
&             jp, kp, ivz)**2
            w(i, j, k, irhoe) = gm1*(w(ip, jp, kp, irhoe)-half*w(ip, jp&
&             , kp, irho)*v2)
            if (w(i, j, k, irhoe) .lt. 1.e-5_realtype*pinf) then
              w(i, j, k, irhoe) = 1.e-5_realtype*pinf
            else
              w(i, j, k, irhoe) = w(i, j, k, irhoe)
            end if
          end do
        end do
      end do
! correct p if a k-equation is present.
      if (kpresent) then
        do k=kbeg,kend
          kp = k + pointeroffset
          do j=jbeg,jend
            jp = j + pointeroffset
            do i=ibeg,iend
              ip = i + pointeroffset
              w(i, j, k, irhoe) = w(i, j, k, irhoe) + factk*w(ip, jp, kp&
&               , irho)*w(ip, jp, kp, itu1)
            end do
          end do
        end do
      end if
    case (cptempcurvefits) 
!        ================================================================
! cp as function of the temperature is given via curve fits.
! store a scale factor when converting the nondimensional
! energy to the units of cpeint
      trefinv = one/tref
      scale = tref/rgas
! loop over the cells to compute the internal energy per
! unit mass. this is stored in w(:,:,:,irhoe) for the moment.
      do k=kbeg,kend
        kp = k + pointeroffset
        do j=jbeg,jend
          jp = j + pointeroffset
          do i=ibeg,iend
            ip = i + pointeroffset
            w(i, j, k, irhoe) = w(ip, jp, kp, irhoe)/w(ip, jp, kp, irho)
            if (kpresent) w(i, j, k, irhoe) = w(i, j, k, irhoe) - w(ip, &
&               jp, kp, itu1)
            v2 = w(ip, jp, kp, ivx)**2 + w(ip, jp, kp, ivy)**2 + w(ip, &
&             jp, kp, ivz)**2
            w(i, j, k, irhoe) = w(i, j, k, irhoe) - half*v2
          end do
        end do
      end do
! newton algorithm to compute the temperature from the known
! value of the internal energy.
      do k=kbeg,kend
        kp = k + pointeroffset
        do j=jbeg,jend
          jp = j + pointeroffset
          do i=ibeg,iend
            ip = i + pointeroffset
! store the internal energy in the same dimensional
! units as cpeint.
            e0 = scale*w(i, j, k, irhoe)
! take care of the exceptional cases.
            if (e0 .le. cpeint(0)) then
! energy smaller than the lowest value of the curve
! fit. use extrapolation using constant cv.
              t = trefinv*(cptrange(0)+(e0-cpeint(0))/cv0)
            else if (e0 .ge. cpeint(cpnparts)) then
! energy larger than the largest value of the curve
! fit. use extrapolation using constant cv.
              t = trefinv*(cptrange(cpnparts)+(e0-cpeint(cpnparts))/cvn)
            else
! the value is in the range of the curve fits.
! a newton algorithm is used to find the temperature.
! first find the curve fit interval to be searched.
              ii = cpnparts
              start = 1
     interval:do 
! next guess for the interval.
                nn = start + ii/2
! determine the situation we are having here.
                if (e0 .gt. cpeint(nn)) then
! energy is larger than the upper boundary of
! the current interval. update the lower
! boundary.
                  start = nn + 1
                  ii = ii - 1
                else if (e0 .ge. cpeint(nn-1)) then
! nn contains the range in which the newton algorithm
! must be applied.
! initial guess of the dimensional temperature.
                  alp = (cpeint(nn)-e0)/(cpeint(nn)-cpeint(nn-1))
                  t = alp*cptrange(nn-1) + (one-alp)*cptrange(nn)
! the actual newton algorithm to compute the
! temperature.
           newton:do 
! compute the internal energy as well as the
! value of cv/r for the given temperature.
! cv/r = cp/r - 1.0
                    cv = -one
! e = integral of cv,
                    e = cptempfit(nn)%eint0 - t
! not of cp.
                    do ii=1,cptempfit(nn)%nterm
! update cv.
                      t2 = t**cptempfit(nn)%exponents(ii)
                      cv = cv + cptempfit(nn)%constants(ii)*t2
! update e, for which this contribution must be
! integrated. take the exceptional case that the
! exponent == -1 into account.
                      if (cptempfit(nn)%exponents(ii) .eq. -1_inttype) &
&                     then
                        e = e + cptempfit(nn)%constants(ii)*log(t)
                      else
                        e = e + cptempfit(nn)%constants(ii)*t2*t/(&
&                         cptempfit(nn)%exponents(ii)+1)
                      end if
                    end do
! compute the update and the new temperature.
                    dt = (e0-e)/cv
                    t = t + dt
                    if (dt .ge. 0.) then
                      abs0 = dt
                    else
                      abs0 = -dt
                    end if
! exit the newton loop if the update is smaller
! than the threshold value.
                    if (abs0 .lt. dtstop) then
! create the nondimensional temperature.
                      t = t*trefinv
                      goto 100
                    end if
                  end do newton
                end if
! this is the correct range. exit the do-loop.
! modify ii for the next branch to search.
                ii = ii/2
              end do interval
            end if
! compute the pressure from the known temperature
! and density. include the correction if a k-equation
! is present.
 100        w(i, j, k, irhoe) = w(ip, jp, kp, irho)*rgas*t
            if (w(i, j, k, irhoe) .lt. 1.e-5_realtype*pinf) then
              w(i, j, k, irhoe) = 1.e-5_realtype*pinf
            else
              w(i, j, k, irhoe) = w(i, j, k, irhoe)
            end if
            if (kpresent) w(i, j, k, irhoe) = w(i, j, k, irhoe) + &
&               twothird*w(ip, jp, kp, irho)*w(ip, jp, kp, itu1)
          end do
        end do
      end do
    end select
  end subroutine computepressure
 
!  differentiation of computelamviscosity in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: muref tref rgas *p *w *rlv
!   with respect to varying inputs: muref tref rgas *p *w *rlv
!   rw status of diff variables: muref:incr tref:incr rgas:incr
!                *p:incr *w:incr *rlv:in-out
!   plus diff mem management of: p:in w:in rlv:in
  subroutine computelamviscosity_b(includehalos)
!
!       computelamviscosity computes the laminar viscosity ratio in
!       the owned cell centers of the given block. sutherland's law is
!       used. it is assumed that the pointes already point to the
!       correct block before entering this subroutine.
!
    use blockpointers
    use constants
    use flowvarrefstate
    use inputphysics
    use iteration
    use utils_b, only : getcorrectfork
    implicit none
! input variables
    logical, intent(in) :: includehalos
!
!      local parameter.
!
    real(kind=realtype), parameter :: twothird=two*third
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii
    real(kind=realtype) :: musuth, tsuth, ssuth, t, pp
    real(kind=realtype) :: musuthd, tsuthd, ssuthd, td, ppd
    logical :: correctfork
    integer(kind=inttype) :: ibeg, iend, isize, jbeg, jend, jsize, kbeg&
&   , kend, ksize
    intrinsic mod
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype) :: tempd
    real(kind=realtype) :: tempd0
    real(kind=realtype) :: temp1
    real(kind=realtype) :: tempd1
! return immediately if no laminar viscosity needs to be computed.
    if (viscous) then
! determine whether or not the pressure must be corrected
! for the presence of the turbulent kinetic energy.
      correctfork = getcorrectfork()
! compute the nondimensional constants in sutherland's law.
      musuth = musuthdim/muref
      tsuth = tsuthdim/tref
      ssuth = ssuthdim/tref
      if (includehalos) then
        ibeg = 1
        jbeg = 1
        kbeg = 1
        iend = ie
        jend = je
        kend = ke
      else
        ibeg = 2
        jbeg = 2
        kbeg = 2
        iend = il
        jend = jl
        kend = kl
      end if
! substract 2/3 rho k, which is a part of the normal turbulent
! stresses, in case the pressure must be corrected.
      if (correctfork) then
        isize = iend - ibeg + 1
        jsize = jend - jbeg + 1
        ksize = kend - kbeg + 1
        ssuthd = 0.0_8
        musuthd = 0.0_8
        tsuthd = 0.0_8
!$bwd-of ii-loop 
        do ii=0,isize*jsize*ksize-1
          i = mod(ii, isize) + ibeg
          j = mod(ii/isize, jsize) + jbeg
          k = ii/(isize*jsize) + kbeg
          pp = p(i, j, k) - twothird*w(i, j, k, irho)*w(i, j, k, itu1)
          t = pp/(rgas*w(i, j, k, irho))
          temp0 = t/tsuth
          temp = musuth*(tsuth+ssuth)/(t+ssuth)
          tempd0 = temp0**1.5_realtype*rlvd(i, j, k)/(t+ssuth)
          tempd = 1.5_realtype*temp0**0.5*temp*rlvd(i, j, k)/tsuth
          rlvd(i, j, k) = 0.0_8
          tsuthd = tsuthd + musuth*tempd0 - temp0*tempd
          musuthd = musuthd + (tsuth+ssuth)*tempd0
          tempd1 = -(temp*tempd0)
          td = tempd + tempd1
          ssuthd = ssuthd + musuth*tempd0 + tempd1
          temp = w(i, j, k, irho)
          temp0 = rgas*temp
          ppd = td/temp0
          tempd = -(pp*td/temp0**2)
          rgasd = rgasd + temp*tempd
          wd(i, j, k, irho) = wd(i, j, k, irho) + rgas*tempd - w(i, j, k&
&           , itu1)*twothird*ppd
          pd(i, j, k) = pd(i, j, k) + ppd
          wd(i, j, k, itu1) = wd(i, j, k, itu1) - w(i, j, k, irho)*&
&           twothird*ppd
        end do
      else
! loop over the owned cells *and* first level halos of this
! block and compute the laminar viscosity ratio.
        isize = iend - ibeg + 1
        jsize = jend - jbeg + 1
        ksize = kend - kbeg + 1
        ssuthd = 0.0_8
        musuthd = 0.0_8
        tsuthd = 0.0_8
!$bwd-of ii-loop 
        do ii=0,isize*jsize*ksize-1
          i = mod(ii, isize) + ibeg
          j = mod(ii/isize, jsize) + jbeg
          k = ii/(isize*jsize) + kbeg
! compute the nondimensional temperature and the
! nondimensional laminar viscosity.
          t = p(i, j, k)/(rgas*w(i, j, k, irho))
          temp1 = t/tsuth
          temp0 = musuth*(tsuth+ssuth)/(t+ssuth)
          tempd = temp1**1.5_realtype*rlvd(i, j, k)/(t+ssuth)
          tempd1 = 1.5_realtype*temp1**0.5*temp0*rlvd(i, j, k)/tsuth
          rlvd(i, j, k) = 0.0_8
          tsuthd = tsuthd + musuth*tempd - temp1*tempd1
          musuthd = musuthd + (tsuth+ssuth)*tempd
          tempd0 = -(temp0*tempd)
          td = tempd1 + tempd0
          ssuthd = ssuthd + musuth*tempd + tempd0
          temp1 = w(i, j, k, irho)
          temp0 = rgas*temp1
          pd(i, j, k) = pd(i, j, k) + td/temp0
          tempd = -(p(i, j, k)*td/temp0**2)
          rgasd = rgasd + temp1*tempd
          wd(i, j, k, irho) = wd(i, j, k, irho) + rgas*tempd
        end do
      end if
      trefd = trefd - ssuthdim*ssuthd/tref**2 - tsuthdim*tsuthd/tref**2
      murefd = murefd - musuthdim*musuthd/muref**2
    end if
  end subroutine computelamviscosity_b
 
  subroutine computelamviscosity(includehalos)
!
!       computelamviscosity computes the laminar viscosity ratio in
!       the owned cell centers of the given block. sutherland's law is
!       used. it is assumed that the pointes already point to the
!       correct block before entering this subroutine.
!
    use blockpointers
    use constants
    use flowvarrefstate
    use inputphysics
    use iteration
    use utils_b, only : getcorrectfork
    implicit none
! input variables
    logical, intent(in) :: includehalos
!
!      local parameter.
!
    real(kind=realtype), parameter :: twothird=two*third
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii
    real(kind=realtype) :: musuth, tsuth, ssuth, t, pp
    logical :: correctfork
    integer(kind=inttype) :: ibeg, iend, isize, jbeg, jend, jsize, kbeg&
&   , kend, ksize
    intrinsic mod
! return immediately if no laminar viscosity needs to be computed.
    if (.not.viscous) then
      return
    else
! determine whether or not the pressure must be corrected
! for the presence of the turbulent kinetic energy.
      correctfork = getcorrectfork()
! compute the nondimensional constants in sutherland's law.
      musuth = musuthdim/muref
      tsuth = tsuthdim/tref
      ssuth = ssuthdim/tref
      if (includehalos) then
        ibeg = 1
        jbeg = 1
        kbeg = 1
        iend = ie
        jend = je
        kend = ke
      else
        ibeg = 2
        jbeg = 2
        kbeg = 2
        iend = il
        jend = jl
        kend = kl
      end if
! substract 2/3 rho k, which is a part of the normal turbulent
! stresses, in case the pressure must be corrected.
      if (correctfork) then
        isize = iend - ibeg + 1
        jsize = jend - jbeg + 1
        ksize = kend - kbeg + 1
!$ad ii-loop
        do ii=0,isize*jsize*ksize-1
          i = mod(ii, isize) + ibeg
          j = mod(ii/isize, jsize) + jbeg
          k = ii/(isize*jsize) + kbeg
          pp = p(i, j, k) - twothird*w(i, j, k, irho)*w(i, j, k, itu1)
          t = pp/(rgas*w(i, j, k, irho))
          rlv(i, j, k) = musuth*((tsuth+ssuth)/(t+ssuth))*(t/tsuth)**&
&           1.5_realtype
        end do
      else
! loop over the owned cells *and* first level halos of this
! block and compute the laminar viscosity ratio.
        isize = iend - ibeg + 1
        jsize = jend - jbeg + 1
        ksize = kend - kbeg + 1
!$ad ii-loop
        do ii=0,isize*jsize*ksize-1
          i = mod(ii, isize) + ibeg
          j = mod(ii/isize, jsize) + jbeg
          k = ii/(isize*jsize) + kbeg
! compute the nondimensional temperature and the
! nondimensional laminar viscosity.
          t = p(i, j, k)/(rgas*w(i, j, k, irho))
          rlv(i, j, k) = musuth*((tsuth+ssuth)/(t+ssuth))*(t/tsuth)**&
&           1.5_realtype
        end do
      end if
    end if
  end subroutine computelamviscosity
 
!  differentiation of adjustinflowangle in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: veldirfreestream dragdirection
!                liftdirection
!   with respect to varying inputs: alpha veldirfreestream beta
!                dragdirection liftdirection
!   rw status of diff variables: alpha:out veldirfreestream:in-out
!                beta:out dragdirection:in-out liftdirection:in-out
  subroutine adjustinflowangle_b()
    use constants
    use inputphysics, only : alpha, alphad, beta, betad, liftindex, &
&   veldirfreestream, veldirfreestreamd, liftdirection, liftdirectiond, &
&   dragdirection, dragdirectiond
    implicit none
!local vars
    real(kind=realtype), dimension(3) :: refdir1, refdir2
! velocity direction given by the rotation of a unit vector
! initially aligned along the positive x-direction (1,0,0)
! 1) rotate alpha radians cw about y or z-axis
! 2) rotate beta radians ccw about z or y-axis
    refdir1(:) = zero
    refdir1(1) = one
! drag direction given by the rotation of a unit vector
! initially aligned along the positive x-direction (1,0,0)
! 1) rotate alpha radians cw about y or z-axis
! 2) rotate beta radians ccw about z or y-axis
! lift direction given by the rotation of a unit vector
! initially aligned along the positive z-direction (0,0,1)
! 1) rotate alpha radians cw about y or z-axis
! 2) rotate beta radians ccw about z or y-axis
    refdir2(:) = zero
    refdir2(liftindex) = one
    alphad = 0.0_8
    betad = 0.0_8
    call getdirvector_b(refdir2, alpha, alphad, beta, betad, &
&                 liftdirection, liftdirectiond, liftindex)
    call getdirvector_b(refdir1, alpha, alphad, beta, betad, &
&                 dragdirection, dragdirectiond, liftindex)
    call getdirvector_b(refdir1, alpha, alphad, beta, betad, &
&                 veldirfreestream, veldirfreestreamd, liftindex)
  end subroutine adjustinflowangle_b
 
  subroutine adjustinflowangle()
    use constants
    use inputphysics, only : alpha, beta, liftindex, veldirfreestream,&
&   liftdirection, dragdirection
    implicit none
!local vars
    real(kind=realtype), dimension(3) :: refdir1, refdir2
! velocity direction given by the rotation of a unit vector
! initially aligned along the positive x-direction (1,0,0)
! 1) rotate alpha radians cw about y or z-axis
! 2) rotate beta radians ccw about z or y-axis
    refdir1(:) = zero
    refdir1(1) = one
    call getdirvector(refdir1, alpha, beta, veldirfreestream, liftindex)
! drag direction given by the rotation of a unit vector
! initially aligned along the positive x-direction (1,0,0)
! 1) rotate alpha radians cw about y or z-axis
! 2) rotate beta radians ccw about z or y-axis
    call getdirvector(refdir1, alpha, beta, dragdirection, liftindex)
! lift direction given by the rotation of a unit vector
! initially aligned along the positive z-direction (0,0,1)
! 1) rotate alpha radians cw about y or z-axis
! 2) rotate beta radians ccw about z or y-axis
    refdir2(:) = zero
    refdir2(liftindex) = one
    call getdirvector(refdir2, alpha, beta, liftdirection, liftindex)
  end subroutine adjustinflowangle
 
!  differentiation of derivativerotmatrixrigid in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: timeref rotationmatrix
!   with respect to varying inputs: timeref
  subroutine derivativerotmatrixrigid_b(rotationmatrix, rotationmatrixd&
&   , rotationpoint, t)
!
!       derivativerotmatrixrigid determines the derivative of the
!       rotation matrix at the given time for the rigid body rotation,
!       such that the grid velocities can be determined analytically.
!       also the rotation point of the current time level is
!       determined. this value can change due to translation of the
!       entire grid.
!
    use constants
    use flowvarrefstate
    use inputmotion
    use monitor
    use utils_b, only : rigidrotangle, derivativerigidrotangle, &
&   derivativerigidrotangle_b
    implicit none
!
!      subroutine arguments.
!
    real(kind=realtype), intent(in) :: t
    real(kind=realtype), dimension(3) :: rotationpoint
    real(kind=realtype), dimension(3, 3) :: rotationmatrix
    real(kind=realtype), dimension(3, 3) :: rotationmatrixd
!
!      local variables.
!
    integer(kind=inttype) :: i, j
    real(kind=realtype) :: phi, dphix, dphiy, dphiz
    real(kind=realtype) :: dphixd, dphiyd, dphizd
    real(kind=realtype) :: cosx, cosy, cosz, sinx, siny, sinz
    real(kind=realtype), dimension(3, 3) :: dm, m
    real(kind=realtype), dimension(3, 3) :: dmd
    intrinsic sin
    intrinsic cos
! determine the rotation angle around the x-axis for the new
! time level and the corresponding values of the sine and cosine.
    phi = rigidrotangle(degreepolxrot, coefpolxrot, degreefourxrot, &
&     omegafourxrot, coscoeffourxrot, sincoeffourxrot, t)
    sinx = sin(phi)
    cosx = cos(phi)
! idem for the y-axis.
    phi = rigidrotangle(degreepolyrot, coefpolyrot, degreefouryrot, &
&     omegafouryrot, coscoeffouryrot, sincoeffouryrot, t)
    siny = sin(phi)
    cosy = cos(phi)
! idem for the z-axis.
    phi = rigidrotangle(degreepolzrot, coefpolzrot, degreefourzrot, &
&     omegafourzrot, coscoeffourzrot, sincoeffourzrot, t)
    sinz = sin(phi)
    cosz = cos(phi)
! compute the time derivative of the rotation angles around the
! x-axis, y-axis and z-axis.
! compute the time derivative of the rotation matrix applied to
! the coordinates at t == 0.
! part 1. derivative of the z-rotation matrix multiplied by the
! x and y rotation matrix, i.e. dmz * my * mx
! part 2: mz * dmy * mx.
! part 3: mz * my * dmx
! determine the rotation matrix at t == t.
    m(1, 1) = cosy*cosz
    m(2, 1) = cosy*sinz
    m(3, 1) = -siny
    m(1, 2) = sinx*siny*cosz - cosx*sinz
    m(2, 2) = sinx*siny*sinz + cosx*cosz
    m(3, 2) = sinx*cosy
    m(1, 3) = cosx*siny*cosz + sinx*sinz
    m(2, 3) = cosx*siny*sinz - sinx*cosz
    m(3, 3) = cosx*cosy
    dmd = 0.0_8
    do j=3,1,-1
      do i=3,1,-1
        dmd(i, 1) = dmd(i, 1) + m(j, 1)*rotationmatrixd(i, j)
        dmd(i, 2) = dmd(i, 2) + m(j, 2)*rotationmatrixd(i, j)
        dmd(i, 3) = dmd(i, 3) + m(j, 3)*rotationmatrixd(i, j)
        rotationmatrixd(i, j) = 0.0_8
      end do
    end do
    dphixd = cosx*cosy*dmd(3, 2) - sinx*cosy*dmd(3, 3) + (cosx*siny*sinz&
&     -sinx*cosz)*dmd(2, 2) - (sinx*siny*sinz+cosx*cosz)*dmd(2, 3) + (&
&     cosx*sinz-sinx*siny*cosz)*dmd(1, 3) + (sinx*sinz+cosx*siny*cosz)*&
&     dmd(1, 2)
    dphiyd = sinz*cosx*cosy*dmd(2, 3) - cosx*siny*dmd(3, 3) - sinx*siny*&
&     dmd(3, 2) - cosy*dmd(3, 1) + sinz*sinx*cosy*dmd(2, 2) + cosz*cosx*&
&     cosy*dmd(1, 3) - siny*sinz*dmd(2, 1) + cosz*sinx*cosy*dmd(1, 2) - &
&     siny*cosz*dmd(1, 1)
    dmd(3, 3) = 0.0_8
    dmd(3, 2) = 0.0_8
    dmd(3, 1) = 0.0_8
    dphizd = (cosx*siny*cosz+sinx*sinz)*dmd(2, 3) + (sinx*siny*cosz-cosx&
&     *sinz)*dmd(2, 2) + cosy*cosz*dmd(2, 1) + (sinx*cosz-cosx*siny*sinz&
&     )*dmd(1, 3) + (-(cosx*cosz)-sinx*siny*sinz)*dmd(1, 2) - cosy*sinz*&
&     dmd(1, 1)
    dmd(2, 3) = 0.0_8
    dmd(2, 2) = 0.0_8
    dmd(2, 1) = 0.0_8
    dmd(1, 3) = 0.0_8
    dmd(1, 2) = 0.0_8
    call derivativerigidrotangle_b(degreepolzrot, coefpolzrot, &
&                            degreefourzrot, omegafourzrot, &
&                            coscoeffourzrot, sincoeffourzrot, t, dphizd&
&                           )
    call derivativerigidrotangle_b(degreepolyrot, coefpolyrot, &
&                            degreefouryrot, omegafouryrot, &
&                            coscoeffouryrot, sincoeffouryrot, t, dphiyd&
&                           )
    call derivativerigidrotangle_b(degreepolxrot, coefpolxrot, &
&                            degreefourxrot, omegafourxrot, &
&                            coscoeffourxrot, sincoeffourxrot, t, dphixd&
&                           )
  end subroutine derivativerotmatrixrigid_b
 
  subroutine derivativerotmatrixrigid(rotationmatrix, rotationpoint, t)
!
!       derivativerotmatrixrigid determines the derivative of the
!       rotation matrix at the given time for the rigid body rotation,
!       such that the grid velocities can be determined analytically.
!       also the rotation point of the current time level is
!       determined. this value can change due to translation of the
!       entire grid.
!
    use constants
    use flowvarrefstate
    use inputmotion
    use monitor
    use utils_b, only : rigidrotangle, derivativerigidrotangle
    implicit none
!
!      subroutine arguments.
!
    real(kind=realtype), intent(in) :: t
    real(kind=realtype), dimension(3), intent(out) :: rotationpoint
    real(kind=realtype), dimension(3, 3), intent(out) :: rotationmatrix
!
!      local variables.
!
    integer(kind=inttype) :: i, j
    real(kind=realtype) :: phi, dphix, dphiy, dphiz
    real(kind=realtype) :: cosx, cosy, cosz, sinx, siny, sinz
    real(kind=realtype), dimension(3, 3) :: dm, m
    intrinsic sin
    intrinsic cos
! determine the rotation angle around the x-axis for the new
! time level and the corresponding values of the sine and cosine.
    phi = rigidrotangle(degreepolxrot, coefpolxrot, degreefourxrot, &
&     omegafourxrot, coscoeffourxrot, sincoeffourxrot, t)
    sinx = sin(phi)
    cosx = cos(phi)
! idem for the y-axis.
    phi = rigidrotangle(degreepolyrot, coefpolyrot, degreefouryrot, &
&     omegafouryrot, coscoeffouryrot, sincoeffouryrot, t)
    siny = sin(phi)
    cosy = cos(phi)
! idem for the z-axis.
    phi = rigidrotangle(degreepolzrot, coefpolzrot, degreefourzrot, &
&     omegafourzrot, coscoeffourzrot, sincoeffourzrot, t)
    sinz = sin(phi)
    cosz = cos(phi)
! compute the time derivative of the rotation angles around the
! x-axis, y-axis and z-axis.
    dphix = derivativerigidrotangle(degreepolxrot, coefpolxrot, &
&     degreefourxrot, omegafourxrot, coscoeffourxrot, sincoeffourxrot, t&
&     )
    dphiy = derivativerigidrotangle(degreepolyrot, coefpolyrot, &
&     degreefouryrot, omegafouryrot, coscoeffouryrot, sincoeffouryrot, t&
&     )
    dphiz = derivativerigidrotangle(degreepolzrot, coefpolzrot, &
&     degreefourzrot, omegafourzrot, coscoeffourzrot, sincoeffourzrot, t&
&     )
! compute the time derivative of the rotation matrix applied to
! the coordinates at t == 0.
! part 1. derivative of the z-rotation matrix multiplied by the
! x and y rotation matrix, i.e. dmz * my * mx
    dm(1, 1) = -(cosy*sinz*dphiz)
    dm(1, 2) = (-(cosx*cosz)-sinx*siny*sinz)*dphiz
    dm(1, 3) = (sinx*cosz-cosx*siny*sinz)*dphiz
    dm(2, 1) = cosy*cosz*dphiz
    dm(2, 2) = (sinx*siny*cosz-cosx*sinz)*dphiz
    dm(2, 3) = (cosx*siny*cosz+sinx*sinz)*dphiz
    dm(3, 1) = zero
    dm(3, 2) = zero
    dm(3, 3) = zero
! part 2: mz * dmy * mx.
    dm(1, 1) = dm(1, 1) - siny*cosz*dphiy
    dm(1, 2) = dm(1, 2) + sinx*cosy*cosz*dphiy
    dm(1, 3) = dm(1, 3) + cosx*cosy*cosz*dphiy
    dm(2, 1) = dm(2, 1) - siny*sinz*dphiy
    dm(2, 2) = dm(2, 2) + sinx*cosy*sinz*dphiy
    dm(2, 3) = dm(2, 3) + cosx*cosy*sinz*dphiy
    dm(3, 1) = dm(3, 1) - cosy*dphiy
    dm(3, 2) = dm(3, 2) - sinx*siny*dphiy
    dm(3, 3) = dm(3, 3) - cosx*siny*dphiy
! part 3: mz * my * dmx
    dm(1, 2) = dm(1, 2) + (sinx*sinz+cosx*siny*cosz)*dphix
    dm(1, 3) = dm(1, 3) + (cosx*sinz-sinx*siny*cosz)*dphix
    dm(2, 2) = dm(2, 2) + (cosx*siny*sinz-sinx*cosz)*dphix
    dm(2, 3) = dm(2, 3) - (sinx*siny*sinz+cosx*cosz)*dphix
    dm(3, 2) = dm(3, 2) + cosx*cosy*dphix
    dm(3, 3) = dm(3, 3) - sinx*cosy*dphix
! determine the rotation matrix at t == t.
    m(1, 1) = cosy*cosz
    m(2, 1) = cosy*sinz
    m(3, 1) = -siny
    m(1, 2) = sinx*siny*cosz - cosx*sinz
    m(2, 2) = sinx*siny*sinz + cosx*cosz
    m(3, 2) = sinx*cosy
    m(1, 3) = cosx*siny*cosz + sinx*sinz
    m(2, 3) = cosx*siny*sinz - sinx*cosz
    m(3, 3) = cosx*cosy
! determine the matrix product dm * inverse(m), which will give
! the derivative of the rotation matrix when applied to the
! current coordinates. note that inverse(m) == transpose(m).
    do j=1,3
      do i=1,3
        rotationmatrix(i, j) = dm(i, 1)*m(j, 1) + dm(i, 2)*m(j, 2) + dm(&
&         i, 3)*m(j, 3)
      end do
    end do
! determine the rotation point at the new time level; it is
! possible that this value changes due to translation of the grid.
!  ainf = sqrt(gammainf*pinf/rhoinf)
!  rotationpoint(1) = lref*rotpoint(1) &
!                    + machgrid(1)*ainf*t/timeref
!  rotationpoint(2) = lref*rotpoint(2) &
!                    + machgrid(2)*ainf*t/timeref
!  rotationpoint(3) = lref*rotpoint(3) &
!                    + machgrid(3)*ainf*t/timeref
    rotationpoint(1) = lref*rotpoint(1)
    rotationpoint(2) = lref*rotpoint(2)
    rotationpoint(3) = lref*rotpoint(3)
  end subroutine derivativerotmatrixrigid
 
!  differentiation of getdirvector in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: alpha beta winddirection
!   with respect to varying inputs: alpha beta winddirection
  subroutine getdirvector_b(refdirection, alpha, alphad, beta, betad, &
&   winddirection, winddirectiond, liftindex)
!(xb,yb,zb,alpha,beta,xw,yw,zw)
!
!      convert the angle of attack and side slip angle to wind axes.
!      the components of the wind direction vector (xw,yw,zw) are
!      computed given the direction angles in radians and the body
!      direction by performing two rotations on the original
!      direction vector:
!        1) rotation about the zb or yb-axis: alpha clockwise (cw)
!           (xb,yb,zb) -> (x1,y1,z1)
!        2) rotation about the yl or z1-axis: beta counter-clockwise
!           (ccw)  (x1,y1,z1) -> (xw,yw,zw)
!         input arguments:
!            alpha    = angle of attack in radians
!            beta     = side slip angle in radians
!            refdirection = reference direction vector
!         output arguments:
!            winddirection = unit wind vector in body axes
!
    use constants
    use utils_b, only : terminate
    implicit none
!
!     subroutine arguments.
!
    real(kind=realtype), dimension(3), intent(in) :: refdirection
    real(kind=realtype) :: alpha, beta
    real(kind=realtype) :: alphad, betad
    real(kind=realtype), dimension(3) :: winddirection
    real(kind=realtype), dimension(3) :: winddirectiond
    integer(kind=inttype) :: liftindex
!
!     local variables.
!
    real(kind=realtype) :: rnorm, x1, y1, z1, xbn, ybn, zbn, xw, yw, zw
    real(kind=realtype) :: x1d, y1d, z1d, xbnd, ybnd, zbnd, xwd, ywd, &
&   zwd
    real(kind=realtype) :: tmp
    real(kind=realtype) :: tmpd
    intrinsic sqrt
    integer :: branch
! normalize the input vector.
    rnorm = sqrt(refdirection(1)**2 + refdirection(2)**2 + refdirection(&
&     3)**2)
    xbn = refdirection(1)/rnorm
    ybn = refdirection(2)/rnorm
    zbn = refdirection(3)/rnorm
!!$      ! compute the wind direction vector.
!!$
!!$      ! 1) rotate alpha radians cw about y-axis
!!$      !    ( <=> rotate y-axis alpha radians ccw)
!!$
!!$      call vectorrotation(x1, y1, z1, 2, alpha, xbn, ybn, zbn)
!!$
!!$      ! 2) rotate beta radians ccw about z-axis
!!$      !    ( <=> rotate z-axis -beta radians ccw)
!!$
!!$      call vectorrotation(xw, yw, zw, 3, -beta, x1, y1, z1)
    if (liftindex .eq. 2) then
! compute the wind direction vector.aerosurf axes different!!
! 1) rotate alpha radians cw about z-axis
!    ( <=> rotate z-axis alpha radians ccw)
      tmp = -alpha
      call vectorrotation(x1, y1, z1, 3, tmp, xbn, ybn, zbn)
! 2) rotate beta radians ccw about y-axis
!    ( <=> rotate z-axis -beta radians ccw)
      tmp = -beta
      call pushcontrol2b(0)
    else if (liftindex .eq. 3) then
! compute the wind direction vector.aerosurf axes different!!
! 1) rotate alpha radians cw about z-axis
!    ( <=> rotate z-axis alpha radians ccw)
      call vectorrotation(x1, y1, z1, 2, alpha, xbn, ybn, zbn)
! 2) rotate beta radians ccw about y-axis
!    ( <=> rotate z-axis -beta radians ccw)
      call pushcontrol2b(1)
    else
      call pushcontrol2b(2)
    end if
    zwd = winddirectiond(3)
    winddirectiond(3) = 0.0_8
    ywd = winddirectiond(2)
    winddirectiond(2) = 0.0_8
    xwd = winddirectiond(1)
    winddirectiond(1) = 0.0_8
    call popcontrol2b(branch)
    if (branch .eq. 0) then
      tmpd = 0.0_8
      call vectorrotation_b(xw, xwd, yw, ywd, zw, zwd, 2, tmp, tmpd, x1&
&                     , x1d, y1, y1d, z1, z1d)
      betad = betad - tmpd
      tmp = -alpha
      tmpd = 0.0_8
      call vectorrotation_b(x1, x1d, y1, y1d, z1, z1d, 3, tmp, tmpd, xbn&
&                     , xbnd, ybn, ybnd, zbn, zbnd)
      alphad = alphad - tmpd
    else if (branch .eq. 1) then
      call vectorrotation_b(xw, xwd, yw, ywd, zw, zwd, 3, beta, betad, &
&                     x1, x1d, y1, y1d, z1, z1d)
      call vectorrotation_b(x1, x1d, y1, y1d, z1, z1d, 2, alpha, alphad&
&                     , xbn, xbnd, ybn, ybnd, zbn, zbnd)
    end if
  end subroutine getdirvector_b
 
  subroutine getdirvector(refdirection, alpha, beta, winddirection, &
&   liftindex)
!(xb,yb,zb,alpha,beta,xw,yw,zw)
!
!      convert the angle of attack and side slip angle to wind axes.
!      the components of the wind direction vector (xw,yw,zw) are
!      computed given the direction angles in radians and the body
!      direction by performing two rotations on the original
!      direction vector:
!        1) rotation about the zb or yb-axis: alpha clockwise (cw)
!           (xb,yb,zb) -> (x1,y1,z1)
!        2) rotation about the yl or z1-axis: beta counter-clockwise
!           (ccw)  (x1,y1,z1) -> (xw,yw,zw)
!         input arguments:
!            alpha    = angle of attack in radians
!            beta     = side slip angle in radians
!            refdirection = reference direction vector
!         output arguments:
!            winddirection = unit wind vector in body axes
!
    use constants
    use utils_b, only : terminate
    implicit none
!
!     subroutine arguments.
!
    real(kind=realtype), dimension(3), intent(in) :: refdirection
    real(kind=realtype) :: alpha, beta
    real(kind=realtype), dimension(3), intent(out) :: winddirection
    integer(kind=inttype) :: liftindex
!
!     local variables.
!
    real(kind=realtype) :: rnorm, x1, y1, z1, xbn, ybn, zbn, xw, yw, zw
    real(kind=realtype) :: tmp
    intrinsic sqrt
! normalize the input vector.
    rnorm = sqrt(refdirection(1)**2 + refdirection(2)**2 + refdirection(&
&     3)**2)
    xbn = refdirection(1)/rnorm
    ybn = refdirection(2)/rnorm
    zbn = refdirection(3)/rnorm
!!$      ! compute the wind direction vector.
!!$
!!$      ! 1) rotate alpha radians cw about y-axis
!!$      !    ( <=> rotate y-axis alpha radians ccw)
!!$
!!$      call vectorrotation(x1, y1, z1, 2, alpha, xbn, ybn, zbn)
!!$
!!$      ! 2) rotate beta radians ccw about z-axis
!!$      !    ( <=> rotate z-axis -beta radians ccw)
!!$
!!$      call vectorrotation(xw, yw, zw, 3, -beta, x1, y1, z1)
    if (liftindex .eq. 2) then
! compute the wind direction vector.aerosurf axes different!!
! 1) rotate alpha radians cw about z-axis
!    ( <=> rotate z-axis alpha radians ccw)
      tmp = -alpha
      call vectorrotation(x1, y1, z1, 3, tmp, xbn, ybn, zbn)
! 2) rotate beta radians ccw about y-axis
!    ( <=> rotate z-axis -beta radians ccw)
      tmp = -beta
      call vectorrotation(xw, yw, zw, 2, tmp, x1, y1, z1)
    else if (liftindex .eq. 3) then
! compute the wind direction vector.aerosurf axes different!!
! 1) rotate alpha radians cw about z-axis
!    ( <=> rotate z-axis alpha radians ccw)
      call vectorrotation(x1, y1, z1, 2, alpha, xbn, ybn, zbn)
! 2) rotate beta radians ccw about y-axis
!    ( <=> rotate z-axis -beta radians ccw)
      call vectorrotation(xw, yw, zw, 3, beta, x1, y1, z1)
    else
      call terminate('getdirvector', 'invalid lift direction')
    end if
    winddirection(1) = xw
    winddirection(2) = yw
    winddirection(3) = zw
  end subroutine getdirvector
 
!  differentiation of vectorrotation in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: xp yp angle zp
!   with respect to varying inputs: x y z angle
  subroutine vectorrotation_b(xp, xpd, yp, ypd, zp, zpd, iaxis, angle, &
&   angled, x, xd, y, yd, z, zd)
!
!      vectorrotation rotates a given vector with respect to a
!      specified axis by a given angle.
!         input arguments:
!            iaxis      = rotation axis (1-x, 2-y, 3-z)
!            angle      = rotation angle (measured ccw in radians)
!            x, y, z    = coordinates in original system
!         output arguments:
!            xp, yp, zp = coordinates in rotated system
!
    use precision
    implicit none
!
!     subroutine arguments.
!
    integer(kind=inttype), intent(in) :: iaxis
    real(kind=realtype), intent(in) :: angle, x, y, z
    real(kind=realtype) :: angled, xd, yd, zd
    real(kind=realtype) :: xp, yp, zp
    real(kind=realtype) :: xpd, ypd, zpd
    intrinsic cos
    intrinsic sin
! rotation about specified axis by specified angle
    select case  (iaxis) 
    case (1) 
      xd = xpd
      angled = angled + (cos(angle)*z-sin(angle)*y)*ypd - (cos(angle)*y+&
&       sin(angle)*z)*zpd
      yd = cos(angle)*ypd - sin(angle)*zpd
      zd = cos(angle)*zpd + sin(angle)*ypd
    case (2) 
      angled = angled + (cos(angle)*x-sin(angle)*z)*zpd - (sin(angle)*x+&
&       cos(angle)*z)*xpd
      xd = sin(angle)*zpd + cos(angle)*xpd
      yd = ypd
      zd = cos(angle)*zpd - sin(angle)*xpd
    case (3) 
      xd = cos(angle)*xpd - sin(angle)*ypd
      yd = cos(angle)*ypd + sin(angle)*xpd
      zd = zpd
      angled = angled + (cos(angle)*y-sin(angle)*x)*xpd - (sin(angle)*y+&
&       cos(angle)*x)*ypd
    case default
      stop
    end select
  end subroutine vectorrotation_b
 
  subroutine vectorrotation(xp, yp, zp, iaxis, angle, x, y, z)
!
!      vectorrotation rotates a given vector with respect to a
!      specified axis by a given angle.
!         input arguments:
!            iaxis      = rotation axis (1-x, 2-y, 3-z)
!            angle      = rotation angle (measured ccw in radians)
!            x, y, z    = coordinates in original system
!         output arguments:
!            xp, yp, zp = coordinates in rotated system
!
    use precision
    implicit none
!
!     subroutine arguments.
!
    integer(kind=inttype), intent(in) :: iaxis
    real(kind=realtype), intent(in) :: angle, x, y, z
    real(kind=realtype), intent(out) :: xp, yp, zp
    intrinsic cos
    intrinsic sin
! rotation about specified axis by specified angle
    select case  (iaxis) 
    case (1) 
! rotation about the x-axis
      xp = 1.*x + 0.*y + 0.*z
      yp = 0.*x + cos(angle)*y + sin(angle)*z
      zp = 0.*x - sin(angle)*y + cos(angle)*z
! rotation about the y-axis
    case (2) 
      xp = cos(angle)*x + 0.*y - sin(angle)*z
      yp = 0.*x + 1.*y + 0.*z
      zp = sin(angle)*x + 0.*y + cos(angle)*z
! rotation about the z-axis
    case (3) 
      xp = cos(angle)*x + sin(angle)*y + 0.*z
      yp = -(sin(angle)*x) + cos(angle)*y + 0.*z
      zp = 0.*x + 0.*y + 1.*z
    case default
      write(*, *) 'vectorrotation called with invalid arguments'
      stop
    end select
  end subroutine vectorrotation
 
!  differentiation of allnodalgradients in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *aa *wx *wy *wz *w *qx *qy
!                *qz *vol *ux *uy *uz *si *sj *sk *vx *vy *vz
!   with respect to varying inputs: *aa *wx *wy *wz *w *qx *qy
!                *qz *vol *ux *uy *uz *si *sj *sk *vx *vy *vz
!   rw status of diff variables: *aa:incr *wx:in-out *wy:in-out
!                *wz:in-out *w:incr *qx:in-out *qy:in-out *qz:in-out
!                *vol:incr *ux:in-out *uy:in-out *uz:in-out *si:incr
!                *sj:incr *sk:incr *vx:in-out *vy:in-out *vz:in-out
!   plus diff mem management of: aa:in wx:in wy:in wz:in w:in qx:in
!                qy:in qz:in vol:in ux:in uy:in uz:in si:in sj:in
!                sk:in vx:in vy:in vz:in
  subroutine allnodalgradients_b()
!
!         nodalgradients computes the nodal velocity gradients and
!         minus the gradient of the speed of sound squared. the minus
!         sign is present, because this is the definition of the heat
!         flux. these gradients are computed for all nodes.
!
    use constants
    use blockpointers
    implicit none
!        local variables.
    integer(kind=inttype) :: i, j, k
    integer(kind=inttype) :: k1, kk
    integer(kind=inttype) :: istart, iend, isize, ii
    integer(kind=inttype) :: jstart, jend, jsize
    integer(kind=inttype) :: kstart, kend, ksize
    real(kind=realtype) :: oneoverv, ubar, vbar, wbar, a2
    real(kind=realtype) :: oneovervd, ubard, vbard, wbard, a2d
    real(kind=realtype) :: sx, sx1, sy, sy1, sz, sz1
    real(kind=realtype) :: sxd, syd, szd
    intrinsic mod
    real(kind=realtype) :: temp
    real(kind=realtype) :: tempd
    integer :: branch
! zero all nodeal gradients:
    ux = zero
    uy = zero
    uz = zero
    vx = zero
    vy = zero
    vz = zero
    wx = zero
    wy = zero
    wz = zero
    qx = zero
    qy = zero
    qz = zero
!$fwd-of ii-loop 
! first part. contribution in the k-direction.
! the contribution is scattered to both the left and right node
! in k-direction.
    do ii=0,il*jl*ke-1
      i = mod(ii, il) + 1
      j = mod(ii/il, jl) + 1
      k = ii/(il*jl) + 1
! compute 8 times the average normal for this part of
! the control volume. the factor 8 is taken care of later
! on when the division by the volume takes place.
      sx = sk(i, j, k-1, 1) + sk(i+1, j, k-1, 1) + sk(i, j+1, k-1, 1) + &
&       sk(i+1, j+1, k-1, 1) + sk(i, j, k, 1) + sk(i+1, j, k, 1) + sk(i&
&       , j+1, k, 1) + sk(i+1, j+1, k, 1)
      sy = sk(i, j, k-1, 2) + sk(i+1, j, k-1, 2) + sk(i, j+1, k-1, 2) + &
&       sk(i+1, j+1, k-1, 2) + sk(i, j, k, 2) + sk(i+1, j, k, 2) + sk(i&
&       , j+1, k, 2) + sk(i+1, j+1, k, 2)
      sz = sk(i, j, k-1, 3) + sk(i+1, j, k-1, 3) + sk(i, j+1, k-1, 3) + &
&       sk(i+1, j+1, k-1, 3) + sk(i, j, k, 3) + sk(i+1, j, k, 3) + sk(i&
&       , j+1, k, 3) + sk(i+1, j+1, k, 3)
! compute the average velocities and speed of sound squared
! for this integration point. node that these variables are
! stored in w(ivx), w(ivy), w(ivz) and p.
      ubar = fourth*(w(i, j, k, ivx)+w(i+1, j, k, ivx)+w(i, j+1, k, ivx)&
&       +w(i+1, j+1, k, ivx))
      vbar = fourth*(w(i, j, k, ivy)+w(i+1, j, k, ivy)+w(i, j+1, k, ivy)&
&       +w(i+1, j+1, k, ivy))
      wbar = fourth*(w(i, j, k, ivz)+w(i+1, j, k, ivz)+w(i, j+1, k, ivz)&
&       +w(i+1, j+1, k, ivz))
      a2 = fourth*(aa(i, j, k)+aa(i+1, j, k)+aa(i, j+1, k)+aa(i+1, j+1, &
&       k))
! add the contributions to the surface integral to the node
! j-1 and substract it from the node j. for the heat flux it
! is reversed, because the negative of the gradient of the
! speed of sound must be computed.
      if (k .gt. 1) then
        ux(i, j, k-1) = ux(i, j, k-1) + ubar*sx
        uy(i, j, k-1) = uy(i, j, k-1) + ubar*sy
        uz(i, j, k-1) = uz(i, j, k-1) + ubar*sz
        vx(i, j, k-1) = vx(i, j, k-1) + vbar*sx
        vy(i, j, k-1) = vy(i, j, k-1) + vbar*sy
        vz(i, j, k-1) = vz(i, j, k-1) + vbar*sz
        wx(i, j, k-1) = wx(i, j, k-1) + wbar*sx
        wy(i, j, k-1) = wy(i, j, k-1) + wbar*sy
        wz(i, j, k-1) = wz(i, j, k-1) + wbar*sz
        qx(i, j, k-1) = qx(i, j, k-1) - a2*sx
        qy(i, j, k-1) = qy(i, j, k-1) - a2*sy
        qz(i, j, k-1) = qz(i, j, k-1) - a2*sz
      end if
      if (k .lt. ke) then
        ux(i, j, k) = ux(i, j, k) - ubar*sx
        uy(i, j, k) = uy(i, j, k) - ubar*sy
        uz(i, j, k) = uz(i, j, k) - ubar*sz
        vx(i, j, k) = vx(i, j, k) - vbar*sx
        vy(i, j, k) = vy(i, j, k) - vbar*sy
        vz(i, j, k) = vz(i, j, k) - vbar*sz
        wx(i, j, k) = wx(i, j, k) - wbar*sx
        wy(i, j, k) = wy(i, j, k) - wbar*sy
        wz(i, j, k) = wz(i, j, k) - wbar*sz
        qx(i, j, k) = qx(i, j, k) + a2*sx
        qy(i, j, k) = qy(i, j, k) + a2*sy
        qz(i, j, k) = qz(i, j, k) + a2*sz
      end if
    end do
!$fwd-of ii-loop 
! second part. contribution in the j-direction.
! the contribution is scattered to both the left and right node
! in j-direction.
    do ii=0,il*je*kl-1
      i = mod(ii, il) + 1
      j = mod(ii/il, je) + 1
      k = ii/(il*je) + 1
! compute 8 times the average normal for this part of
! the control volume. the factor 8 is taken care of later
! on when the division by the volume takes place.
      sx = sj(i, j-1, k, 1) + sj(i+1, j-1, k, 1) + sj(i, j-1, k+1, 1) + &
&       sj(i+1, j-1, k+1, 1) + sj(i, j, k, 1) + sj(i+1, j, k, 1) + sj(i&
&       , j, k+1, 1) + sj(i+1, j, k+1, 1)
      sy = sj(i, j-1, k, 2) + sj(i+1, j-1, k, 2) + sj(i, j-1, k+1, 2) + &
&       sj(i+1, j-1, k+1, 2) + sj(i, j, k, 2) + sj(i+1, j, k, 2) + sj(i&
&       , j, k+1, 2) + sj(i+1, j, k+1, 2)
      sz = sj(i, j-1, k, 3) + sj(i+1, j-1, k, 3) + sj(i, j-1, k+1, 3) + &
&       sj(i+1, j-1, k+1, 3) + sj(i, j, k, 3) + sj(i+1, j, k, 3) + sj(i&
&       , j, k+1, 3) + sj(i+1, j, k+1, 3)
! compute the average velocities and speed of sound squared
! for this integration point. node that these variables are
! stored in w(ivx), w(ivy), w(ivz) and p.
      ubar = fourth*(w(i, j, k, ivx)+w(i+1, j, k, ivx)+w(i, j, k+1, ivx)&
&       +w(i+1, j, k+1, ivx))
      vbar = fourth*(w(i, j, k, ivy)+w(i+1, j, k, ivy)+w(i, j, k+1, ivy)&
&       +w(i+1, j, k+1, ivy))
      wbar = fourth*(w(i, j, k, ivz)+w(i+1, j, k, ivz)+w(i, j, k+1, ivz)&
&       +w(i+1, j, k+1, ivz))
      a2 = fourth*(aa(i, j, k)+aa(i+1, j, k)+aa(i, j, k+1)+aa(i+1, j, k+&
&       1))
! add the contributions to the surface integral to the node
! j-1 and substract it from the node j. for the heat flux it
! is reversed, because the negative of the gradient of the
! speed of sound must be computed.
      if (j .gt. 1) then
        ux(i, j-1, k) = ux(i, j-1, k) + ubar*sx
        uy(i, j-1, k) = uy(i, j-1, k) + ubar*sy
        uz(i, j-1, k) = uz(i, j-1, k) + ubar*sz
        vx(i, j-1, k) = vx(i, j-1, k) + vbar*sx
        vy(i, j-1, k) = vy(i, j-1, k) + vbar*sy
        vz(i, j-1, k) = vz(i, j-1, k) + vbar*sz
        wx(i, j-1, k) = wx(i, j-1, k) + wbar*sx
        wy(i, j-1, k) = wy(i, j-1, k) + wbar*sy
        wz(i, j-1, k) = wz(i, j-1, k) + wbar*sz
        qx(i, j-1, k) = qx(i, j-1, k) - a2*sx
        qy(i, j-1, k) = qy(i, j-1, k) - a2*sy
        qz(i, j-1, k) = qz(i, j-1, k) - a2*sz
      end if
      if (j .lt. je) then
        ux(i, j, k) = ux(i, j, k) - ubar*sx
        uy(i, j, k) = uy(i, j, k) - ubar*sy
        uz(i, j, k) = uz(i, j, k) - ubar*sz
        vx(i, j, k) = vx(i, j, k) - vbar*sx
        vy(i, j, k) = vy(i, j, k) - vbar*sy
        vz(i, j, k) = vz(i, j, k) - vbar*sz
        wx(i, j, k) = wx(i, j, k) - wbar*sx
        wy(i, j, k) = wy(i, j, k) - wbar*sy
        wz(i, j, k) = wz(i, j, k) - wbar*sz
        qx(i, j, k) = qx(i, j, k) + a2*sx
        qy(i, j, k) = qy(i, j, k) + a2*sy
        qz(i, j, k) = qz(i, j, k) + a2*sz
      end if
    end do
!$fwd-of ii-loop 
! third part. contribution in the i-direction.
! the contribution is scattered to both the left and right node
! in i-direction.
    do ii=0,ie*jl*kl-1
      i = mod(ii, ie) + 1
      j = mod(ii/ie, jl) + 1
      k = ii/(ie*jl) + 1
! compute 8 times the average normal for this part of
! the control volume. the factor 8 is taken care of later
! on when the division by the volume takes place.
      sx = si(i-1, j, k, 1) + si(i-1, j+1, k, 1) + si(i-1, j, k+1, 1) + &
&       si(i-1, j+1, k+1, 1) + si(i, j, k, 1) + si(i, j+1, k, 1) + si(i&
&       , j, k+1, 1) + si(i, j+1, k+1, 1)
      sy = si(i-1, j, k, 2) + si(i-1, j+1, k, 2) + si(i-1, j, k+1, 2) + &
&       si(i-1, j+1, k+1, 2) + si(i, j, k, 2) + si(i, j+1, k, 2) + si(i&
&       , j, k+1, 2) + si(i, j+1, k+1, 2)
      sz = si(i-1, j, k, 3) + si(i-1, j+1, k, 3) + si(i-1, j, k+1, 3) + &
&       si(i-1, j+1, k+1, 3) + si(i, j, k, 3) + si(i, j+1, k, 3) + si(i&
&       , j, k+1, 3) + si(i, j+1, k+1, 3)
! compute the average velocities and speed of sound squared
! for this integration point. node that these variables are
! stored in w(ivx), w(ivy), w(ivz) and p.
      ubar = fourth*(w(i, j, k, ivx)+w(i, j+1, k, ivx)+w(i, j, k+1, ivx)&
&       +w(i, j+1, k+1, ivx))
      vbar = fourth*(w(i, j, k, ivy)+w(i, j+1, k, ivy)+w(i, j, k+1, ivy)&
&       +w(i, j+1, k+1, ivy))
      wbar = fourth*(w(i, j, k, ivz)+w(i, j+1, k, ivz)+w(i, j, k+1, ivz)&
&       +w(i, j+1, k+1, ivz))
      a2 = fourth*(aa(i, j, k)+aa(i, j+1, k)+aa(i, j, k+1)+aa(i, j+1, k+&
&       1))
! add the contributions to the surface integral to the node
! j-1 and substract it from the node j. for the heat flux it
! is reversed, because the negative of the gradient of the
! speed of sound must be computed.
      if (i .gt. 1) then
        ux(i-1, j, k) = ux(i-1, j, k) + ubar*sx
        uy(i-1, j, k) = uy(i-1, j, k) + ubar*sy
        uz(i-1, j, k) = uz(i-1, j, k) + ubar*sz
        vx(i-1, j, k) = vx(i-1, j, k) + vbar*sx
        vy(i-1, j, k) = vy(i-1, j, k) + vbar*sy
        vz(i-1, j, k) = vz(i-1, j, k) + vbar*sz
        wx(i-1, j, k) = wx(i-1, j, k) + wbar*sx
        wy(i-1, j, k) = wy(i-1, j, k) + wbar*sy
        wz(i-1, j, k) = wz(i-1, j, k) + wbar*sz
        qx(i-1, j, k) = qx(i-1, j, k) - a2*sx
        qy(i-1, j, k) = qy(i-1, j, k) - a2*sy
        qz(i-1, j, k) = qz(i-1, j, k) - a2*sz
      end if
      if (i .lt. ie) then
        ux(i, j, k) = ux(i, j, k) - ubar*sx
        uy(i, j, k) = uy(i, j, k) - ubar*sy
        uz(i, j, k) = uz(i, j, k) - ubar*sz
        vx(i, j, k) = vx(i, j, k) - vbar*sx
        vy(i, j, k) = vy(i, j, k) - vbar*sy
        vz(i, j, k) = vz(i, j, k) - vbar*sz
        wx(i, j, k) = wx(i, j, k) - wbar*sx
        wy(i, j, k) = wy(i, j, k) - wbar*sy
        wz(i, j, k) = wz(i, j, k) - wbar*sz
        qx(i, j, k) = qx(i, j, k) + a2*sx
        qy(i, j, k) = qy(i, j, k) + a2*sy
        qz(i, j, k) = qz(i, j, k) + a2*sz
      end if
    end do
!$bwd-of ii-loop 
    do ii=0,il*jl*kl-1
      i = mod(ii, il) + 1
      j = mod(ii/il, jl) + 1
      k = ii/(il*jl) + 1
! compute the inverse of 8 times the volume for this node.
      oneoverv = one/(vol(i, j, k)+vol(i, j, k+1)+vol(i+1, j, k)+vol(i+1&
&       , j, k+1)+vol(i, j+1, k)+vol(i, j+1, k+1)+vol(i+1, j+1, k)+vol(i&
&       +1, j+1, k+1))
! compute the correct velocity gradients and "unit" heat
! fluxes. the velocity gradients are stored in ux, etc.
      oneovervd = qz(i, j, k)*qzd(i, j, k) + qy(i, j, k)*qyd(i, j, k) + &
&       qx(i, j, k)*qxd(i, j, k) + wz(i, j, k)*wzd(i, j, k) + wy(i, j, k&
&       )*wyd(i, j, k) + wx(i, j, k)*wxd(i, j, k) + vz(i, j, k)*vzd(i, j&
&       , k) + vy(i, j, k)*vyd(i, j, k) + vx(i, j, k)*vxd(i, j, k) + uz(&
&       i, j, k)*uzd(i, j, k) + uy(i, j, k)*uyd(i, j, k) + ux(i, j, k)*&
&       uxd(i, j, k)
      qzd(i, j, k) = oneoverv*qzd(i, j, k)
      qyd(i, j, k) = oneoverv*qyd(i, j, k)
      qxd(i, j, k) = oneoverv*qxd(i, j, k)
      wzd(i, j, k) = oneoverv*wzd(i, j, k)
      wyd(i, j, k) = oneoverv*wyd(i, j, k)
      wxd(i, j, k) = oneoverv*wxd(i, j, k)
      vzd(i, j, k) = oneoverv*vzd(i, j, k)
      vyd(i, j, k) = oneoverv*vyd(i, j, k)
      vxd(i, j, k) = oneoverv*vxd(i, j, k)
      uzd(i, j, k) = oneoverv*uzd(i, j, k)
      uyd(i, j, k) = oneoverv*uyd(i, j, k)
      uxd(i, j, k) = oneoverv*uxd(i, j, k)
      temp = vol(i, j, k) + vol(i, j, k+1) + vol(i+1, j, k) + vol(i+1, j&
&       , k+1) + vol(i, j+1, k) + vol(i, j+1, k+1) + vol(i+1, j+1, k) + &
&       vol(i+1, j+1, k+1)
      tempd = -(one*oneovervd/temp**2)
      vold(i, j, k) = vold(i, j, k) + tempd
      vold(i, j, k+1) = vold(i, j, k+1) + tempd
      vold(i+1, j, k) = vold(i+1, j, k) + tempd
      vold(i+1, j, k+1) = vold(i+1, j, k+1) + tempd
      vold(i, j+1, k) = vold(i, j+1, k) + tempd
      vold(i, j+1, k+1) = vold(i, j+1, k+1) + tempd
      vold(i+1, j+1, k) = vold(i+1, j+1, k) + tempd
      vold(i+1, j+1, k+1) = vold(i+1, j+1, k+1) + tempd
    end do
!$bwd-of ii-loop 
    do ii=0,ie*jl*kl-1
      i = mod(ii, ie) + 1
      j = mod(ii/ie, jl) + 1
      k = ii/(ie*jl) + 1
! compute 8 times the average normal for this part of
! the control volume. the factor 8 is taken care of later
! on when the division by the volume takes place.
      sx = si(i-1, j, k, 1) + si(i-1, j+1, k, 1) + si(i-1, j, k+1, 1) + &
&       si(i-1, j+1, k+1, 1) + si(i, j, k, 1) + si(i, j+1, k, 1) + si(i&
&       , j, k+1, 1) + si(i, j+1, k+1, 1)
      sy = si(i-1, j, k, 2) + si(i-1, j+1, k, 2) + si(i-1, j, k+1, 2) + &
&       si(i-1, j+1, k+1, 2) + si(i, j, k, 2) + si(i, j+1, k, 2) + si(i&
&       , j, k+1, 2) + si(i, j+1, k+1, 2)
      sz = si(i-1, j, k, 3) + si(i-1, j+1, k, 3) + si(i-1, j, k+1, 3) + &
&       si(i-1, j+1, k+1, 3) + si(i, j, k, 3) + si(i, j+1, k, 3) + si(i&
&       , j, k+1, 3) + si(i, j+1, k+1, 3)
! compute the average velocities and speed of sound squared
! for this integration point. node that these variables are
! stored in w(ivx), w(ivy), w(ivz) and p.
      ubar = fourth*(w(i, j, k, ivx)+w(i, j+1, k, ivx)+w(i, j, k+1, ivx)&
&       +w(i, j+1, k+1, ivx))
      vbar = fourth*(w(i, j, k, ivy)+w(i, j+1, k, ivy)+w(i, j, k+1, ivy)&
&       +w(i, j+1, k+1, ivy))
      wbar = fourth*(w(i, j, k, ivz)+w(i, j+1, k, ivz)+w(i, j, k+1, ivz)&
&       +w(i, j+1, k+1, ivz))
      a2 = fourth*(aa(i, j, k)+aa(i, j+1, k)+aa(i, j, k+1)+aa(i, j+1, k+&
&       1))
! add the contributions to the surface integral to the node
! j-1 and substract it from the node j. for the heat flux it
! is reversed, because the negative of the gradient of the
! speed of sound must be computed.
      if (i .gt. 1) then
        call pushcontrol1b(0)
      else
        call pushcontrol1b(1)
      end if
      if (i .lt. ie) then
        a2d = sz*qzd(i, j, k) + sy*qyd(i, j, k) + sx*qxd(i, j, k)
        szd = a2*qzd(i, j, k) - wbar*wzd(i, j, k) - vbar*vzd(i, j, k) - &
&         ubar*uzd(i, j, k)
        syd = a2*qyd(i, j, k) - wbar*wyd(i, j, k) - vbar*vyd(i, j, k) - &
&         ubar*uyd(i, j, k)
        sxd = a2*qxd(i, j, k) - wbar*wxd(i, j, k) - vbar*vxd(i, j, k) - &
&         ubar*uxd(i, j, k)
        wbard = -(sz*wzd(i, j, k)) - sy*wyd(i, j, k) - sx*wxd(i, j, k)
        vbard = -(sz*vzd(i, j, k)) - sy*vyd(i, j, k) - sx*vxd(i, j, k)
        ubard = -(sz*uzd(i, j, k)) - sy*uyd(i, j, k) - sx*uxd(i, j, k)
      else
        wbard = 0.0_8
        vbard = 0.0_8
        ubard = 0.0_8
        a2d = 0.0_8
        sxd = 0.0_8
        syd = 0.0_8
        szd = 0.0_8
      end if
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        a2d = a2d - sz*qzd(i-1, j, k) - sy*qyd(i-1, j, k) - sx*qxd(i-1, &
&         j, k)
        szd = szd + wbar*wzd(i-1, j, k) - a2*qzd(i-1, j, k) + vbar*vzd(i&
&         -1, j, k) + ubar*uzd(i-1, j, k)
        syd = syd + wbar*wyd(i-1, j, k) - a2*qyd(i-1, j, k) + vbar*vyd(i&
&         -1, j, k) + ubar*uyd(i-1, j, k)
        sxd = sxd + wbar*wxd(i-1, j, k) - a2*qxd(i-1, j, k) + vbar*vxd(i&
&         -1, j, k) + ubar*uxd(i-1, j, k)
        wbard = wbard + sz*wzd(i-1, j, k) + sy*wyd(i-1, j, k) + sx*wxd(i&
&         -1, j, k)
        vbard = vbard + sz*vzd(i-1, j, k) + sy*vyd(i-1, j, k) + sx*vxd(i&
&         -1, j, k)
        ubard = ubard + sz*uzd(i-1, j, k) + sy*uyd(i-1, j, k) + sx*uxd(i&
&         -1, j, k)
      end if
      tempd = fourth*a2d
      aad(i, j, k) = aad(i, j, k) + tempd
      aad(i, j+1, k) = aad(i, j+1, k) + tempd
      aad(i, j, k+1) = aad(i, j, k+1) + tempd
      aad(i, j+1, k+1) = aad(i, j+1, k+1) + tempd
      tempd = fourth*wbard
      wd(i, j, k, ivz) = wd(i, j, k, ivz) + tempd
      wd(i, j+1, k, ivz) = wd(i, j+1, k, ivz) + tempd
      wd(i, j, k+1, ivz) = wd(i, j, k+1, ivz) + tempd
      wd(i, j+1, k+1, ivz) = wd(i, j+1, k+1, ivz) + tempd
      tempd = fourth*vbard
      wd(i, j, k, ivy) = wd(i, j, k, ivy) + tempd
      wd(i, j+1, k, ivy) = wd(i, j+1, k, ivy) + tempd
      wd(i, j, k+1, ivy) = wd(i, j, k+1, ivy) + tempd
      wd(i, j+1, k+1, ivy) = wd(i, j+1, k+1, ivy) + tempd
      tempd = fourth*ubard
      wd(i, j, k, ivx) = wd(i, j, k, ivx) + tempd
      wd(i, j+1, k, ivx) = wd(i, j+1, k, ivx) + tempd
      wd(i, j, k+1, ivx) = wd(i, j, k+1, ivx) + tempd
      wd(i, j+1, k+1, ivx) = wd(i, j+1, k+1, ivx) + tempd
      sid(i-1, j, k, 3) = sid(i-1, j, k, 3) + szd
      sid(i-1, j+1, k, 3) = sid(i-1, j+1, k, 3) + szd
      sid(i-1, j, k+1, 3) = sid(i-1, j, k+1, 3) + szd
      sid(i-1, j+1, k+1, 3) = sid(i-1, j+1, k+1, 3) + szd
      sid(i, j, k, 3) = sid(i, j, k, 3) + szd
      sid(i, j+1, k, 3) = sid(i, j+1, k, 3) + szd
      sid(i, j, k+1, 3) = sid(i, j, k+1, 3) + szd
      sid(i, j+1, k+1, 3) = sid(i, j+1, k+1, 3) + szd
      sid(i-1, j, k, 2) = sid(i-1, j, k, 2) + syd
      sid(i-1, j+1, k, 2) = sid(i-1, j+1, k, 2) + syd
      sid(i-1, j, k+1, 2) = sid(i-1, j, k+1, 2) + syd
      sid(i-1, j+1, k+1, 2) = sid(i-1, j+1, k+1, 2) + syd
      sid(i, j, k, 2) = sid(i, j, k, 2) + syd
      sid(i, j+1, k, 2) = sid(i, j+1, k, 2) + syd
      sid(i, j, k+1, 2) = sid(i, j, k+1, 2) + syd
      sid(i, j+1, k+1, 2) = sid(i, j+1, k+1, 2) + syd
      sid(i-1, j, k, 1) = sid(i-1, j, k, 1) + sxd
      sid(i-1, j+1, k, 1) = sid(i-1, j+1, k, 1) + sxd
      sid(i-1, j, k+1, 1) = sid(i-1, j, k+1, 1) + sxd
      sid(i-1, j+1, k+1, 1) = sid(i-1, j+1, k+1, 1) + sxd
      sid(i, j, k, 1) = sid(i, j, k, 1) + sxd
      sid(i, j+1, k, 1) = sid(i, j+1, k, 1) + sxd
      sid(i, j, k+1, 1) = sid(i, j, k+1, 1) + sxd
      sid(i, j+1, k+1, 1) = sid(i, j+1, k+1, 1) + sxd
    end do
!$bwd-of ii-loop 
    do ii=0,il*je*kl-1
      i = mod(ii, il) + 1
      j = mod(ii/il, je) + 1
      k = ii/(il*je) + 1
! compute 8 times the average normal for this part of
! the control volume. the factor 8 is taken care of later
! on when the division by the volume takes place.
      sx = sj(i, j-1, k, 1) + sj(i+1, j-1, k, 1) + sj(i, j-1, k+1, 1) + &
&       sj(i+1, j-1, k+1, 1) + sj(i, j, k, 1) + sj(i+1, j, k, 1) + sj(i&
&       , j, k+1, 1) + sj(i+1, j, k+1, 1)
      sy = sj(i, j-1, k, 2) + sj(i+1, j-1, k, 2) + sj(i, j-1, k+1, 2) + &
&       sj(i+1, j-1, k+1, 2) + sj(i, j, k, 2) + sj(i+1, j, k, 2) + sj(i&
&       , j, k+1, 2) + sj(i+1, j, k+1, 2)
      sz = sj(i, j-1, k, 3) + sj(i+1, j-1, k, 3) + sj(i, j-1, k+1, 3) + &
&       sj(i+1, j-1, k+1, 3) + sj(i, j, k, 3) + sj(i+1, j, k, 3) + sj(i&
&       , j, k+1, 3) + sj(i+1, j, k+1, 3)
! compute the average velocities and speed of sound squared
! for this integration point. node that these variables are
! stored in w(ivx), w(ivy), w(ivz) and p.
      ubar = fourth*(w(i, j, k, ivx)+w(i+1, j, k, ivx)+w(i, j, k+1, ivx)&
&       +w(i+1, j, k+1, ivx))
      vbar = fourth*(w(i, j, k, ivy)+w(i+1, j, k, ivy)+w(i, j, k+1, ivy)&
&       +w(i+1, j, k+1, ivy))
      wbar = fourth*(w(i, j, k, ivz)+w(i+1, j, k, ivz)+w(i, j, k+1, ivz)&
&       +w(i+1, j, k+1, ivz))
      a2 = fourth*(aa(i, j, k)+aa(i+1, j, k)+aa(i, j, k+1)+aa(i+1, j, k+&
&       1))
! add the contributions to the surface integral to the node
! j-1 and substract it from the node j. for the heat flux it
! is reversed, because the negative of the gradient of the
! speed of sound must be computed.
      if (j .gt. 1) then
        call pushcontrol1b(0)
      else
        call pushcontrol1b(1)
      end if
      if (j .lt. je) then
        a2d = sz*qzd(i, j, k) + sy*qyd(i, j, k) + sx*qxd(i, j, k)
        szd = a2*qzd(i, j, k) - wbar*wzd(i, j, k) - vbar*vzd(i, j, k) - &
&         ubar*uzd(i, j, k)
        syd = a2*qyd(i, j, k) - wbar*wyd(i, j, k) - vbar*vyd(i, j, k) - &
&         ubar*uyd(i, j, k)
        sxd = a2*qxd(i, j, k) - wbar*wxd(i, j, k) - vbar*vxd(i, j, k) - &
&         ubar*uxd(i, j, k)
        wbard = -(sz*wzd(i, j, k)) - sy*wyd(i, j, k) - sx*wxd(i, j, k)
        vbard = -(sz*vzd(i, j, k)) - sy*vyd(i, j, k) - sx*vxd(i, j, k)
        ubard = -(sz*uzd(i, j, k)) - sy*uyd(i, j, k) - sx*uxd(i, j, k)
      else
        wbard = 0.0_8
        vbard = 0.0_8
        ubard = 0.0_8
        a2d = 0.0_8
        sxd = 0.0_8
        syd = 0.0_8
        szd = 0.0_8
      end if
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        a2d = a2d - sz*qzd(i, j-1, k) - sy*qyd(i, j-1, k) - sx*qxd(i, j-&
&         1, k)
        szd = szd + wbar*wzd(i, j-1, k) - a2*qzd(i, j-1, k) + vbar*vzd(i&
&         , j-1, k) + ubar*uzd(i, j-1, k)
        syd = syd + wbar*wyd(i, j-1, k) - a2*qyd(i, j-1, k) + vbar*vyd(i&
&         , j-1, k) + ubar*uyd(i, j-1, k)
        sxd = sxd + wbar*wxd(i, j-1, k) - a2*qxd(i, j-1, k) + vbar*vxd(i&
&         , j-1, k) + ubar*uxd(i, j-1, k)
        wbard = wbard + sz*wzd(i, j-1, k) + sy*wyd(i, j-1, k) + sx*wxd(i&
&         , j-1, k)
        vbard = vbard + sz*vzd(i, j-1, k) + sy*vyd(i, j-1, k) + sx*vxd(i&
&         , j-1, k)
        ubard = ubard + sz*uzd(i, j-1, k) + sy*uyd(i, j-1, k) + sx*uxd(i&
&         , j-1, k)
      end if
      tempd = fourth*a2d
      aad(i, j, k) = aad(i, j, k) + tempd
      aad(i+1, j, k) = aad(i+1, j, k) + tempd
      aad(i, j, k+1) = aad(i, j, k+1) + tempd
      aad(i+1, j, k+1) = aad(i+1, j, k+1) + tempd
      tempd = fourth*wbard
      wd(i, j, k, ivz) = wd(i, j, k, ivz) + tempd
      wd(i+1, j, k, ivz) = wd(i+1, j, k, ivz) + tempd
      wd(i, j, k+1, ivz) = wd(i, j, k+1, ivz) + tempd
      wd(i+1, j, k+1, ivz) = wd(i+1, j, k+1, ivz) + tempd
      tempd = fourth*vbard
      wd(i, j, k, ivy) = wd(i, j, k, ivy) + tempd
      wd(i+1, j, k, ivy) = wd(i+1, j, k, ivy) + tempd
      wd(i, j, k+1, ivy) = wd(i, j, k+1, ivy) + tempd
      wd(i+1, j, k+1, ivy) = wd(i+1, j, k+1, ivy) + tempd
      tempd = fourth*ubard
      wd(i, j, k, ivx) = wd(i, j, k, ivx) + tempd
      wd(i+1, j, k, ivx) = wd(i+1, j, k, ivx) + tempd
      wd(i, j, k+1, ivx) = wd(i, j, k+1, ivx) + tempd
      wd(i+1, j, k+1, ivx) = wd(i+1, j, k+1, ivx) + tempd
      sjd(i, j-1, k, 3) = sjd(i, j-1, k, 3) + szd
      sjd(i+1, j-1, k, 3) = sjd(i+1, j-1, k, 3) + szd
      sjd(i, j-1, k+1, 3) = sjd(i, j-1, k+1, 3) + szd
      sjd(i+1, j-1, k+1, 3) = sjd(i+1, j-1, k+1, 3) + szd
      sjd(i, j, k, 3) = sjd(i, j, k, 3) + szd
      sjd(i+1, j, k, 3) = sjd(i+1, j, k, 3) + szd
      sjd(i, j, k+1, 3) = sjd(i, j, k+1, 3) + szd
      sjd(i+1, j, k+1, 3) = sjd(i+1, j, k+1, 3) + szd
      sjd(i, j-1, k, 2) = sjd(i, j-1, k, 2) + syd
      sjd(i+1, j-1, k, 2) = sjd(i+1, j-1, k, 2) + syd
      sjd(i, j-1, k+1, 2) = sjd(i, j-1, k+1, 2) + syd
      sjd(i+1, j-1, k+1, 2) = sjd(i+1, j-1, k+1, 2) + syd
      sjd(i, j, k, 2) = sjd(i, j, k, 2) + syd
      sjd(i+1, j, k, 2) = sjd(i+1, j, k, 2) + syd
      sjd(i, j, k+1, 2) = sjd(i, j, k+1, 2) + syd
      sjd(i+1, j, k+1, 2) = sjd(i+1, j, k+1, 2) + syd
      sjd(i, j-1, k, 1) = sjd(i, j-1, k, 1) + sxd
      sjd(i+1, j-1, k, 1) = sjd(i+1, j-1, k, 1) + sxd
      sjd(i, j-1, k+1, 1) = sjd(i, j-1, k+1, 1) + sxd
      sjd(i+1, j-1, k+1, 1) = sjd(i+1, j-1, k+1, 1) + sxd
      sjd(i, j, k, 1) = sjd(i, j, k, 1) + sxd
      sjd(i+1, j, k, 1) = sjd(i+1, j, k, 1) + sxd
      sjd(i, j, k+1, 1) = sjd(i, j, k+1, 1) + sxd
      sjd(i+1, j, k+1, 1) = sjd(i+1, j, k+1, 1) + sxd
    end do
!$bwd-of ii-loop 
    do ii=0,il*jl*ke-1
      i = mod(ii, il) + 1
      j = mod(ii/il, jl) + 1
      k = ii/(il*jl) + 1
! compute 8 times the average normal for this part of
! the control volume. the factor 8 is taken care of later
! on when the division by the volume takes place.
      sx = sk(i, j, k-1, 1) + sk(i+1, j, k-1, 1) + sk(i, j+1, k-1, 1) + &
&       sk(i+1, j+1, k-1, 1) + sk(i, j, k, 1) + sk(i+1, j, k, 1) + sk(i&
&       , j+1, k, 1) + sk(i+1, j+1, k, 1)
      sy = sk(i, j, k-1, 2) + sk(i+1, j, k-1, 2) + sk(i, j+1, k-1, 2) + &
&       sk(i+1, j+1, k-1, 2) + sk(i, j, k, 2) + sk(i+1, j, k, 2) + sk(i&
&       , j+1, k, 2) + sk(i+1, j+1, k, 2)
      sz = sk(i, j, k-1, 3) + sk(i+1, j, k-1, 3) + sk(i, j+1, k-1, 3) + &
&       sk(i+1, j+1, k-1, 3) + sk(i, j, k, 3) + sk(i+1, j, k, 3) + sk(i&
&       , j+1, k, 3) + sk(i+1, j+1, k, 3)
! compute the average velocities and speed of sound squared
! for this integration point. node that these variables are
! stored in w(ivx), w(ivy), w(ivz) and p.
      ubar = fourth*(w(i, j, k, ivx)+w(i+1, j, k, ivx)+w(i, j+1, k, ivx)&
&       +w(i+1, j+1, k, ivx))
      vbar = fourth*(w(i, j, k, ivy)+w(i+1, j, k, ivy)+w(i, j+1, k, ivy)&
&       +w(i+1, j+1, k, ivy))
      wbar = fourth*(w(i, j, k, ivz)+w(i+1, j, k, ivz)+w(i, j+1, k, ivz)&
&       +w(i+1, j+1, k, ivz))
      a2 = fourth*(aa(i, j, k)+aa(i+1, j, k)+aa(i, j+1, k)+aa(i+1, j+1, &
&       k))
! add the contributions to the surface integral to the node
! j-1 and substract it from the node j. for the heat flux it
! is reversed, because the negative of the gradient of the
! speed of sound must be computed.
      if (k .gt. 1) then
        call pushcontrol1b(0)
      else
        call pushcontrol1b(1)
      end if
      if (k .lt. ke) then
        a2d = sz*qzd(i, j, k) + sy*qyd(i, j, k) + sx*qxd(i, j, k)
        szd = a2*qzd(i, j, k) - wbar*wzd(i, j, k) - vbar*vzd(i, j, k) - &
&         ubar*uzd(i, j, k)
        syd = a2*qyd(i, j, k) - wbar*wyd(i, j, k) - vbar*vyd(i, j, k) - &
&         ubar*uyd(i, j, k)
        sxd = a2*qxd(i, j, k) - wbar*wxd(i, j, k) - vbar*vxd(i, j, k) - &
&         ubar*uxd(i, j, k)
        wbard = -(sz*wzd(i, j, k)) - sy*wyd(i, j, k) - sx*wxd(i, j, k)
        vbard = -(sz*vzd(i, j, k)) - sy*vyd(i, j, k) - sx*vxd(i, j, k)
        ubard = -(sz*uzd(i, j, k)) - sy*uyd(i, j, k) - sx*uxd(i, j, k)
      else
        wbard = 0.0_8
        vbard = 0.0_8
        ubard = 0.0_8
        a2d = 0.0_8
        sxd = 0.0_8
        syd = 0.0_8
        szd = 0.0_8
      end if
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        a2d = a2d - sz*qzd(i, j, k-1) - sy*qyd(i, j, k-1) - sx*qxd(i, j&
&         , k-1)
        szd = szd + wbar*wzd(i, j, k-1) - a2*qzd(i, j, k-1) + vbar*vzd(i&
&         , j, k-1) + ubar*uzd(i, j, k-1)
        syd = syd + wbar*wyd(i, j, k-1) - a2*qyd(i, j, k-1) + vbar*vyd(i&
&         , j, k-1) + ubar*uyd(i, j, k-1)
        sxd = sxd + wbar*wxd(i, j, k-1) - a2*qxd(i, j, k-1) + vbar*vxd(i&
&         , j, k-1) + ubar*uxd(i, j, k-1)
        wbard = wbard + sz*wzd(i, j, k-1) + sy*wyd(i, j, k-1) + sx*wxd(i&
&         , j, k-1)
        vbard = vbard + sz*vzd(i, j, k-1) + sy*vyd(i, j, k-1) + sx*vxd(i&
&         , j, k-1)
        ubard = ubard + sz*uzd(i, j, k-1) + sy*uyd(i, j, k-1) + sx*uxd(i&
&         , j, k-1)
      end if
      tempd = fourth*a2d
      aad(i, j, k) = aad(i, j, k) + tempd
      aad(i+1, j, k) = aad(i+1, j, k) + tempd
      aad(i, j+1, k) = aad(i, j+1, k) + tempd
      aad(i+1, j+1, k) = aad(i+1, j+1, k) + tempd
      tempd = fourth*wbard
      wd(i, j, k, ivz) = wd(i, j, k, ivz) + tempd
      wd(i+1, j, k, ivz) = wd(i+1, j, k, ivz) + tempd
      wd(i, j+1, k, ivz) = wd(i, j+1, k, ivz) + tempd
      wd(i+1, j+1, k, ivz) = wd(i+1, j+1, k, ivz) + tempd
      tempd = fourth*vbard
      wd(i, j, k, ivy) = wd(i, j, k, ivy) + tempd
      wd(i+1, j, k, ivy) = wd(i+1, j, k, ivy) + tempd
      wd(i, j+1, k, ivy) = wd(i, j+1, k, ivy) + tempd
      wd(i+1, j+1, k, ivy) = wd(i+1, j+1, k, ivy) + tempd
      tempd = fourth*ubard
      wd(i, j, k, ivx) = wd(i, j, k, ivx) + tempd
      wd(i+1, j, k, ivx) = wd(i+1, j, k, ivx) + tempd
      wd(i, j+1, k, ivx) = wd(i, j+1, k, ivx) + tempd
      wd(i+1, j+1, k, ivx) = wd(i+1, j+1, k, ivx) + tempd
      skd(i, j, k-1, 3) = skd(i, j, k-1, 3) + szd
      skd(i+1, j, k-1, 3) = skd(i+1, j, k-1, 3) + szd
      skd(i, j+1, k-1, 3) = skd(i, j+1, k-1, 3) + szd
      skd(i+1, j+1, k-1, 3) = skd(i+1, j+1, k-1, 3) + szd
      skd(i, j, k, 3) = skd(i, j, k, 3) + szd
      skd(i+1, j, k, 3) = skd(i+1, j, k, 3) + szd
      skd(i, j+1, k, 3) = skd(i, j+1, k, 3) + szd
      skd(i+1, j+1, k, 3) = skd(i+1, j+1, k, 3) + szd
      skd(i, j, k-1, 2) = skd(i, j, k-1, 2) + syd
      skd(i+1, j, k-1, 2) = skd(i+1, j, k-1, 2) + syd
      skd(i, j+1, k-1, 2) = skd(i, j+1, k-1, 2) + syd
      skd(i+1, j+1, k-1, 2) = skd(i+1, j+1, k-1, 2) + syd
      skd(i, j, k, 2) = skd(i, j, k, 2) + syd
      skd(i+1, j, k, 2) = skd(i+1, j, k, 2) + syd
      skd(i, j+1, k, 2) = skd(i, j+1, k, 2) + syd
      skd(i+1, j+1, k, 2) = skd(i+1, j+1, k, 2) + syd
      skd(i, j, k-1, 1) = skd(i, j, k-1, 1) + sxd
      skd(i+1, j, k-1, 1) = skd(i+1, j, k-1, 1) + sxd
      skd(i, j+1, k-1, 1) = skd(i, j+1, k-1, 1) + sxd
      skd(i+1, j+1, k-1, 1) = skd(i+1, j+1, k-1, 1) + sxd
      skd(i, j, k, 1) = skd(i, j, k, 1) + sxd
      skd(i+1, j, k, 1) = skd(i+1, j, k, 1) + sxd
      skd(i, j+1, k, 1) = skd(i, j+1, k, 1) + sxd
      skd(i+1, j+1, k, 1) = skd(i+1, j+1, k, 1) + sxd
    end do
    qzd = 0.0_8
    qyd = 0.0_8
    qxd = 0.0_8
    wzd = 0.0_8
    wyd = 0.0_8
    wxd = 0.0_8
    vzd = 0.0_8
    vyd = 0.0_8
    vxd = 0.0_8
    uzd = 0.0_8
    uyd = 0.0_8
    uxd = 0.0_8
  end subroutine allnodalgradients_b
 
  subroutine allnodalgradients()
!
!         nodalgradients computes the nodal velocity gradients and
!         minus the gradient of the speed of sound squared. the minus
!         sign is present, because this is the definition of the heat
!         flux. these gradients are computed for all nodes.
!
    use constants
    use blockpointers
    implicit none
!        local variables.
    integer(kind=inttype) :: i, j, k
    integer(kind=inttype) :: k1, kk
    integer(kind=inttype) :: istart, iend, isize, ii
    integer(kind=inttype) :: jstart, jend, jsize
    integer(kind=inttype) :: kstart, kend, ksize
    real(kind=realtype) :: oneoverv, ubar, vbar, wbar, a2
    real(kind=realtype) :: sx, sx1, sy, sy1, sz, sz1
    intrinsic mod
! zero all nodeal gradients:
    ux = zero
    uy = zero
    uz = zero
    vx = zero
    vy = zero
    vz = zero
    wx = zero
    wy = zero
    wz = zero
    qx = zero
    qy = zero
    qz = zero
!$ad ii-loop
! first part. contribution in the k-direction.
! the contribution is scattered to both the left and right node
! in k-direction.
    do ii=0,il*jl*ke-1
      i = mod(ii, il) + 1
      j = mod(ii/il, jl) + 1
      k = ii/(il*jl) + 1
! compute 8 times the average normal for this part of
! the control volume. the factor 8 is taken care of later
! on when the division by the volume takes place.
      sx = sk(i, j, k-1, 1) + sk(i+1, j, k-1, 1) + sk(i, j+1, k-1, 1) + &
&       sk(i+1, j+1, k-1, 1) + sk(i, j, k, 1) + sk(i+1, j, k, 1) + sk(i&
&       , j+1, k, 1) + sk(i+1, j+1, k, 1)
      sy = sk(i, j, k-1, 2) + sk(i+1, j, k-1, 2) + sk(i, j+1, k-1, 2) + &
&       sk(i+1, j+1, k-1, 2) + sk(i, j, k, 2) + sk(i+1, j, k, 2) + sk(i&
&       , j+1, k, 2) + sk(i+1, j+1, k, 2)
      sz = sk(i, j, k-1, 3) + sk(i+1, j, k-1, 3) + sk(i, j+1, k-1, 3) + &
&       sk(i+1, j+1, k-1, 3) + sk(i, j, k, 3) + sk(i+1, j, k, 3) + sk(i&
&       , j+1, k, 3) + sk(i+1, j+1, k, 3)
! compute the average velocities and speed of sound squared
! for this integration point. node that these variables are
! stored in w(ivx), w(ivy), w(ivz) and p.
      ubar = fourth*(w(i, j, k, ivx)+w(i+1, j, k, ivx)+w(i, j+1, k, ivx)&
&       +w(i+1, j+1, k, ivx))
      vbar = fourth*(w(i, j, k, ivy)+w(i+1, j, k, ivy)+w(i, j+1, k, ivy)&
&       +w(i+1, j+1, k, ivy))
      wbar = fourth*(w(i, j, k, ivz)+w(i+1, j, k, ivz)+w(i, j+1, k, ivz)&
&       +w(i+1, j+1, k, ivz))
      a2 = fourth*(aa(i, j, k)+aa(i+1, j, k)+aa(i, j+1, k)+aa(i+1, j+1, &
&       k))
! add the contributions to the surface integral to the node
! j-1 and substract it from the node j. for the heat flux it
! is reversed, because the negative of the gradient of the
! speed of sound must be computed.
      if (k .gt. 1) then
        ux(i, j, k-1) = ux(i, j, k-1) + ubar*sx
        uy(i, j, k-1) = uy(i, j, k-1) + ubar*sy
        uz(i, j, k-1) = uz(i, j, k-1) + ubar*sz
        vx(i, j, k-1) = vx(i, j, k-1) + vbar*sx
        vy(i, j, k-1) = vy(i, j, k-1) + vbar*sy
        vz(i, j, k-1) = vz(i, j, k-1) + vbar*sz
        wx(i, j, k-1) = wx(i, j, k-1) + wbar*sx
        wy(i, j, k-1) = wy(i, j, k-1) + wbar*sy
        wz(i, j, k-1) = wz(i, j, k-1) + wbar*sz
        qx(i, j, k-1) = qx(i, j, k-1) - a2*sx
        qy(i, j, k-1) = qy(i, j, k-1) - a2*sy
        qz(i, j, k-1) = qz(i, j, k-1) - a2*sz
      end if
      if (k .lt. ke) then
        ux(i, j, k) = ux(i, j, k) - ubar*sx
        uy(i, j, k) = uy(i, j, k) - ubar*sy
        uz(i, j, k) = uz(i, j, k) - ubar*sz
        vx(i, j, k) = vx(i, j, k) - vbar*sx
        vy(i, j, k) = vy(i, j, k) - vbar*sy
        vz(i, j, k) = vz(i, j, k) - vbar*sz
        wx(i, j, k) = wx(i, j, k) - wbar*sx
        wy(i, j, k) = wy(i, j, k) - wbar*sy
        wz(i, j, k) = wz(i, j, k) - wbar*sz
        qx(i, j, k) = qx(i, j, k) + a2*sx
        qy(i, j, k) = qy(i, j, k) + a2*sy
        qz(i, j, k) = qz(i, j, k) + a2*sz
      end if
    end do
!$ad ii-loop
! second part. contribution in the j-direction.
! the contribution is scattered to both the left and right node
! in j-direction.
    do ii=0,il*je*kl-1
      i = mod(ii, il) + 1
      j = mod(ii/il, je) + 1
      k = ii/(il*je) + 1
! compute 8 times the average normal for this part of
! the control volume. the factor 8 is taken care of later
! on when the division by the volume takes place.
      sx = sj(i, j-1, k, 1) + sj(i+1, j-1, k, 1) + sj(i, j-1, k+1, 1) + &
&       sj(i+1, j-1, k+1, 1) + sj(i, j, k, 1) + sj(i+1, j, k, 1) + sj(i&
&       , j, k+1, 1) + sj(i+1, j, k+1, 1)
      sy = sj(i, j-1, k, 2) + sj(i+1, j-1, k, 2) + sj(i, j-1, k+1, 2) + &
&       sj(i+1, j-1, k+1, 2) + sj(i, j, k, 2) + sj(i+1, j, k, 2) + sj(i&
&       , j, k+1, 2) + sj(i+1, j, k+1, 2)
      sz = sj(i, j-1, k, 3) + sj(i+1, j-1, k, 3) + sj(i, j-1, k+1, 3) + &
&       sj(i+1, j-1, k+1, 3) + sj(i, j, k, 3) + sj(i+1, j, k, 3) + sj(i&
&       , j, k+1, 3) + sj(i+1, j, k+1, 3)
! compute the average velocities and speed of sound squared
! for this integration point. node that these variables are
! stored in w(ivx), w(ivy), w(ivz) and p.
      ubar = fourth*(w(i, j, k, ivx)+w(i+1, j, k, ivx)+w(i, j, k+1, ivx)&
&       +w(i+1, j, k+1, ivx))
      vbar = fourth*(w(i, j, k, ivy)+w(i+1, j, k, ivy)+w(i, j, k+1, ivy)&
&       +w(i+1, j, k+1, ivy))
      wbar = fourth*(w(i, j, k, ivz)+w(i+1, j, k, ivz)+w(i, j, k+1, ivz)&
&       +w(i+1, j, k+1, ivz))
      a2 = fourth*(aa(i, j, k)+aa(i+1, j, k)+aa(i, j, k+1)+aa(i+1, j, k+&
&       1))
! add the contributions to the surface integral to the node
! j-1 and substract it from the node j. for the heat flux it
! is reversed, because the negative of the gradient of the
! speed of sound must be computed.
      if (j .gt. 1) then
        ux(i, j-1, k) = ux(i, j-1, k) + ubar*sx
        uy(i, j-1, k) = uy(i, j-1, k) + ubar*sy
        uz(i, j-1, k) = uz(i, j-1, k) + ubar*sz
        vx(i, j-1, k) = vx(i, j-1, k) + vbar*sx
        vy(i, j-1, k) = vy(i, j-1, k) + vbar*sy
        vz(i, j-1, k) = vz(i, j-1, k) + vbar*sz
        wx(i, j-1, k) = wx(i, j-1, k) + wbar*sx
        wy(i, j-1, k) = wy(i, j-1, k) + wbar*sy
        wz(i, j-1, k) = wz(i, j-1, k) + wbar*sz
        qx(i, j-1, k) = qx(i, j-1, k) - a2*sx
        qy(i, j-1, k) = qy(i, j-1, k) - a2*sy
        qz(i, j-1, k) = qz(i, j-1, k) - a2*sz
      end if
      if (j .lt. je) then
        ux(i, j, k) = ux(i, j, k) - ubar*sx
        uy(i, j, k) = uy(i, j, k) - ubar*sy
        uz(i, j, k) = uz(i, j, k) - ubar*sz
        vx(i, j, k) = vx(i, j, k) - vbar*sx
        vy(i, j, k) = vy(i, j, k) - vbar*sy
        vz(i, j, k) = vz(i, j, k) - vbar*sz
        wx(i, j, k) = wx(i, j, k) - wbar*sx
        wy(i, j, k) = wy(i, j, k) - wbar*sy
        wz(i, j, k) = wz(i, j, k) - wbar*sz
        qx(i, j, k) = qx(i, j, k) + a2*sx
        qy(i, j, k) = qy(i, j, k) + a2*sy
        qz(i, j, k) = qz(i, j, k) + a2*sz
      end if
    end do
!$ad ii-loop
! third part. contribution in the i-direction.
! the contribution is scattered to both the left and right node
! in i-direction.
    do ii=0,ie*jl*kl-1
      i = mod(ii, ie) + 1
      j = mod(ii/ie, jl) + 1
      k = ii/(ie*jl) + 1
! compute 8 times the average normal for this part of
! the control volume. the factor 8 is taken care of later
! on when the division by the volume takes place.
      sx = si(i-1, j, k, 1) + si(i-1, j+1, k, 1) + si(i-1, j, k+1, 1) + &
&       si(i-1, j+1, k+1, 1) + si(i, j, k, 1) + si(i, j+1, k, 1) + si(i&
&       , j, k+1, 1) + si(i, j+1, k+1, 1)
      sy = si(i-1, j, k, 2) + si(i-1, j+1, k, 2) + si(i-1, j, k+1, 2) + &
&       si(i-1, j+1, k+1, 2) + si(i, j, k, 2) + si(i, j+1, k, 2) + si(i&
&       , j, k+1, 2) + si(i, j+1, k+1, 2)
      sz = si(i-1, j, k, 3) + si(i-1, j+1, k, 3) + si(i-1, j, k+1, 3) + &
&       si(i-1, j+1, k+1, 3) + si(i, j, k, 3) + si(i, j+1, k, 3) + si(i&
&       , j, k+1, 3) + si(i, j+1, k+1, 3)
! compute the average velocities and speed of sound squared
! for this integration point. node that these variables are
! stored in w(ivx), w(ivy), w(ivz) and p.
      ubar = fourth*(w(i, j, k, ivx)+w(i, j+1, k, ivx)+w(i, j, k+1, ivx)&
&       +w(i, j+1, k+1, ivx))
      vbar = fourth*(w(i, j, k, ivy)+w(i, j+1, k, ivy)+w(i, j, k+1, ivy)&
&       +w(i, j+1, k+1, ivy))
      wbar = fourth*(w(i, j, k, ivz)+w(i, j+1, k, ivz)+w(i, j, k+1, ivz)&
&       +w(i, j+1, k+1, ivz))
      a2 = fourth*(aa(i, j, k)+aa(i, j+1, k)+aa(i, j, k+1)+aa(i, j+1, k+&
&       1))
! add the contributions to the surface integral to the node
! j-1 and substract it from the node j. for the heat flux it
! is reversed, because the negative of the gradient of the
! speed of sound must be computed.
      if (i .gt. 1) then
        ux(i-1, j, k) = ux(i-1, j, k) + ubar*sx
        uy(i-1, j, k) = uy(i-1, j, k) + ubar*sy
        uz(i-1, j, k) = uz(i-1, j, k) + ubar*sz
        vx(i-1, j, k) = vx(i-1, j, k) + vbar*sx
        vy(i-1, j, k) = vy(i-1, j, k) + vbar*sy
        vz(i-1, j, k) = vz(i-1, j, k) + vbar*sz
        wx(i-1, j, k) = wx(i-1, j, k) + wbar*sx
        wy(i-1, j, k) = wy(i-1, j, k) + wbar*sy
        wz(i-1, j, k) = wz(i-1, j, k) + wbar*sz
        qx(i-1, j, k) = qx(i-1, j, k) - a2*sx
        qy(i-1, j, k) = qy(i-1, j, k) - a2*sy
        qz(i-1, j, k) = qz(i-1, j, k) - a2*sz
      end if
      if (i .lt. ie) then
        ux(i, j, k) = ux(i, j, k) - ubar*sx
        uy(i, j, k) = uy(i, j, k) - ubar*sy
        uz(i, j, k) = uz(i, j, k) - ubar*sz
        vx(i, j, k) = vx(i, j, k) - vbar*sx
        vy(i, j, k) = vy(i, j, k) - vbar*sy
        vz(i, j, k) = vz(i, j, k) - vbar*sz
        wx(i, j, k) = wx(i, j, k) - wbar*sx
        wy(i, j, k) = wy(i, j, k) - wbar*sy
        wz(i, j, k) = wz(i, j, k) - wbar*sz
        qx(i, j, k) = qx(i, j, k) + a2*sx
        qy(i, j, k) = qy(i, j, k) + a2*sy
        qz(i, j, k) = qz(i, j, k) + a2*sz
      end if
    end do
!$ad ii-loop
! divide by 8 times the volume to obtain the correct gradients.
    do ii=0,il*jl*kl-1
      i = mod(ii, il) + 1
      j = mod(ii/il, jl) + 1
      k = ii/(il*jl) + 1
! compute the inverse of 8 times the volume for this node.
      oneoverv = one/(vol(i, j, k)+vol(i, j, k+1)+vol(i+1, j, k)+vol(i+1&
&       , j, k+1)+vol(i, j+1, k)+vol(i, j+1, k+1)+vol(i+1, j+1, k)+vol(i&
&       +1, j+1, k+1))
! compute the correct velocity gradients and "unit" heat
! fluxes. the velocity gradients are stored in ux, etc.
      ux(i, j, k) = ux(i, j, k)*oneoverv
      uy(i, j, k) = uy(i, j, k)*oneoverv
      uz(i, j, k) = uz(i, j, k)*oneoverv
      vx(i, j, k) = vx(i, j, k)*oneoverv
      vy(i, j, k) = vy(i, j, k)*oneoverv
      vz(i, j, k) = vz(i, j, k)*oneoverv
      wx(i, j, k) = wx(i, j, k)*oneoverv
      wy(i, j, k) = wy(i, j, k)*oneoverv
      wz(i, j, k) = wz(i, j, k)*oneoverv
      qx(i, j, k) = qx(i, j, k)*oneoverv
      qy(i, j, k) = qy(i, j, k)*oneoverv
      qz(i, j, k) = qz(i, j, k)*oneoverv
    end do
  end subroutine allnodalgradients
! ----------------------------------------------------------------------
!                                                                      |
!                    no tapenade routine below this line               |
!                                                                      |
! ----------------------------------------------------------------------
 
end module flowutils_b