flowUtils.F90

Go to the documentation of this file.

module flowutils
contains
 
    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, 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
 
            !        ================================================================
#ifndef USE_TAPENADE
        case (cptempcurvefits)
 
            ! Cp as function of the temperature is given via curve fits.
 
            do i = 1, mm
 
                ! Determine the case we are having here.
 
                if (t(i) <= cptrange(0)) then
 
                    ! Temperature is less than the smallest temperature
                    ! in the curve fits. Use cv0 to compute gamma.
 
                    gamma(i) = (cv0 + one) / cv0
 
                else if (t(i) >= cptrange(cpnparts)) then
 
                    ! Temperature is larger than the largest temperature
                    ! in the curve fits. Use cvn to compute gamma.
 
                    gamma(i) = (cvn + one) / cvn
 
                else
 
                    ! Temperature is in the curve fit range.
                    ! First find the valid range.
 
                    ii = cpnparts
                    start = 1
                    interval: do
 
                        ! Next guess for the interval.
 
                        nn = start + ii / 2
 
                        ! Determine the situation we are having here.
 
                        if (t(i) > cptrange(nn)) then
 
                            ! Temperature is larger than the upper boundary of
                            ! the current interval. Update the lower boundary.
 
                            start = nn + 1
                            ii = ii - 1
 
                        else if (t(i) >= cptrange(nn - 1)) then
 
                            ! This is the correct range. Exit the do-loop.
 
                            exit
 
                        end if
 
                        ! Modify ii for the next branch to search.
 
                        ii = ii / 2
 
                    end do interval
 
                    ! Nn contains the correct curve fit interval.
                    ! Compute the value of cp.
 
                    cp = zero
                    do ii = 1, cptempfit(nn)%nterm
                        t2 = t(i)**(cptempfit(nn)%exponents(ii))
                        cp = cp + cptempfit(nn)%constants(ii) * t2
                    end do
 
                    ! Compute the corresponding value of gamma.
 
                    gamma(i) = cp / (cp - one)
 
                end if
 
            end do
#endif
        end select
 
    end subroutine computegamma
 
    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)
 
            !===============================================================
