flowUtils_d.f90

Go to the documentation of this file.

!        generated by tapenade     (inria, ecuador team)
!  tapenade 3.16 (develop) - 22 aug 2023 15:51
!
module flowutils_d
  implicit none
 
contains
!  differentiation of computettot in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: ttot
!   with respect to varying inputs: rgas p u v w ttot rho
  subroutine computettot_d(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_d, only : terminate
    implicit none
!
!      subroutine arguments.
!
    real(kind=realtype), intent(in) :: rho, p, u, v, w
    real(kind=realtype), intent(in) :: rhod, pd, ud, vd, wd
    real(kind=realtype), intent(out) :: ttot
    real(kind=realtype), intent(out) :: ttotd
!
!      local variables.
!
    integer(kind=inttype) :: i
    real(kind=realtype) :: govgm1, t, kin
    real(kind=realtype) :: td, kind0
    real(kind=realtype) :: temp
! 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)
      temp = p/(rho*rgas)
      td = (pd-temp*(rgas*rhod+rho*rgasd))/(rho*rgas)
      t = temp
      kind0 = half*(2*u*ud+2*v*vd+2*w*wd)
      kin = half*(u*u+v*v+w*w)
      temp = rho*kin/(govgm1*p)
      ttotd = (one+temp)*td + t*(kin*rhod+rho*kind0-temp*govgm1*pd)/(&
&       govgm1*p)
      ttot = t*(one+temp)
!===============================================================
    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_d
 
  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_d, 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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: ptot
!   with respect to varying inputs: p u v w ptot rho
  subroutine computeptot_d(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), intent(in) :: rhod, pd, ud, vd, wd
    real(kind=realtype), intent(out) :: ptot
    real(kind=realtype), intent(out) :: 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) :: temp0
    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)
      kind0 = half*(2*u*ud+2*v*vd+2*w*wd)
      kin = half*(u*u+v*v+w*w)
      temp = rho*kin/(govgm1*p)
      temp0 = (one+temp)**govgm1
      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)*(kin*rhod+rho*kind0-temp*govgm1*&
&         pd)/p
      end if
      ptotd = temp0*pd + p*tempd
      ptot = p*temp0
!===============================================================
    end select
  end subroutine computeptot_d
 
  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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *aa
!   with respect to varying inputs: *p *w
!   rw status of diff variables: *aa:out *p:in *w:in
!   plus diff mem management of: aa:in p:in w:in
  subroutine computespeedofsoundsquared_d()
!
!       computespeedofsoundsquared does what it says.
!
    use constants
    use blockpointers, only : ie, je, ke, w, wd, p, pd, aa, aad, gamma
    use utils_d, 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
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
! determine if we need to correct for k
    correctfork = getcorrectfork()
    if (correctfork) then
      if (associated(aad)) aad = 0.0_8
      do k=1,ke
        do j=1,je
          do i=1,ie
            temp = w(i, j, k, itu1)
            temp0 = w(i, j, k, irho)
            ppd = pd(i, j, k) - twothird*(temp*wd(i, j, k, irho)+temp0*&
&             wd(i, j, k, itu1))
            pp = p(i, j, k) - twothird*(temp0*temp)
            temp0 = w(i, j, k, irho)
            aad(i, j, k) = gamma(i, j, k)*(ppd-pp*wd(i, j, k, irho)/&
&             temp0)/temp0
            aa(i, j, k) = gamma(i, j, k)*(pp/temp0)
          end do
        end do
      end do
    else
      if (associated(aad)) aad = 0.0_8
      do k=1,ke
        do j=1,je
          do i=1,ie
            temp0 = w(i, j, k, irho)
            temp = p(i, j, k)/temp0
            aad(i, j, k) = gamma(i, j, k)*(pd(i, j, k)-temp*wd(i, j, k, &
&             irho))/temp0
            aa(i, j, k) = gamma(i, j, k)*temp
          end do
        end do
      end do
    end if
  end subroutine computespeedofsoundsquared_d
 
  subroutine computespeedofsoundsquared()
!
!       computespeedofsoundsquared does what it says.
!
    use constants
    use blockpointers, only : ie, je, ke, w, p, aa, gamma
    use utils_d, 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
! determine if we need to correct for k
    correctfork = getcorrectfork()
    if (correctfork) then
      do k=1,ke
        do j=1,je
          do i=1,ie
            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
        end do
      end do
    else
      do k=1,ke
        do j=1,je
          do i=1,ie
            aa(i, j, k) = gamma(i, j, k)*p(i, j, k)/w(i, j, k, irho)
          end do
        end do
      end do
    end if
  end subroutine computespeedofsoundsquared
 
!  differentiation of computeetotblock in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *w
!   with respect to varying inputs: *p *w
!   rw status of diff variables: *p:in *w:in-out
!   plus diff mem management of: p:in w:in
  subroutine computeetotblock_d(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
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype) :: temp1
    real(kind=realtype) :: temp2
    real(kind=realtype) :: temp3
    real(kind=realtype) :: temp4
    real(kind=realtype) :: temp5
!
! 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
        do k=kstart,kend
          do j=jstart,jend
            do i=istart,iend
              temp = w(i, j, k, ivz)
              temp0 = w(i, j, k, ivy)
              temp1 = w(i, j, k, ivx)
              temp2 = temp1*temp1 + temp0*temp0 + temp*temp
              temp3 = w(i, j, k, irho)
              temp4 = w(i, j, k, itu1)
              temp5 = w(i, j, k, irho)
              wd(i, j, k, irhoe) = ovgm1*pd(i, j, k) + half*(temp2*wd(i&
&               , j, k, irho)+temp3*(2*temp1*wd(i, j, k, ivx)+2*temp0*wd&
&               (i, j, k, ivy)+2*temp*wd(i, j, k, ivz))) - factk*(temp4*&
&               wd(i, j, k, irho)+temp5*wd(i, j, k, itu1))
              w(i, j, k, irhoe) = ovgm1*p(i, j, k) + half*(temp3*temp2) &
&               - factk*(temp5*temp4)
            end do
          end do
        end do
      else
        do k=kstart,kend
          do j=jstart,jend
            do i=istart,iend
              temp5 = w(i, j, k, ivz)
              temp4 = w(i, j, k, ivy)
              temp3 = w(i, j, k, ivx)
              temp2 = temp3*temp3 + temp4*temp4 + temp5*temp5
              temp1 = w(i, j, k, irho)
              wd(i, j, k, irhoe) = ovgm1*pd(i, j, k) + half*(temp2*wd(i&
&               , j, k, irho)+temp1*(2*temp3*wd(i, j, k, ivx)+2*temp4*wd&
&               (i, j, k, ivy)+2*temp5*wd(i, j, k, ivz)))
              w(i, j, k, irhoe) = ovgm1*p(i, j, k) + half*(temp1*temp2)
            end do
          end do
        end do
      end if
    end select
  end subroutine computeetotblock_d
 
  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
!
! 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
        do k=kstart,kend
          do j=jstart,jend
            do i=istart,iend
              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
          end do
        end do
      else
        do k=kstart,kend
          do j=jstart,jend
            do i=istart,iend
              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 do
        end do
      end if
    end select
  end subroutine computeetotblock
 
!  differentiation of etot in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: etotal
!   with respect to varying inputs: k p u v w etotal rho
  subroutine etot_d(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), intent(in) :: rhod, pd, kd
    real(kind=realtype), intent(in) :: u, v, w
    real(kind=realtype), intent(in) :: ud, vd, wd
    real(kind=realtype), intent(out) :: etotal
    real(kind=realtype), intent(out) :: etotald
    logical, intent(in) :: correctfork
!
!      local variables.
!
    integer(kind=inttype) :: i
    real(kind=realtype) :: temp
! compute the internal energy for unit mass.
    call eint_d(rho, rhod, p, pd, k, kd, etotal, etotald, correctfork)
    temp = etotal + half*(u*u+v*v+w*w)
    etotald = temp*rhod + rho*(etotald+half*(2*u*ud+2*v*vd+2*w*wd))
    etotal = rho*temp
  end subroutine etot_d
 
  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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: einternal
!   with respect to varying inputs: k p einternal rho
!      ==================================================================
  subroutine eint_d(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), intent(in) :: rhod, pd, kd
    real(kind=realtype), intent(out) :: einternal
    real(kind=realtype), intent(out) :: 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
! 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.
      einternald = ovgm1*(pd-p*rhod/rho)/rho
      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)
        einternald = einternald - factk*kd
        einternal = einternal - factk*k
      end if
    end select
  end subroutine eint_d
 
!      ==================================================================
  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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *p
!   with respect to varying inputs: pinfcorr *w
!   rw status of diff variables: pinfcorr:in *p:out *w:in
!   plus diff mem management of: p:in w:in
  subroutine computepressuresimple_d(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 max
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype) :: temp1
! compute the pressures
    gm1 = gammaconstant - one
    if (includehalos) then
      ibeg = 0
      jbeg = 0
      kbeg = 0
      iend = ib
      jend = jb
      kend = kb
      if (associated(pd)) pd = 0.0_8
    else
      ibeg = 2
      jbeg = 2
      kbeg = 2
      iend = il
      jend = jl
      kend = kl
      if (associated(pd)) pd = 0.0_8
    end if
    do k=kbeg,kend
      do j=jbeg,jend
        do i=ibeg,iend
          temp = w(i, j, k, ivx)
          temp0 = w(i, j, k, ivy)
          temp1 = w(i, j, k, ivz)
          v2d = 2*temp*wd(i, j, k, ivx) + 2*temp0*wd(i, j, k, ivy) + 2*&
&           temp1*wd(i, j, k, ivz)
          v2 = temp*temp + temp0*temp0 + temp1*temp1
          temp1 = w(i, j, k, irho)
          pd(i, j, k) = gm1*(wd(i, j, k, irhoe)-half*(v2*wd(i, j, k, &
&           irho)+temp1*v2d))
          p(i, j, k) = gm1*(w(i, j, k, irhoe)-half*(temp1*v2))
          if (p(i, j, k) .lt. 1.e-4_realtype*pinfcorr) then
            pd(i, j, k) = 1.e-4_realtype*pinfcorrd
            p(i, j, k) = 1.e-4_realtype*pinfcorr
          else
            p(i, j, k) = p(i, j, k)
          end if
        end do
      end do
    end do
  end subroutine computepressuresimple_d
 
  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 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
    do k=kbeg,kend
      do j=jbeg,jend
        do i=ibeg,iend
          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 do
    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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *rlv
!   with respect to varying inputs: muref tref rgas *p *w
!   rw status of diff variables: muref:in tref:in rgas:in *p:in
!                *w:in *rlv:out
!   plus diff mem management of: p:in w:in rlv:in
  subroutine computelamviscosity_d(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_d, 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
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype) :: temp1
! return immediately if no laminar viscosity needs to be computed.
    if (.not.viscous) then
      if (associated(rlvd)) rlvd = 0.0_8
      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.
      musuthd = -(musuthdim*murefd/muref**2)
      musuth = musuthdim/muref
      tsuthd = -(tsuthdim*trefd/tref**2)
      tsuth = tsuthdim/tref
      ssuthd = -(ssuthdim*trefd/tref**2)
      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
        if (associated(rlvd)) rlvd = 0.0_8
        do k=kbeg,kend
          do j=jbeg,jend
            do i=ibeg,iend
              temp = w(i, j, k, itu1)
              temp0 = w(i, j, k, irho)
              ppd = pd(i, j, k) - twothird*(temp*wd(i, j, k, irho)+temp0&
&               *wd(i, j, k, itu1))
              pp = p(i, j, k) - twothird*(temp0*temp)
              temp0 = w(i, j, k, irho)
              temp = rgas*temp0
              td = (ppd-pp*(temp0*rgasd+rgas*wd(i, j, k, irho))/temp)/&