#ifndef USE_TAPENADE
        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.
 
            ! Compute the dimensional temperature and the scale
            ! factor to convert the integral of cp to the correct
            ! nonDimensional value.
 
            t = tref * p / (rgas * rho)
            scale = rgas / tref
 
            ! Compute the enthalpy and the integrand of cp/(r*t) at the
            ! given temperature. Take care of the exceptional situations.
 
            if (t <= cptrange(0)) then
 
                ! Temperature is smaller than the smallest temperature in
                ! the curve fit range. Use extrapolation using constant cp.
 
                nn = 0
                cp = cv0 + one
                h = scale * (cphint(0) + cp * (t - cptrange(0)))
 
                intcportt = cp * log(t)
 
            else if (t >= cptrange(cpnparts)) then
 
                ! Temperature is larger than the largest temperature in the
                ! curve fit range. Use extrapolation using constant cp.
 
                nn = cpnparts + 1
                cp = cvn + one
                h = scale * (cphint(cpnparts) + cp * (t - cptrange(cpnparts)))
 
                intcportt = cp * log(t)
 
            else
 
                ! Temperature lies in the curve fit range. Find the correct
                ! interval.
 
                ii = cpnparts
                start = 1
                interval: do
 
                    ! Next guess for the interval.
 
                    nn = start + ii / 2
 
                    ! Determine the situation we are having here.
 
                    if (t > cptrange(nn)) then
 
                        ! Temperature is larger than the upper boundary of
                        ! the current interval. Update the lower boundary.
 
                        start = nn + 1
                        ii = ii - 1
 
                    else if (t >= cptrange(nn - 1)) then
 
                        ! This is the correct range. Exit the do-loop.
 
                        exit
 
                    end if
 
                    ! Modify ii for the next branch to search.
 
                    ii = ii / 2
 
                end do interval
 
                ! nn contains the correct curve fit interval.
                ! Integrate cp to compute h and the integrand of cp/(r*t)
 
                h = cptempfit(nn)%eint0
                intcportt = zero
 
                do ii = 1, cptempfit(nn)%nterm
                    if (cptempfit(nn)%exponents(ii) == -1_inttype) then
                        h = h + cptempfit(nn)%constants(ii) * log(t)
                    else
                        mm = cptempfit(nn)%exponents(ii) + 1
                        t2 = t**mm
                        h = h + cptempfit(nn)%constants(ii) * t2 / mm
                    end if
 
                    if (cptempfit(nn)%exponents(ii) == 0_inttype) then
                        intcportt = intcportt &
                                    + cptempfit(nn)%constants(ii) * log(t)
                    else
                        mm = cptempfit(nn)%exponents(ii)
                        t2 = t**mm
                        intcportt = intcportt &
                                    + cptempfit(nn)%constants(ii) * t2 / mm
                    end if
                end do
 
                h = scale * h
 
            end if
 
            ! Compute the total enthalpy. Divide by scale to get the same
            ! dimensions as for the integral of cp/r.
 
            htot = (h + half * (u * u + v * v + w * w)) / scale
 
            ! Compute the corresponding total temperature. First determine
            ! the situation we are having here.
 
            if (htot <= cphint(0)) then
 
                ! Total enthalpy is smaller than the lowest value of the
                ! curve fit. Use extrapolation using constant cp.
 
                nnt = 0
                tt = cptrange(0) + (htot - cphint(0)) / (cv0 + one)
 
            else if (htot >= cphint(cpnparts)) then
 
                ! Total enthalpy is larger than the largest value of the
                ! curve fit. Use extrapolation using constant cp.
 
                nnt = cpnparts + 1
                tt = cptrange(cpnparts) &
                     + (htot - cphint(cpnparts)) / (cvn + one)
 
            else
 
                ! Total temperature is in the range of the curve fits.
                ! Use a newton algorithm to find the correct temperature.
                ! First find the correct interval.
 
                ii = cpnparts
                start = 1
                intervaltt: do
 
                    ! Next guess for the interval.
 
                    nnt = start + ii / 2
 
                    ! Determine the situation we are having here.
 
                    if (htot > cphint(nnt)) then
 
                        ! Enthalpy is larger than the upper boundary of
                        ! the current interval. Update the lower boundary.
 
                        start = nnt + 1
                        ii = ii - 1
 
                    else if (htot >= cphint(nnt - 1)) then
 
                        ! This is the correct range. Exit the do-loop.
 
                        exit
 
                    end if
 
                    ! Modify ii for the next branch to search.
 
                    ii = ii / 2
 
                end do intervaltt
 
                ! Nnt contains the range in which the newton algorithm must
                ! be applied. Initial guess of the total temperature.
 
                alp = (cphint(nnt) - htot) / (cphint(nnt) - cphint(nnt - 1))
                tt = alp * cptrange(nnt - 1) + (one - alp) * cptrange(nnt)
 
                ! The actual newton algorithm to compute the total
                ! temperature.
 
                newton: do
 
                    ! Compute the energy as well as the value of cv/r for the
                    ! given temperature.
 
                    cp = zero
                    h = cptempfit(nnt)%eint0
 
                    do ii = 1, cptempfit(nnt)%nterm
 
                        ! Update cp.
 
                        t2 = tt**(cptempfit(nnt)%exponents(ii))
                        cp = cp + cptempfit(nnt)%constants(ii) * t2
 
                        ! Update h, for which this contribution must be
                        ! integrated. Take the exceptional case that the
                        ! exponent == -1 into account.
 
                        if (cptempfit(nnt)%exponents(ii) == -1_inttype) then
                            h = h + cptempfit(nnt)%constants(ii) * log(tt)
                        else
                            h = h + cptempfit(nnt)%constants(ii) * t2 * tt &
                                / (cptempfit(nnt)%exponents(ii) + 1)
                        end if
 
                    end do
 
                    ! Compute the update and the new total temperature.
 
                    dt = (htot - h) / cp
                    tt = tt + dt
 
                    ! Exit the newton loop if the update is smaller than the
                    ! threshold value.
 
                    if (abs(dt) < dtstop) exit
 
                end do newton
 
            end if
 
            ! To compute the total pressure, the integral of cp/(r*T)
            ! must be computed from T = T to T = Tt. Compute the integrand
            ! at T = Tt; take care of the exceptional situations.
 
            if (tt <= cptrange(0)) then
                intcporttt = (cv0 + one) * log(tt)
            else if (tt >= cptrange(cpnparts)) then
                intcporttt = (cvn + one) * log(tt)
            else
 
                intcporttt = zero
                do ii = 1, cptempfit(nnt)%nterm
                    if (cptempfit(nnt)%exponents(ii) == 0_inttype) then
                        intcporttt = intcporttt &
                                     + cptempfit(nnt)%constants(ii) * log(tt)
                    else
                        mm = cptempfit(nnt)%exponents(ii)
                        t2 = tt**mm
                        intcporttt = intcporttt &
                                     + cptempfit(nnt)%constants(ii) * t2 / mm
                    end if
                end do
 
            end if
 
            ! Compute the integral of cp/(r*T) from T to Tt. First
            ! substract the lower boundary from the upper boundary.
 
            intcport = intcporttt - intcportt
 
            ! Add the contributions from the possible internal curve fit
            ! boundaries if Tt and T are in different curve fit intervals.
 
            do mm = (nn + 1), nnt
                ii = mm - 1
 
                if (ii == 0_inttype) then
                    intcport = intcport + (cv0 + one) * log(cptrange(0))
                else
                    intcport = intcport + cptempfit(ii)%intCpovrt_2
                end if
 
                if (mm > cpnparts) then
                    intcport = intcport - (cvn + one) * log(cptrange(cpnparts))
                else
                    intcport = intcport - cptempfit(mm)%intCpovrt_1
                end if
            end do
 
            ! And finally, compute the total pressure.
 
            ptot = p * exp(intcport)