&               temp
              t = pp/temp
              temp0 = t/tsuth
              temp = temp0**1.5_realtype
              temp1 = musuth*(tsuth+ssuth)/(t+ssuth)
              rlvd(i, j, k) = temp*((tsuth+ssuth)*musuthd+musuth*(tsuthd&
&               +ssuthd)-temp1*(td+ssuthd))/(t+ssuth) + temp1*&
&               1.5_realtype*temp0**0.5*(td-temp0*tsuthd)/tsuth
              rlv(i, j, k) = temp1*temp
            end do
          end do
        end do
      else
        if (associated(rlvd)) rlvd = 0.0_8
! loop over the owned cells *and* first level halos of this
! block and compute the laminar viscosity ratio.
        do k=kbeg,kend
          do j=jbeg,jend
            do i=ibeg,iend
! compute the nondimensional temperature and the
! nondimensional laminar viscosity.
              temp1 = w(i, j, k, irho)
              temp0 = p(i, j, k)/(rgas*temp1)
              td = (pd(i, j, k)-temp0*(temp1*rgasd+rgas*wd(i, j, k, irho&
&               )))/(rgas*temp1)
              t = temp0
              temp1 = t/tsuth
              temp0 = temp1**1.5_realtype
              temp = musuth*(tsuth+ssuth)/(t+ssuth)
              rlvd(i, j, k) = temp0*((tsuth+ssuth)*musuthd+musuth*(&
&               tsuthd+ssuthd)-temp*(td+ssuthd))/(t+ssuth) + temp*&
&               1.5_realtype*temp1**0.5*(td-temp1*tsuthd)/tsuth
              rlv(i, j, k) = temp*temp0
            end do
          end do
        end do
      end if
    end if
  end subroutine computelamviscosity_d
 
  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_d, 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
! 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
        do k=kbeg,kend
          do j=jbeg,jend
            do i=ibeg,iend
              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
          end do
        end do
      else
! loop over the owned cells *and* first level halos of this
! block and compute the laminar viscosity ratio.
        do k=kbeg,kend
          do j=jbeg,jend
            do i=ibeg,iend
! 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 do
        end do
      end if
    end if
  end subroutine computelamviscosity
 
!  differentiation of adjustinflowangle in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: veldirfreestream dragdirection
!                liftdirection
!   with respect to varying inputs: alpha beta
!   rw status of diff variables: alpha:in veldirfreestream:out
!                beta:in dragdirection:out liftdirection:out
  subroutine adjustinflowangle_d()
    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
    call getdirvector_d(refdir1, alpha, alphad, beta, betad, &
&                 veldirfreestream, veldirfreestreamd, 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_d(refdir1, alpha, alphad, beta, betad, &
&                 dragdirection, dragdirectiond, 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_d(refdir2, alpha, alphad, beta, betad, &
&                 liftdirection, liftdirectiond, liftindex)
  end subroutine adjustinflowangle_d
 
  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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: rotationmatrix
!   with respect to varying inputs: timeref
  subroutine derivativerotmatrixrigid_d(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_d, only : rigidrotangle, derivativerigidrotangle, &
&   derivativerigidrotangle_d
    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
    real(kind=realtype), dimension(3, 3), intent(out) :: 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
    real(kind=realtype) :: temp
! 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.
    dphixd = derivativerigidrotangle_d(degreepolxrot, coefpolxrot, &
&     degreefourxrot, omegafourxrot, coscoeffourxrot, sincoeffourxrot, t&
&     , dphix)
    dphiyd = derivativerigidrotangle_d(degreepolyrot, coefpolyrot, &
&     degreefouryrot, omegafouryrot, coscoeffouryrot, sincoeffouryrot, t&
&     , dphiy)
    dphizd = derivativerigidrotangle_d(degreepolzrot, coefpolzrot, &
&     degreefourzrot, omegafourzrot, coscoeffourzrot, sincoeffourzrot, t&
&     , dphiz)
! 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
    dmd = 0.0_8
    dmd(1, 1) = -(cosy*sinz*dphizd)
    dm(1, 1) = -(cosy*sinz*dphiz)
    temp = -(cosx*cosz) - sinx*siny*sinz
    dmd(1, 2) = temp*dphizd
    dm(1, 2) = temp*dphiz
    temp = sinx*cosz - cosx*siny*sinz
    dmd(1, 3) = temp*dphizd
    dm(1, 3) = temp*dphiz
    dmd(2, 1) = cosy*cosz*dphizd
    dm(2, 1) = cosy*cosz*dphiz
    temp = sinx*siny*cosz - cosx*sinz
    dmd(2, 2) = temp*dphizd
    dm(2, 2) = temp*dphiz
    temp = cosx*siny*cosz + sinx*sinz
    dmd(2, 3) = temp*dphizd
    dm(2, 3) = temp*dphiz
    dmd(3, 1) = 0.0_8
    dm(3, 1) = zero
    dmd(3, 2) = 0.0_8
    dm(3, 2) = zero
    dmd(3, 3) = 0.0_8
    dm(3, 3) = zero
! part 2: mz * dmy * mx.
    dmd(1, 1) = dmd(1, 1) - siny*cosz*dphiyd
    dm(1, 1) = dm(1, 1) - siny*cosz*dphiy
    dmd(1, 2) = dmd(1, 2) + sinx*cosy*cosz*dphiyd
    dm(1, 2) = dm(1, 2) + sinx*cosy*cosz*dphiy
    dmd(1, 3) = dmd(1, 3) + cosx*cosy*cosz*dphiyd
    dm(1, 3) = dm(1, 3) + cosx*cosy*cosz*dphiy
    dmd(2, 1) = dmd(2, 1) - siny*sinz*dphiyd
    dm(2, 1) = dm(2, 1) - siny*sinz*dphiy
    dmd(2, 2) = dmd(2, 2) + sinx*cosy*sinz*dphiyd
    dm(2, 2) = dm(2, 2) + sinx*cosy*sinz*dphiy
    dmd(2, 3) = dmd(2, 3) + cosx*cosy*sinz*dphiyd
    dm(2, 3) = dm(2, 3) + cosx*cosy*sinz*dphiy
    dmd(3, 1) = dmd(3, 1) - cosy*dphiyd
    dm(3, 1) = dm(3, 1) - cosy*dphiy
    dmd(3, 2) = dmd(3, 2) - sinx*siny*dphiyd
    dm(3, 2) = dm(3, 2) - sinx*siny*dphiy
    dmd(3, 3) = dmd(3, 3) - cosx*siny*dphiyd
    dm(3, 3) = dm(3, 3) - cosx*siny*dphiy
! part 3: mz * my * dmx
    temp = sinx*sinz + cosx*siny*cosz
    dmd(1, 2) = dmd(1, 2) + temp*dphixd
    dm(1, 2) = dm(1, 2) + temp*dphix
    temp = cosx*sinz - sinx*siny*cosz
    dmd(1, 3) = dmd(1, 3) + temp*dphixd
    dm(1, 3) = dm(1, 3) + temp*dphix
    temp = cosx*siny*sinz - sinx*cosz
    dmd(2, 2) = dmd(2, 2) + temp*dphixd
    dm(2, 2) = dm(2, 2) + temp*dphix
    temp = sinx*siny*sinz + cosx*cosz
    dmd(2, 3) = dmd(2, 3) - temp*dphixd
    dm(2, 3) = dm(2, 3) - temp*dphix
    dmd(3, 2) = dmd(3, 2) + cosx*cosy*dphixd
    dm(3, 2) = dm(3, 2) + cosx*cosy*dphix
    dmd(3, 3) = dmd(3, 3) - sinx*cosy*dphixd
    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
    rotationmatrixd = 0.0_8
! 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
        rotationmatrixd(i, j) = m(j, 1)*dmd(i, 1) + m(j, 2)*dmd(i, 2) + &
&         m(j, 3)*dmd(i, 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_d
 
  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_d, 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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: winddirection
!   with respect to varying inputs: alpha beta
  subroutine getdirvector_d(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_d, 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), intent(out) :: winddirection
    real(kind=realtype), dimension(3), intent(out) :: 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
    real(kind=realtype) :: arg1
! normalize the input vector.
    arg1 = refdirection(1)**2 + refdirection(2)**2 + refdirection(3)**2
    rnorm = sqrt(arg1)
    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)
      tmpd = -alphad
      tmp = -alpha
      zbnd = 0.0_8
      ybnd = 0.0_8
      xbnd = 0.0_8
      call vectorrotation_d(x1, x1d, y1, y1d, z1, z1d, 3, tmp, tmpd, xbn&
&                     , xbnd, ybn, ybnd, zbn, zbnd)
! 2) rotate beta radians ccw about y-axis
!    ( <=> rotate z-axis -beta radians ccw)
      tmpd = -betad
      tmp = -beta
      call vectorrotation_d(xw, xwd, yw, ywd, zw, zwd, 2, tmp, tmpd, x1&
&                     , x1d, y1, y1d, z1, z1d)
    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)
      zbnd = 0.0_8
      ybnd = 0.0_8
      xbnd = 0.0_8
      call vectorrotation_d(x1, x1d, y1, y1d, z1, z1d, 2, alpha, alphad&
&                     , xbn, xbnd, ybn, ybnd, zbn, zbnd)
! 2) rotate beta radians ccw about y-axis
!    ( <=> rotate z-axis -beta radians ccw)
      call vectorrotation_d(xw, xwd, yw, ywd, zw, zwd, 3, beta, betad, &
&                     x1, x1d, y1, y1d, z1, z1d)
    else
      call terminate('getdirvector', 'invalid lift direction')
      zwd = 0.0_8
      xwd = 0.0_8
      ywd = 0.0_8
    end if
    winddirectiond = 0.0_8
    winddirectiond(1) = xwd
    winddirection(1) = xw
    winddirectiond(2) = ywd
    winddirection(2) = yw
    winddirectiond(3) = zwd
    winddirection(3) = zw
  end subroutine getdirvector_d
 
  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_d, 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
    real(kind=realtype) :: arg1
! normalize the input vector.
    arg1 = refdirection(1)**2 + refdirection(2)**2 + refdirection(3)**2
    rnorm = sqrt(arg1)
    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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: xp yp zp
!   with respect to varying inputs: x y z angle
  subroutine vectorrotation_d(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), intent(in) :: angled, xd, yd, zd
    real(kind=realtype), intent(out) :: xp, yp, zp
    real(kind=realtype), intent(out) :: xpd, ypd, zpd
    intrinsic cos
    intrinsic sin
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
! rotation about specified axis by specified angle
    select case  (iaxis) 
    case (1) 
! rotation about the x-axis
      xpd = xd
      xp = 1.*x + 0.*y + 0.*z
      temp = cos(angle)
      temp0 = sin(angle)
      ypd = temp*yd - (y*sin(angle)-z*cos(angle))*angled + temp0*zd
      yp = 0.*x + temp*y + temp0*z
      temp0 = sin(angle)
      temp = cos(angle)
      zpd = temp*zd - temp0*yd - (y*cos(angle)+z*sin(angle))*angled
      zp = 0.*x - temp0*y + temp*z
! rotation about the y-axis
    case (2) 
      temp0 = cos(angle)
      temp = sin(angle)
      xpd = temp0*xd - (x*sin(angle)+z*cos(angle))*angled - temp*zd
      xp = temp0*x + 0.*y - temp*z
      ypd = yd
      yp = 0.*x + 1.*y + 0.*z
      temp0 = sin(angle)
      temp = cos(angle)
      zpd = (x*cos(angle)-z*sin(angle))*angled + temp0*xd + temp*zd
      zp = temp0*x + 0.*y + temp*z