#endif
        end select
 
    end subroutine computeptot
 
    subroutine computespeedofsoundsquared
 
        !
        !       computeSpeedOfSoundSquared does what it says.
        !
        use constants
        use blockpointers, only: ie, je, ke, w, p, aa, gamma
        use utils, 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
#ifdef TAPENADE_REVERSE
            !$AD II-LOOP
            do ii = 0, ie * je * ke - 1
                i = mod(ii, ie) + 1
                j = mod(ii / ie, je) + 1
                k = ii / (ie * je) + 1
#else
                do k = 1, ke
                    do j = 1, je
                        do i = 1, ie
#endif
                            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)
#ifdef TAPENADE_REVERSE
                        end do
#else
                    end do
                end do
            end do
#endif
        else
#ifdef TAPENADE_REVERSE
            !$AD II-LOOP
            do ii = 0, ie * je * ke - 1
                i = mod(ii, ie) + 1
                j = mod(ii / ie, je) + 1
                k = ii / (ie * je) + 1
#else
                do k = 1, ke
                    do j = 1, je
                        do i = 1, ie
#endif
                            aa(i, j, k) = gamma(i, j, k) * p(i, j, k) / w(i, j, k, irho)
#ifdef TAPENADE_REVERSE
                        end do
#else
                    end do
                end do
            end do
#endif
        end if
    end subroutine computespeedofsoundsquared
    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
#ifdef TAPENADE_REVERSE
                !$AD II-LOOP
                do ii = 0, isize * jsize * ksize - 1
                    i = mod(ii, isize) + istart
                    j = mod(ii / (isize), jsize) + jstart
                    k = ii / ((isize * jsize)) + kstart
#else
                    do k = kstart, kend
                        do j = jstart, jend
                            do i = istart, iend
#endif
                                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)
 
#ifdef TAPENADE_REVERSE
                            end do
#else
                        end do
                    end do
                end do
#endif
 
            else