! rotation about the z-axis
    case (3) 
      temp0 = cos(angle)
      temp = sin(angle)
      xpd = temp0*xd - (x*sin(angle)-y*cos(angle))*angled + temp*yd
      xp = temp0*x + temp*y + 0.*z
      temp0 = cos(angle)
      temp = sin(angle)
      ypd = temp0*yd - (y*sin(angle)+x*cos(angle))*angled - temp*xd
      yp = temp0*y - temp*x + 0.*z
      zpd = zd
      zp = 0.*x + 0.*y + 1.*z
    case default
      write(*, *) 'vectorrotation called with invalid arguments'
      stop
    end select
  end subroutine vectorrotation_d
 
  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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *wx *wy *wz *qx *qy *qz *ux
!                *uy *uz *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:in *wx:in-out *wy:in-out *wz:in-out
!                *w:in *qx:in-out *qy:in-out *qz:in-out *vol:in
!                *ux:in-out *uy:in-out *uz:in-out *si:in *sj:in
!                *sk:in *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_d()
!
!         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
    real(kind=realtype) :: temp
! zero all nodeal gradients:
    uxd = 0.0_8
    ux = zero
    uyd = 0.0_8
    uy = zero
    uzd = 0.0_8
    uz = zero
    vxd = 0.0_8
    vx = zero
    vyd = 0.0_8
    vy = zero
    vzd = 0.0_8
    vz = zero
    wxd = 0.0_8
    wx = zero
    wyd = 0.0_8
    wy = zero
    wzd = 0.0_8
    wz = zero
    qxd = 0.0_8
    qx = zero
    qyd = 0.0_8
    qy = zero
    qzd = 0.0_8
    qz = zero
! first part. contribution in the k-direction.
! the contribution is scattered to both the left and right node
! in k-direction.
    do k=1,ke
      do j=1,jl
        do i=1,il
! 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.
          sxd = skd(i, j, k-1, 1) + skd(i+1, j, k-1, 1) + skd(i, j+1, k-&
&           1, 1) + skd(i+1, j+1, k-1, 1) + skd(i, j, k, 1) + skd(i+1, j&
&           , k, 1) + skd(i, j+1, k, 1) + skd(i+1, j+1, k, 1)
          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)
          syd = skd(i, j, k-1, 2) + skd(i+1, j, k-1, 2) + skd(i, j+1, k-&
&           1, 2) + skd(i+1, j+1, k-1, 2) + skd(i, j, k, 2) + skd(i+1, j&
&           , k, 2) + skd(i, j+1, k, 2) + skd(i+1, j+1, k, 2)
          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)
          szd = skd(i, j, k-1, 3) + skd(i+1, j, k-1, 3) + skd(i, j+1, k-&
&           1, 3) + skd(i+1, j+1, k-1, 3) + skd(i, j, k, 3) + skd(i+1, j&
&           , k, 3) + skd(i, j+1, k, 3) + skd(i+1, j+1, k, 3)
          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.
          ubard = fourth*(wd(i, j, k, ivx)+wd(i+1, j, k, ivx)+wd(i, j+1&
&           , k, ivx)+wd(i+1, j+1, k, ivx))
          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))
          vbard = fourth*(wd(i, j, k, ivy)+wd(i+1, j, k, ivy)+wd(i, j+1&
&           , k, ivy)+wd(i+1, j+1, k, ivy))
          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))
          wbard = fourth*(wd(i, j, k, ivz)+wd(i+1, j, k, ivz)+wd(i, j+1&
&           , k, ivz)+wd(i+1, j+1, k, ivz))
          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))
          a2d = fourth*(aad(i, j, k)+aad(i+1, j, k)+aad(i, j+1, k)+aad(i&
&           +1, j+1, k))
          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
            uxd(i, j, k-1) = uxd(i, j, k-1) + sx*ubard + ubar*sxd
            ux(i, j, k-1) = ux(i, j, k-1) + ubar*sx
            uyd(i, j, k-1) = uyd(i, j, k-1) + sy*ubard + ubar*syd
            uy(i, j, k-1) = uy(i, j, k-1) + ubar*sy
            uzd(i, j, k-1) = uzd(i, j, k-1) + sz*ubard + ubar*szd
            uz(i, j, k-1) = uz(i, j, k-1) + ubar*sz
            vxd(i, j, k-1) = vxd(i, j, k-1) + sx*vbard + vbar*sxd
            vx(i, j, k-1) = vx(i, j, k-1) + vbar*sx
            vyd(i, j, k-1) = vyd(i, j, k-1) + sy*vbard + vbar*syd
            vy(i, j, k-1) = vy(i, j, k-1) + vbar*sy
            vzd(i, j, k-1) = vzd(i, j, k-1) + sz*vbard + vbar*szd
            vz(i, j, k-1) = vz(i, j, k-1) + vbar*sz
            wxd(i, j, k-1) = wxd(i, j, k-1) + sx*wbard + wbar*sxd
            wx(i, j, k-1) = wx(i, j, k-1) + wbar*sx
            wyd(i, j, k-1) = wyd(i, j, k-1) + sy*wbard + wbar*syd
            wy(i, j, k-1) = wy(i, j, k-1) + wbar*sy
            wzd(i, j, k-1) = wzd(i, j, k-1) + sz*wbard + wbar*szd
            wz(i, j, k-1) = wz(i, j, k-1) + wbar*sz
            qxd(i, j, k-1) = qxd(i, j, k-1) - sx*a2d - a2*sxd
            qx(i, j, k-1) = qx(i, j, k-1) - a2*sx
            qyd(i, j, k-1) = qyd(i, j, k-1) - sy*a2d - a2*syd
            qy(i, j, k-1) = qy(i, j, k-1) - a2*sy
            qzd(i, j, k-1) = qzd(i, j, k-1) - sz*a2d - a2*szd
            qz(i, j, k-1) = qz(i, j, k-1) - a2*sz
          end if
          if (k .lt. ke) then
            uxd(i, j, k) = uxd(i, j, k) - sx*ubard - ubar*sxd
            ux(i, j, k) = ux(i, j, k) - ubar*sx
            uyd(i, j, k) = uyd(i, j, k) - sy*ubard - ubar*syd
            uy(i, j, k) = uy(i, j, k) - ubar*sy
            uzd(i, j, k) = uzd(i, j, k) - sz*ubard - ubar*szd
            uz(i, j, k) = uz(i, j, k) - ubar*sz
            vxd(i, j, k) = vxd(i, j, k) - sx*vbard - vbar*sxd
            vx(i, j, k) = vx(i, j, k) - vbar*sx
            vyd(i, j, k) = vyd(i, j, k) - sy*vbard - vbar*syd
            vy(i, j, k) = vy(i, j, k) - vbar*sy
            vzd(i, j, k) = vzd(i, j, k) - sz*vbard - vbar*szd
            vz(i, j, k) = vz(i, j, k) - vbar*sz
            wxd(i, j, k) = wxd(i, j, k) - sx*wbard - wbar*sxd
            wx(i, j, k) = wx(i, j, k) - wbar*sx
            wyd(i, j, k) = wyd(i, j, k) - sy*wbard - wbar*syd
            wy(i, j, k) = wy(i, j, k) - wbar*sy
            wzd(i, j, k) = wzd(i, j, k) - sz*wbard - wbar*szd
            wz(i, j, k) = wz(i, j, k) - wbar*sz
            qxd(i, j, k) = qxd(i, j, k) + sx*a2d + a2*sxd
            qx(i, j, k) = qx(i, j, k) + a2*sx
            qyd(i, j, k) = qyd(i, j, k) + sy*a2d + a2*syd
            qy(i, j, k) = qy(i, j, k) + a2*sy
            qzd(i, j, k) = qzd(i, j, k) + sz*a2d + a2*szd
            qz(i, j, k) = qz(i, j, k) + a2*sz
          end if
        end do
      end do
    end do