#ifdef TAPENADE_REVERSE
                !$AD II-LOOP
                do ii = 0, isize * jsize * ksize - 1
                    i = mod(ii, isize) + istart
                    j = mod(ii / (isize), jsize) + jstart
                    k = ii / ((isize * jsize)) + kstart
#else
                    do k = kstart, kend
                        do j = jstart, jend
                            do i = istart, iend
#endif
                                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)
#ifdef TAPENADE_REVERSE
                            end do
#else
                        end do
                    end do
                end do
#endif
            end if
 
#ifndef USE_TAPENADE
        case (cptempcurvefits)
 
            ! Cp as function of the temperature is given via curve fits.
 
            ! Store a scale factor to compute the nonDimensional
            ! internal energy.
 
            scale = rgas / tref
 
            ! Loop over the given range of the block.
 
            do k = kstart, kend
                do j = jstart, jend
                    do i = istart, iend
                        call computeetotcellcpfit(i, j, k, scale, &
                                                  correctfork)
                    end do
                end do
            end do
#endif
 
        end select
    end subroutine computeetotblock
 
    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
 
    !      ==================================================================
 
    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
 
#ifndef USE_TAPENADE
        case (cptempcurvefits)
 
            ! Cp as function of the temperature is given via curve fits.
 
            ! Store a scale factor to compute the nonDimensional
            ! internal energy.
 
            scale = rgas / tref
 
            ! Loop over the number of elements of the array
 
            ! Compute the dimensional temperature.
 
            pp = p
            if (correctfork) pp = pp - twothird * rho * k
            t = tref * pp / (rgas * rho)
 
            ! Determine the case we are having here.
 
            if (t <= cptrange(0)) then
 
                ! Temperature is less than the smallest temperature
                ! in the curve fits. Use extrapolation using
                ! constant cv.
 
                einternal = scale * (cpeint(0) + cv0 * (t - cptrange(0)))
 
            else if (t >= cptrange(cpnparts)) then
 
                ! Temperature is larger than the largest temperature
                ! in the curve fits. Use extrapolation using
                ! constant cv.
 
                einternal = scale * (cpeint(cpnparts) &
                                     + cvn * (t - cptrange(cpnparts)))
 
            else
 
                ! Temperature is in the curve fit range.
                ! First find the valid range.
 
                ii = cpnparts
                start = 1
                interval: do
 
                    ! Next guess for the interval.
 
                    nn = start + ii / 2
 
                    ! Determine the situation we are having here.
 
                    if (t > cptrange(nn)) then
 
                        ! Temperature is larger than the upper boundary of
                        ! the current interval. Update the lower boundary.
 
                        start = nn + 1
                        ii = ii - 1
 
                    else if (t >= cptrange(nn - 1)) then
 
                        ! This is the correct range. Exit the do-loop.
 
                        exit
 
                    end if
 
                    ! Modify ii for the next branch to search.
 
                    ii = ii / 2
 
                end do interval
 
                ! Nn contains the correct curve fit interval.
                ! Integrate cv to compute eInternal.
 
                einternal = cptempfit(nn)%eint0 - t
                do ii = 1, cptempfit(nn)%nterm
                    if (cptempfit(nn)%exponents(ii) == -1) then
                        einternal = einternal &
                                    + cptempfit(nn)%constants(ii) * log(t)
                    else
                        mm = cptempfit(nn)%exponents(ii) + 1
                        t2 = t**mm
                        einternal = einternal &
                                    + cptempfit(nn)%constants(ii) * t2 / mm
                    end if
                end do
 
                einternal = scale * einternal
 
            end if
 
            ! Add the turbulent energy if needed.
 
            if (correctfork) einternal = einternal + k
 
#endif
        end select
 
    end subroutine eint
 
    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
        ! 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
 
#ifdef TAPENADE_REVERSE
        isize = (iend - ibeg) + 1
        jsize = (jend - jbeg) + 1
        ksize = (kend - kbeg) + 1
 
        !$AD II-LOOP
        do ii = 0, isize * jsize * ksize - 1
            i = mod(ii, isize) + ibeg
            j = mod(ii / (isize), jsize) + jbeg
            k = ii / ((isize * jsize)) + kbeg