! second part. contribution in the j-direction.
! the contribution is scattered to both the left and right node
! in j-direction.
    do k=1,kl
      do j=1,je
        do i=1,il
! 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.
          sxd = sjd(i, j-1, k, 1) + sjd(i+1, j-1, k, 1) + sjd(i, j-1, k+&
&           1, 1) + sjd(i+1, j-1, k+1, 1) + sjd(i, j, k, 1) + sjd(i+1, j&
&           , k, 1) + sjd(i, j, k+1, 1) + sjd(i+1, j, k+1, 1)
          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)
          syd = sjd(i, j-1, k, 2) + sjd(i+1, j-1, k, 2) + sjd(i, j-1, k+&
&           1, 2) + sjd(i+1, j-1, k+1, 2) + sjd(i, j, k, 2) + sjd(i+1, j&
&           , k, 2) + sjd(i, j, k+1, 2) + sjd(i+1, j, k+1, 2)
          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)
          szd = sjd(i, j-1, k, 3) + sjd(i+1, j-1, k, 3) + sjd(i, j-1, k+&
&           1, 3) + sjd(i+1, j-1, k+1, 3) + sjd(i, j, k, 3) + sjd(i+1, j&
&           , k, 3) + sjd(i, j, k+1, 3) + sjd(i+1, j, k+1, 3)
          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.
          ubard = fourth*(wd(i, j, k, ivx)+wd(i+1, j, k, ivx)+wd(i, j, k&
&           +1, ivx)+wd(i+1, j, k+1, ivx))
          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))
          vbard = fourth*(wd(i, j, k, ivy)+wd(i+1, j, k, ivy)+wd(i, j, k&
&           +1, ivy)+wd(i+1, j, k+1, ivy))
          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))
          wbard = fourth*(wd(i, j, k, ivz)+wd(i+1, j, k, ivz)+wd(i, j, k&
&           +1, ivz)+wd(i+1, j, k+1, ivz))
          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))
          a2d = fourth*(aad(i, j, k)+aad(i+1, j, k)+aad(i, j, k+1)+aad(i&
&           +1, j, k+1))
          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
            uxd(i, j-1, k) = uxd(i, j-1, k) + sx*ubard + ubar*sxd
            ux(i, j-1, k) = ux(i, j-1, k) + ubar*sx
            uyd(i, j-1, k) = uyd(i, j-1, k) + sy*ubard + ubar*syd
            uy(i, j-1, k) = uy(i, j-1, k) + ubar*sy
            uzd(i, j-1, k) = uzd(i, j-1, k) + sz*ubard + ubar*szd
            uz(i, j-1, k) = uz(i, j-1, k) + ubar*sz
            vxd(i, j-1, k) = vxd(i, j-1, k) + sx*vbard + vbar*sxd
            vx(i, j-1, k) = vx(i, j-1, k) + vbar*sx
            vyd(i, j-1, k) = vyd(i, j-1, k) + sy*vbard + vbar*syd
            vy(i, j-1, k) = vy(i, j-1, k) + vbar*sy
            vzd(i, j-1, k) = vzd(i, j-1, k) + sz*vbard + vbar*szd
            vz(i, j-1, k) = vz(i, j-1, k) + vbar*sz
            wxd(i, j-1, k) = wxd(i, j-1, k) + sx*wbard + wbar*sxd
            wx(i, j-1, k) = wx(i, j-1, k) + wbar*sx
            wyd(i, j-1, k) = wyd(i, j-1, k) + sy*wbard + wbar*syd
            wy(i, j-1, k) = wy(i, j-1, k) + wbar*sy
            wzd(i, j-1, k) = wzd(i, j-1, k) + sz*wbard + wbar*szd
            wz(i, j-1, k) = wz(i, j-1, k) + wbar*sz
            qxd(i, j-1, k) = qxd(i, j-1, k) - sx*a2d - a2*sxd
            qx(i, j-1, k) = qx(i, j-1, k) - a2*sx
            qyd(i, j-1, k) = qyd(i, j-1, k) - sy*a2d - a2*syd
            qy(i, j-1, k) = qy(i, j-1, k) - a2*sy
            qzd(i, j-1, k) = qzd(i, j-1, k) - sz*a2d - a2*szd
            qz(i, j-1, k) = qz(i, j-1, k) - a2*sz
          end if
          if (j .lt. je) then
            uxd(i, j, k) = uxd(i, j, k) - sx*ubard - ubar*sxd
            ux(i, j, k) = ux(i, j, k) - ubar*sx
            uyd(i, j, k) = uyd(i, j, k) - sy*ubard - ubar*syd
            uy(i, j, k) = uy(i, j, k) - ubar*sy
            uzd(i, j, k) = uzd(i, j, k) - sz*ubard - ubar*szd
            uz(i, j, k) = uz(i, j, k) - ubar*sz
            vxd(i, j, k) = vxd(i, j, k) - sx*vbard - vbar*sxd
            vx(i, j, k) = vx(i, j, k) - vbar*sx
            vyd(i, j, k) = vyd(i, j, k) - sy*vbard - vbar*syd
            vy(i, j, k) = vy(i, j, k) - vbar*sy
            vzd(i, j, k) = vzd(i, j, k) - sz*vbard - vbar*szd
            vz(i, j, k) = vz(i, j, k) - vbar*sz
            wxd(i, j, k) = wxd(i, j, k) - sx*wbard - wbar*sxd
            wx(i, j, k) = wx(i, j, k) - wbar*sx
            wyd(i, j, k) = wyd(i, j, k) - sy*wbard - wbar*syd
            wy(i, j, k) = wy(i, j, k) - wbar*sy
            wzd(i, j, k) = wzd(i, j, k) - sz*wbard - wbar*szd
            wz(i, j, k) = wz(i, j, k) - wbar*sz
            qxd(i, j, k) = qxd(i, j, k) + sx*a2d + a2*sxd
            qx(i, j, k) = qx(i, j, k) + a2*sx
            qyd(i, j, k) = qyd(i, j, k) + sy*a2d + a2*syd
            qy(i, j, k) = qy(i, j, k) + a2*sy
            qzd(i, j, k) = qzd(i, j, k) + sz*a2d + a2*szd
            qz(i, j, k) = qz(i, j, k) + a2*sz
          end if
        end do
      end do
    end do
! third part. contribution in the i-direction.
! the contribution is scattered to both the left and right node
! in i-direction.
    do k=1,kl
      do j=1,jl
        do i=1,ie
! 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.
          sxd = sid(i-1, j, k, 1) + sid(i-1, j+1, k, 1) + sid(i-1, j, k+&
&           1, 1) + sid(i-1, j+1, k+1, 1) + sid(i, j, k, 1) + sid(i, j+1&
&           , k, 1) + sid(i, j, k+1, 1) + sid(i, j+1, k+1, 1)
          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)
          syd = sid(i-1, j, k, 2) + sid(i-1, j+1, k, 2) + sid(i-1, j, k+&
&           1, 2) + sid(i-1, j+1, k+1, 2) + sid(i, j, k, 2) + sid(i, j+1&
&           , k, 2) + sid(i, j, k+1, 2) + sid(i, j+1, k+1, 2)
          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)
          szd = sid(i-1, j, k, 3) + sid(i-1, j+1, k, 3) + sid(i-1, j, k+&
&           1, 3) + sid(i-1, j+1, k+1, 3) + sid(i, j, k, 3) + sid(i, j+1&
&           , k, 3) + sid(i, j, k+1, 3) + sid(i, j+1, k+1, 3)
          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.
          ubard = fourth*(wd(i, j, k, ivx)+wd(i, j+1, k, ivx)+wd(i, j, k&
&           +1, ivx)+wd(i, j+1, k+1, ivx))
          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))
          vbard = fourth*(wd(i, j, k, ivy)+wd(i, j+1, k, ivy)+wd(i, j, k&
&           +1, ivy)+wd(i, j+1, k+1, ivy))
          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))
          wbard = fourth*(wd(i, j, k, ivz)+wd(i, j+1, k, ivz)+wd(i, j, k&
&           +1, ivz)+wd(i, j+1, k+1, ivz))
          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))
          a2d = fourth*(aad(i, j, k)+aad(i, j+1, k)+aad(i, j, k+1)+aad(i&
&           , j+1, k+1))
          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
            uxd(i-1, j, k) = uxd(i-1, j, k) + sx*ubard + ubar*sxd
            ux(i-1, j, k) = ux(i-1, j, k) + ubar*sx
            uyd(i-1, j, k) = uyd(i-1, j, k) + sy*ubard + ubar*syd
            uy(i-1, j, k) = uy(i-1, j, k) + ubar*sy
            uzd(i-1, j, k) = uzd(i-1, j, k) + sz*ubard + ubar*szd
            uz(i-1, j, k) = uz(i-1, j, k) + ubar*sz
            vxd(i-1, j, k) = vxd(i-1, j, k) + sx*vbard + vbar*sxd
            vx(i-1, j, k) = vx(i-1, j, k) + vbar*sx
            vyd(i-1, j, k) = vyd(i-1, j, k) + sy*vbard + vbar*syd
            vy(i-1, j, k) = vy(i-1, j, k) + vbar*sy
            vzd(i-1, j, k) = vzd(i-1, j, k) + sz*vbard + vbar*szd
            vz(i-1, j, k) = vz(i-1, j, k) + vbar*sz
            wxd(i-1, j, k) = wxd(i-1, j, k) + sx*wbard + wbar*sxd
            wx(i-1, j, k) = wx(i-1, j, k) + wbar*sx
            wyd(i-1, j, k) = wyd(i-1, j, k) + sy*wbard + wbar*syd
            wy(i-1, j, k) = wy(i-1, j, k) + wbar*sy
            wzd(i-1, j, k) = wzd(i-1, j, k) + sz*wbard + wbar*szd
            wz(i-1, j, k) = wz(i-1, j, k) + wbar*sz
            qxd(i-1, j, k) = qxd(i-1, j, k) - sx*a2d - a2*sxd
            qx(i-1, j, k) = qx(i-1, j, k) - a2*sx
            qyd(i-1, j, k) = qyd(i-1, j, k) - sy*a2d - a2*syd
            qy(i-1, j, k) = qy(i-1, j, k) - a2*sy
            qzd(i-1, j, k) = qzd(i-1, j, k) - sz*a2d - a2*szd
            qz(i-1, j, k) = qz(i-1, j, k) - a2*sz
          end if
          if (i .lt. ie) then
            uxd(i, j, k) = uxd(i, j, k) - sx*ubard - ubar*sxd
            ux(i, j, k) = ux(i, j, k) - ubar*sx
            uyd(i, j, k) = uyd(i, j, k) - sy*ubard - ubar*syd
            uy(i, j, k) = uy(i, j, k) - ubar*sy
            uzd(i, j, k) = uzd(i, j, k) - sz*ubard - ubar*szd
            uz(i, j, k) = uz(i, j, k) - ubar*sz
            vxd(i, j, k) = vxd(i, j, k) - sx*vbard - vbar*sxd
            vx(i, j, k) = vx(i, j, k) - vbar*sx
            vyd(i, j, k) = vyd(i, j, k) - sy*vbard - vbar*syd
            vy(i, j, k) = vy(i, j, k) - vbar*sy
            vzd(i, j, k) = vzd(i, j, k) - sz*vbard - vbar*szd
            vz(i, j, k) = vz(i, j, k) - vbar*sz
            wxd(i, j, k) = wxd(i, j, k) - sx*wbard - wbar*sxd
            wx(i, j, k) = wx(i, j, k) - wbar*sx
            wyd(i, j, k) = wyd(i, j, k) - sy*wbard - wbar*syd
            wy(i, j, k) = wy(i, j, k) - wbar*sy
            wzd(i, j, k) = wzd(i, j, k) - sz*wbard - wbar*szd
            wz(i, j, k) = wz(i, j, k) - wbar*sz
            qxd(i, j, k) = qxd(i, j, k) + sx*a2d + a2*sxd
            qx(i, j, k) = qx(i, j, k) + a2*sx
            qyd(i, j, k) = qyd(i, j, k) + sy*a2d + a2*syd
            qy(i, j, k) = qy(i, j, k) + a2*sy
            qzd(i, j, k) = qzd(i, j, k) + sz*a2d + a2*szd
            qz(i, j, k) = qz(i, j, k) + a2*sz
          end if
        end do
      end do
    end do