#else
            do k = kbeg, kend
                do j = jbeg, jend
                    do i = ibeg, iend
#endif
                        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)
                        p(i, j, k) = max(p(i, j, k), 1.e-4_realtype * pinfcorr)
 
#ifdef TAPENADE_REVERSE
                    end do
#else
                end do
            end do
        end do
#endif
    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
 
        ! 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)
                        w(i, j, k, irhoe) = max(w(i, j, k, irhoe), &
                                                1.e-5_realtype * pinf)
                    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 <= 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 >= 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 > 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 >= cpeint(nn - 1)) then
 
                                    ! This is the correct range. Exit the do-loop.
 
                                    exit
 
                                end if
 
                                ! Modify ii for the next branch to search.
 
                                ii = ii / 2
 
                            end do interval
 
                            ! 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 = -one                      ! cv/r = cp/r - 1.0
                                e = cptempfit(nn)%eint0 - t  ! e = integral of cv,
                                ! 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) == -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
 
                                ! Exit the Newton loop if the update is smaller
                                ! than the threshold value.
 
                                if (abs(dt) < dtstop) exit
 
                            end do newton
 
                            ! Create the nonDimensional temperature.
 
                            t = t * trefinv
 
                        end if
 
                        ! Compute the pressure from the known temperature
                        ! and density. Include the correction if a k-equation
                        ! is present.
 
                        w(i, j, k, irhoe) = w(ip, jp, kp, irho) * rgas * t
                        w(i, j, k, irhoe) = max(w(i, j, k, irhoe), &
                                                1.e-5_realtype * pinf)
                        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
 
    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, 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) return
 
        ! 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
#ifdef TAPENADE_REVERSE
            isize = (iend - ibeg) + 1
            jsize = (jend - jbeg) + 1
            ksize = (kend - kbeg) + 1
 
            !$AD II-LOOP
            do ii = 0, isize * jsize * ksize - 1
                i = mod(ii, isize) + ibeg
                j = mod(ii / (isize), jsize) + jbeg
                k = ii / ((isize * jsize)) + kbeg
#else
                do k = kbeg, kend
                    do j = jbeg, jend
                        do i = ibeg, iend
#endif
                            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)
#ifdef TAPENADE_REVERSE
                        end do
#else
                    end do
                end do
            end do
#endif
        else
            ! Loop over the owned cells *AND* first level halos of this
            ! block and compute the laminar viscosity ratio.
#ifdef TAPENADE_REVERSE
            isize = (iend - ibeg) + 1
            jsize = (jend - jbeg) + 1
            ksize = (kend - kbeg) + 1
 
            !$AD II-LOOP
            do ii = 0, isize * jsize * ksize - 1
                i = mod(ii, isize) + ibeg
                j = mod(ii / (isize), jsize) + jbeg
                k = ii / ((isize * jsize)) + kbeg
#else
                do k = kbeg, kend
                    do j = jbeg, jend
                        do i = ibeg, iend
#endif
 
                            ! 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)
#ifdef TAPENADE_REVERSE
                        end do
#else
                    end do
                end do
            end do
#endif
        end if
    end subroutine computelamviscosity
 
    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
 
    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, 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
 
        ! 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
 
    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, 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
 
        ! Normalize the input vector.
 
        rnorm = sqrt(refdirection(1)**2 + refdirection(2)**2 + refdirection(3)**2)
        xbn = refdirection(1) / rnorm
        ybn = refdirection(2) / rnorm
        zbn = refdirection(3) / rnorm
 