! divide by 8 times the volume to obtain the correct gradients.
    do k=1,kl
      do j=1,jl
        do i=1,il
! compute the inverse of 8 times the volume for this node.
          temp = 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))
          oneovervd = -(temp*(vold(i, j, k)+vold(i, j, k+1)+vold(i+1, j&
&           , k)+vold(i+1, j, k+1)+vold(i, j+1, k)+vold(i, j+1, k+1)+&
&           vold(i+1, j+1, k)+vold(i+1, j+1, k+1))/(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)))
          oneoverv = temp
! compute the correct velocity gradients and "unit" heat
! fluxes. the velocity gradients are stored in ux, etc.
          uxd(i, j, k) = oneoverv*uxd(i, j, k) + ux(i, j, k)*oneovervd
          ux(i, j, k) = ux(i, j, k)*oneoverv
          uyd(i, j, k) = oneoverv*uyd(i, j, k) + uy(i, j, k)*oneovervd
          uy(i, j, k) = uy(i, j, k)*oneoverv
          uzd(i, j, k) = oneoverv*uzd(i, j, k) + uz(i, j, k)*oneovervd
          uz(i, j, k) = uz(i, j, k)*oneoverv
          vxd(i, j, k) = oneoverv*vxd(i, j, k) + vx(i, j, k)*oneovervd
          vx(i, j, k) = vx(i, j, k)*oneoverv
          vyd(i, j, k) = oneoverv*vyd(i, j, k) + vy(i, j, k)*oneovervd
          vy(i, j, k) = vy(i, j, k)*oneoverv
          vzd(i, j, k) = oneoverv*vzd(i, j, k) + vz(i, j, k)*oneovervd
          vz(i, j, k) = vz(i, j, k)*oneoverv
          wxd(i, j, k) = oneoverv*wxd(i, j, k) + wx(i, j, k)*oneovervd
          wx(i, j, k) = wx(i, j, k)*oneoverv
          wyd(i, j, k) = oneoverv*wyd(i, j, k) + wy(i, j, k)*oneovervd
          wy(i, j, k) = wy(i, j, k)*oneoverv
          wzd(i, j, k) = oneoverv*wzd(i, j, k) + wz(i, j, k)*oneovervd
          wz(i, j, k) = wz(i, j, k)*oneoverv
          qxd(i, j, k) = oneoverv*qxd(i, j, k) + qx(i, j, k)*oneovervd
          qx(i, j, k) = qx(i, j, k)*oneoverv
          qyd(i, j, k) = oneoverv*qyd(i, j, k) + qy(i, j, k)*oneovervd
          qy(i, j, k) = qy(i, j, k)*oneoverv
          qzd(i, j, k) = oneoverv*qzd(i, j, k) + qz(i, j, k)*oneovervd
          qz(i, j, k) = qz(i, j, k)*oneoverv
        end do
      end do
    end do
  end subroutine allnodalgradients_d
 
  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
! 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
! first part. contribution in the k-direction.
! the contribution is scattered to both the left and right node
! in k-direction.
    do k=1,ke
      do j=1,jl
        do i=1,il
! 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
      end do
    end do
! second part. contribution in the j-direction.
! the contribution is scattered to both the left and right node
! in j-direction.
    do k=1,kl
      do j=1,je
        do i=1,il
! 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
      end do
    end do
! third part. contribution in the i-direction.
! the contribution is scattered to both the left and right node
! in i-direction.
    do k=1,kl
      do j=1,jl
        do i=1,ie
! 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
      end do
    end do
! divide by 8 times the volume to obtain the correct gradients.
    do k=1,kl
      do j=1,jl
        do i=1,il
! 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 do
    end do
  end subroutine allnodalgradients
! ----------------------------------------------------------------------
!                                                                      |
!                    no tapenade routine below this line               |
!                                                                      |
! ----------------------------------------------------------------------
 
end module flowutils_d