!!$      ! Compute the wind direction vector.
!!$
!!$      ! 1) rotate alpha radians cw about y-axis
!!$      !    ( <=> rotate y-axis alpha radians ccw)
!!$
!!$      call vectorRotation(x1, y1, z1, 2, alpha, xbn, ybn, zbn)
!!$
!!$      ! 2) rotate beta radians ccw about z-axis
!!$      !    ( <=> rotate z-axis -beta radians ccw)
!!$
!!$      call vectorRotation(xw, yw, zw, 3, -beta, x1, y1, z1)
 
        if (liftindex == 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)
 
        elseif (liftindex == 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
    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
 
        ! rotation about specified axis by specified angle
 
        select case (iaxis)
 
            ! rotation about the x-axis
 
        case (1)
            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
 
    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.
 
#ifdef TAPENADE_REVERSE
        !$AD II-LOOP
        do ii = 0, il * jl * ke - 1
            i = mod(ii, il) + 1
            j = mod(ii / il, jl) + 1
            k = ii / (il * jl) + 1
#else
            do k = 1, ke
                do j = 1, jl
                    do i = 1, il
#endif
                        ! 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 > 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 < 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
#ifdef TAPENADE_REVERSE
                    end do
#else
                end do
            end do
        end do
#endif
 
        ! Second part. Contribution in the j-direction.
        ! The contribution is scattered to both the left and right node
        ! in j-direction.
 
#ifdef TAPENADE_REVERSE
        !$AD II-LOOP
        do ii = 0, il * je * kl - 1
            i = mod(ii, il) + 1
            j = mod(ii / il, je) + 1
            k = ii / (il * je) + 1
#else
            do k = 1, kl
                do j = 1, je
                    do i = 1, il
#endif
 
                        ! 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 > 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 < 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
#ifdef TAPENADE_REVERSE
                    end do
#else
                end do
            end do
        end do
#endif
 
        ! Third part. Contribution in the i-direction.
        ! The contribution is scattered to both the left and right node
        ! in i-direction.
 
#ifdef TAPENADE_REVERSE
        !$AD II-LOOP
        do ii = 0, ie * jl * kl - 1
            i = mod(ii, ie) + 1
            j = mod(ii / ie, jl) + 1
            k = ii / (ie * jl) + 1
#else
            do k = 1, kl
                do j = 1, jl
                    do i = 1, ie
#endif
 
                        ! 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 > 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 < 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
#ifdef TAPENADE_REVERSE
                    end do
#else
                end do
            end do
        end do
#endif
        ! Divide by 8 times the volume to obtain the correct gradients.
 
#ifdef TAPENADE_REVERSE
        !$AD II-LOOP
        do ii = 0, il * jl * kl - 1
            i = mod(ii, il) + 1
            j = mod(ii / il, jl) + 1
            k = ii / (il * jl) + 1
#else
            do k = 1, kl
                do j = 1, jl
                    do i = 1, il
#endif
                        ! 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
#ifdef TAPENADE_REVERSE
                    end do
#else
                end do
            end do
        end do
#endif
    end subroutine allnodalgradients
 
    ! ----------------------------------------------------------------------
    !                                                                      |
    !                    No Tapenade Routine below this line               |
    !                                                                      |
    ! ----------------------------------------------------------------------
 
#ifndef  USE_TAPENADE
 
    subroutine fixallnodalgradientsfromad
 
        use constants
        use blockpointers, only: il, jl, kl, vol, ux, uy, uz, vx, vy, vz, &
                                 wx, wy, wz, qx, qy, qz, vol
        implicit none
 
        ! WorkingVariables
        integer(kind=intType) :: i, j, k
        real(kind=realtype) :: oneoverv
 
        ! So the all nodal gradients doesnt' perform the final
        ! scalaing by the volume since it isn't necessary for the
        ! derivative. We need to fix that here:
 
        ! 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 fixallnodalgradientsfromad
 
    subroutine computeetotcellcpfit(i, j, k, scale, correctForK)
        !
        !       ComputeEtotCellCpfit will compute the total energy for the
        !       given cell of the block given by the current pointers with the
        !       cp temperature curve fit model.
        !
        use constants
        use cpcurvefits
        use blockpointers, only: gamma, w, p
        use flowvarrefstate, only: rgas, tref
        implicit none
        !
        !      Local parameter.
        !
        real(kind=realtype), parameter :: twothird = two * third
        !
        !      Subroutine arguments.
        !
        integer(kind=intType), intent(in) :: i, j, k
        real(kind=realtype), intent(in) :: scale
        logical, intent(in) :: correctForK
        !
        !      Local variables.
        !
        integer(kind=intType) :: nn, mm, ii, start
 
        real(kind=realtype) :: pp, t, t2, cv, eint
 
        ! Compute the dimensional temperature.
 
        pp = p(i, j, k)
        if (correctfork) pp = pp - twothird &
                              * w(i, j, k, irho) * w(i, j, k, itu1)
        t = tref * pp / (rgas * w(i, j, k, irho))
 
        ! Determine the case we are having here.
 
        if (t <= cptrange(0)) then
 
            ! Temperature is less than the smallest temperature
            ! in the curve fits. Use extrapolation using
            ! constant cv.
 
            eint = scale * (cpeint(0) + cv0 * (t - cptrange(0)))
            gamma(i, j, k) = (cv0 + one) / cv0
 
        else if (t >= cptrange(cpnparts)) then
 
            ! Temperature is larger than the largest temperature
            ! in the curve fits. Use extrapolation using
            ! constant cv.
 
            eint = scale * (cpeint(cpnparts) &
                            + cvn * (t - cptrange(cpnparts)))
 
            gamma(i, j, k) = (cvn + one) / cvn
 
        else
 
            ! Temperature is in the curve fit range.
            ! First find the valid range.
 
            ii = cpnparts
            start = 1
            interval: do
 
                ! Next guess for the interval.
 
                nn = start + ii / 2
 
                ! Determine the situation we are having here.
 
                if (t > cptrange(nn)) then
 
                    ! Temperature is larger than the upper boundary of
                    ! the current interval. Update the lower boundary.
 
                    start = nn + 1
                    ii = ii - 1
 
                else if (t >= cptrange(nn - 1)) then
 
                    ! This is the correct range. Exit the do-loop.
 
                    exit
 
                end if
 
                ! Modify ii for the next branch to search.
 
                ii = ii / 2
 
            end do interval
 
            ! Nn contains the correct curve fit interval.
            ! Integrate cv to compute eint.
 
            eint = cptempfit(nn)%eint0 - t
            cv = -one
 
            do ii = 1, cptempfit(nn)%nterm
                t2 = t**cptempfit(nn)%exponents(ii)
                cv = cv + cptempfit(nn)%constants(ii) * t2
 
                if (cptempfit(nn)%exponents(ii) == -1) then
                    eint = eint + cptempfit(nn)%constants(ii) * log(t)
                else
                    mm = cptempfit(nn)%exponents(ii) + 1
                    t2 = t * t2
                    eint = eint + cptempfit(nn)%constants(ii) * t2 / mm
                end if
            end do
 
            eint = scale * eint
            gamma(i, j, k) = (cv + one) / cv
 
        end if
 
        ! Compute the total energy per unit volume.
 
        w(i, j, k, irhoe) = w(i, j, k, irho) * (eint &
                                                + half * (w(i, j, k, ivx)**2 &
                                                          + w(i, j, k, ivy)**2 &
                                                          + w(i, j, k, ivz)**2))
 
        if (correctfork) &
            w(i, j, k, irhoe) = w(i, j, k, irhoe) &
                                + w(i, j, k, irho) * w(i, j, k, itu1)
 
    end subroutine computeetotcellcpfit
 
    subroutine updategamma
        !
        !       This is a utility routine to update the gamma variable from
        !       from gammaConstant if gammaConstant has changed.
        !
        use constants
        use blockpointers, only: ndom, gamma
        use inputtimespectral, only: ntimeintervalsspectral
        use inputphysics, only: gammaconstant
        use utils, only: setpointers
        implicit none
 
        integer :: nn, sps
 
        do sps = 1, ntimeintervalsspectral
            do nn = 1, ndom
                call setpointers(nn, 1_inttype, sps)
                gamma = gammaconstant
            end do
        end do
    end subroutine updategamma
#endif
end module flowutils