turbUtils_b.f90

Go to the documentation of this file.

!        generated by tapenade     (inria, ecuador team)
!  tapenade 3.16 (develop) - 22 aug 2023 15:51
!
module turbutils_b
  implicit none
 
contains
  subroutine prodkatolaunder()
!
!       prodkatolaunder computes the turbulent production term using
!       the kato-launder formulation.
!
    use constants
    use blockpointers, only : nx, ny, nz, il, jl, kl, w, si, sj, sk, &
&   vol, sectionid, scratch
    use flowvarrefstate, only : timeref
    use section, only : sections
    use turbmod, only : prod
    implicit none
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii
    real(kind=realtype) :: uux, uuy, uuz, vvx, vvy, vvz, wwx, wwy, wwz
    real(kind=realtype) :: qxx, qyy, qzz, qxy, qxz, qyz, sijsij
    real(kind=realtype) :: oxy, oxz, oyz, oijoij
    real(kind=realtype) :: fact, omegax, omegay, omegaz
    intrinsic mod
    intrinsic sqrt
! determine the non-dimensional wheel speed of this block.
! the vorticity term, which appears in kato-launder is of course
! not frame invariant. to approximate frame invariance the wheel
! speed should be substracted from oxy, oxz and oyz, which results
! in the vorticity in the rotating frame. however some people
! claim that the absolute vorticity should be used to obtain the
! best results. in that omega should be set to zero.
    omegax = timeref*sections(sectionid)%rotrate(1)
    omegay = timeref*sections(sectionid)%rotrate(2)
    omegaz = timeref*sections(sectionid)%rotrate(3)
!$ad ii-loop
! loop over the cell centers of the given block. it may be more
! efficient to loop over the faces and to scatter the gradient,
! but in that case the gradients for u, v and w must be stored.
! in the current approach no extra memory is needed.
    do ii=0,nx*ny*nz-1
      i = mod(ii, nx) + 2
      j = mod(ii/nx, ny) + 2
      k = ii/(nx*ny) + 2
! compute the gradient of u in the cell center. use is made
! of the fact that the surrounding normals sum up to zero,
! such that the cell i,j,k does not give a contribution.
! the gradient is scaled by a factor 2*vol.
      uux = w(i+1, j, k, ivx)*si(i, j, k, 1) - w(i-1, j, k, ivx)*si(i-1&
&       , j, k, 1) + w(i, j+1, k, ivx)*sj(i, j, k, 1) - w(i, j-1, k, ivx&
&       )*sj(i, j-1, k, 1) + w(i, j, k+1, ivx)*sk(i, j, k, 1) - w(i, j, &
&       k-1, ivx)*sk(i, j, k-1, 1)
      uuy = w(i+1, j, k, ivx)*si(i, j, k, 2) - w(i-1, j, k, ivx)*si(i-1&
&       , j, k, 2) + w(i, j+1, k, ivx)*sj(i, j, k, 2) - w(i, j-1, k, ivx&
&       )*sj(i, j-1, k, 2) + w(i, j, k+1, ivx)*sk(i, j, k, 2) - w(i, j, &
&       k-1, ivx)*sk(i, j, k-1, 2)
      uuz = w(i+1, j, k, ivx)*si(i, j, k, 3) - w(i-1, j, k, ivx)*si(i-1&
&       , j, k, 3) + w(i, j+1, k, ivx)*sj(i, j, k, 3) - w(i, j-1, k, ivx&
&       )*sj(i, j-1, k, 3) + w(i, j, k+1, ivx)*sk(i, j, k, 3) - w(i, j, &
&       k-1, ivx)*sk(i, j, k-1, 3)
! idem for the gradient of v.
      vvx = w(i+1, j, k, ivy)*si(i, j, k, 1) - w(i-1, j, k, ivy)*si(i-1&
&       , j, k, 1) + w(i, j+1, k, ivy)*sj(i, j, k, 1) - w(i, j-1, k, ivy&
&       )*sj(i, j-1, k, 1) + w(i, j, k+1, ivy)*sk(i, j, k, 1) - w(i, j, &
&       k-1, ivy)*sk(i, j, k-1, 1)
      vvy = w(i+1, j, k, ivy)*si(i, j, k, 2) - w(i-1, j, k, ivy)*si(i-1&
&       , j, k, 2) + w(i, j+1, k, ivy)*sj(i, j, k, 2) - w(i, j-1, k, ivy&
&       )*sj(i, j-1, k, 2) + w(i, j, k+1, ivy)*sk(i, j, k, 2) - w(i, j, &
&       k-1, ivy)*sk(i, j, k-1, 2)
      vvz = w(i+1, j, k, ivy)*si(i, j, k, 3) - w(i-1, j, k, ivy)*si(i-1&
&       , j, k, 3) + w(i, j+1, k, ivy)*sj(i, j, k, 3) - w(i, j-1, k, ivy&
&       )*sj(i, j-1, k, 3) + w(i, j, k+1, ivy)*sk(i, j, k, 3) - w(i, j, &
&       k-1, ivy)*sk(i, j, k-1, 3)
! and for the gradient of w.
      wwx = w(i+1, j, k, ivz)*si(i, j, k, 1) - w(i-1, j, k, ivz)*si(i-1&
&       , j, k, 1) + w(i, j+1, k, ivz)*sj(i, j, k, 1) - w(i, j-1, k, ivz&
&       )*sj(i, j-1, k, 1) + w(i, j, k+1, ivz)*sk(i, j, k, 1) - w(i, j, &
&       k-1, ivz)*sk(i, j, k-1, 1)
      wwy = w(i+1, j, k, ivz)*si(i, j, k, 2) - w(i-1, j, k, ivz)*si(i-1&
&       , j, k, 2) + w(i, j+1, k, ivz)*sj(i, j, k, 2) - w(i, j-1, k, ivz&
&       )*sj(i, j-1, k, 2) + w(i, j, k+1, ivz)*sk(i, j, k, 2) - w(i, j, &
&       k-1, ivz)*sk(i, j, k-1, 2)
      wwz = w(i+1, j, k, ivz)*si(i, j, k, 3) - w(i-1, j, k, ivz)*si(i-1&
&       , j, k, 3) + w(i, j+1, k, ivz)*sj(i, j, k, 3) - w(i, j-1, k, ivz&
&       )*sj(i, j-1, k, 3) + w(i, j, k+1, ivz)*sk(i, j, k, 3) - w(i, j, &
&       k-1, ivz)*sk(i, j, k-1, 3)
! compute the strain and vorticity terms. the multiplication
! is present to obtain the correct gradients. note that
! the wheel speed is substracted from the vorticity terms.
      fact = half/vol(i, j, k)
      qxx = fact*uux
      qyy = fact*vvy
      qzz = fact*wwz
      qxy = fact*half*(uuy+vvx)
      qxz = fact*half*(uuz+wwx)
      qyz = fact*half*(vvz+wwy)
      oxy = fact*half*(vvx-uuy) - omegaz
      oxz = fact*half*(uuz-wwx) - omegay
      oyz = fact*half*(wwy-vvz) - omegax
! compute the summation of the strain and vorticity tensors.
      sijsij = two*(qxy**2+qxz**2+qyz**2) + qxx**2 + qyy**2 + qzz**2
      oijoij = two*(oxy**2+oxz**2+oyz**2)
! compute the production term.
      scratch(i, j, k, iprod) = two*sqrt(sijsij*oijoij)
    end do
  end subroutine prodkatolaunder
 
  subroutine prodsmag2()
!
!       prodsmag2 computes the term:
!              2*sij*sij - 2/3 div(u)**2 with  sij=0.5*(duidxj+dujdxi)
!       which is used for the turbulence equations.
!       it is assumed that the pointer prod, stored in turbmod, is
!       already set to the correct entry.
!
    use constants
    use blockpointers, only : nx, ny, nz, il, jl, kl, w, si, sj, sk, &
&   vol, sectionid, scratch
    implicit none
!
!      local parameter
!
    real(kind=realtype), parameter :: f23=two*third
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii
    real(kind=realtype) :: uux, uuy, uuz, vvx, vvy, vvz, wwx, wwy, wwz
    real(kind=realtype) :: div2, fact, sxx, syy, szz, sxy, sxz, syz
    intrinsic mod
!$ad ii-loop
! loop over the cell centers of the given block. it may be more
! efficient to loop over the faces and to scatter the gradient,
! but in that case the gradients for u, v and w must be stored.
! in the current approach no extra memory is needed.
    do ii=0,nx*ny*nz-1
      i = mod(ii, nx) + 2
      j = mod(ii/nx, ny) + 2
      k = ii/(nx*ny) + 2
! compute the gradient of u in the cell center. use is made
! of the fact that the surrounding normals sum up to zero,
! such that the cell i,j,k does not give a contribution.
! the gradient is scaled by the factor 2*vol.
      uux = w(i+1, j, k, ivx)*si(i, j, k, 1) - w(i-1, j, k, ivx)*si(i-1&
&       , j, k, 1) + w(i, j+1, k, ivx)*sj(i, j, k, 1) - w(i, j-1, k, ivx&
&       )*sj(i, j-1, k, 1) + w(i, j, k+1, ivx)*sk(i, j, k, 1) - w(i, j, &
&       k-1, ivx)*sk(i, j, k-1, 1)
      uuy = w(i+1, j, k, ivx)*si(i, j, k, 2) - w(i-1, j, k, ivx)*si(i-1&
&       , j, k, 2) + w(i, j+1, k, ivx)*sj(i, j, k, 2) - w(i, j-1, k, ivx&
&       )*sj(i, j-1, k, 2) + w(i, j, k+1, ivx)*sk(i, j, k, 2) - w(i, j, &
&       k-1, ivx)*sk(i, j, k-1, 2)
      uuz = w(i+1, j, k, ivx)*si(i, j, k, 3) - w(i-1, j, k, ivx)*si(i-1&
&       , j, k, 3) + w(i, j+1, k, ivx)*sj(i, j, k, 3) - w(i, j-1, k, ivx&
&       )*sj(i, j-1, k, 3) + w(i, j, k+1, ivx)*sk(i, j, k, 3) - w(i, j, &
&       k-1, ivx)*sk(i, j, k-1, 3)
! idem for the gradient of v.
      vvx = w(i+1, j, k, ivy)*si(i, j, k, 1) - w(i-1, j, k, ivy)*si(i-1&
&       , j, k, 1) + w(i, j+1, k, ivy)*sj(i, j, k, 1) - w(i, j-1, k, ivy&
&       )*sj(i, j-1, k, 1) + w(i, j, k+1, ivy)*sk(i, j, k, 1) - w(i, j, &
&       k-1, ivy)*sk(i, j, k-1, 1)
      vvy = w(i+1, j, k, ivy)*si(i, j, k, 2) - w(i-1, j, k, ivy)*si(i-1&
&       , j, k, 2) + w(i, j+1, k, ivy)*sj(i, j, k, 2) - w(i, j-1, k, ivy&
&       )*sj(i, j-1, k, 2) + w(i, j, k+1, ivy)*sk(i, j, k, 2) - w(i, j, &
&       k-1, ivy)*sk(i, j, k-1, 2)
      vvz = w(i+1, j, k, ivy)*si(i, j, k, 3) - w(i-1, j, k, ivy)*si(i-1&
&       , j, k, 3) + w(i, j+1, k, ivy)*sj(i, j, k, 3) - w(i, j-1, k, ivy&
&       )*sj(i, j-1, k, 3) + w(i, j, k+1, ivy)*sk(i, j, k, 3) - w(i, j, &
&       k-1, ivy)*sk(i, j, k-1, 3)
! and for the gradient of w.
      wwx = w(i+1, j, k, ivz)*si(i, j, k, 1) - w(i-1, j, k, ivz)*si(i-1&
&       , j, k, 1) + w(i, j+1, k, ivz)*sj(i, j, k, 1) - w(i, j-1, k, ivz&
&       )*sj(i, j-1, k, 1) + w(i, j, k+1, ivz)*sk(i, j, k, 1) - w(i, j, &
&       k-1, ivz)*sk(i, j, k-1, 1)
      wwy = w(i+1, j, k, ivz)*si(i, j, k, 2) - w(i-1, j, k, ivz)*si(i-1&
&       , j, k, 2) + w(i, j+1, k, ivz)*sj(i, j, k, 2) - w(i, j-1, k, ivz&
&       )*sj(i, j-1, k, 2) + w(i, j, k+1, ivz)*sk(i, j, k, 2) - w(i, j, &
&       k-1, ivz)*sk(i, j, k-1, 2)
      wwz = w(i+1, j, k, ivz)*si(i, j, k, 3) - w(i-1, j, k, ivz)*si(i-1&
&       , j, k, 3) + w(i, j+1, k, ivz)*sj(i, j, k, 3) - w(i, j-1, k, ivz&
&       )*sj(i, j-1, k, 3) + w(i, j, k+1, ivz)*sk(i, j, k, 3) - w(i, j, &
&       k-1, ivz)*sk(i, j, k-1, 3)
! compute the components of the stress tensor.
! the combination of the current scaling of the velocity
! gradients (2*vol) and the definition of the stress tensor,
! leads to the factor 1/(4*vol).
      fact = fourth/vol(i, j, k)
      sxx = two*fact*uux
      syy = two*fact*vvy
      szz = two*fact*wwz
      sxy = fact*(uuy+vvx)
      sxz = fact*(uuz+wwx)
      syz = fact*(vvz+wwy)
! compute 2/3 * divergence of velocity squared
      div2 = f23*(sxx+syy+szz)**2
! store the square of strain as the production term.
      scratch(i, j, k, iprod) = two*(two*(sxy**2+sxz**2+syz**2)+sxx**2+&
&       syy**2+szz**2) - div2
    end do
  end subroutine prodsmag2
 
  subroutine prodwmag2()
!
!       prodwmag2 computes the term:
!          2*oij*oij  with oij=0.5*(duidxj - dujdxi).
!       this is equal to the magnitude squared of the vorticity.
!       it is assumed that the pointer vort, stored in turbmod, is
!       already set to the correct entry.
!
    use constants
    use blockpointers, only : nx, ny, nz, il, jl, kl, w, si, sj, sk, &
&   vol, sectionid, scratch
    use flowvarrefstate, only : timeref
    use section, only : sections
    implicit none
!
!      local variables.
!
    integer :: i, j, k, ii
    real(kind=realtype) :: uuy, uuz, vvx, vvz, wwx, wwy
    real(kind=realtype) :: fact, vortx, vorty, vortz
    real(kind=realtype) :: omegax, omegay, omegaz
    intrinsic mod
! determine the non-dimensional wheel speed of this block.
    omegax = timeref*sections(sectionid)%rotrate(1)
    omegay = timeref*sections(sectionid)%rotrate(2)
    omegaz = timeref*sections(sectionid)%rotrate(3)
!$ad ii-loop
! loop over the cell centers of the given block. it may be more
! efficient to loop over the faces and to scatter the gradient,
! but in that case the gradients for u, v and w must be stored.
! in the current approach no extra memory is needed.
    do ii=0,nx*ny*nz-1
      i = mod(ii, nx) + 2
      j = mod(ii/nx, ny) + 2
      k = ii/(nx*ny) + 2
! compute the necessary derivatives of u in the cell center.
! use is made of the fact that the surrounding normals sum up
! to zero, such that the cell i,j,k does not give a
! contribution. the gradient is scaled by a factor 2*vol.
      uuy = w(i+1, j, k, ivx)*si(i, j, k, 2) - w(i-1, j, k, ivx)*si(i-1&
&       , j, k, 2) + w(i, j+1, k, ivx)*sj(i, j, k, 2) - w(i, j-1, k, ivx&
&       )*sj(i, j-1, k, 2) + w(i, j, k+1, ivx)*sk(i, j, k, 2) - w(i, j, &
&       k-1, ivx)*sk(i, j, k-1, 2)
      uuz = w(i+1, j, k, ivx)*si(i, j, k, 3) - w(i-1, j, k, ivx)*si(i-1&
&       , j, k, 3) + w(i, j+1, k, ivx)*sj(i, j, k, 3) - w(i, j-1, k, ivx&
&       )*sj(i, j-1, k, 3) + w(i, j, k+1, ivx)*sk(i, j, k, 3) - w(i, j, &
&       k-1, ivx)*sk(i, j, k-1, 3)
! idem for the gradient of v.
      vvx = w(i+1, j, k, ivy)*si(i, j, k, 1) - w(i-1, j, k, ivy)*si(i-1&
&       , j, k, 1) + w(i, j+1, k, ivy)*sj(i, j, k, 1) - w(i, j-1, k, ivy&
&       )*sj(i, j-1, k, 1) + w(i, j, k+1, ivy)*sk(i, j, k, 1) - w(i, j, &
&       k-1, ivy)*sk(i, j, k-1, 1)
      vvz = w(i+1, j, k, ivy)*si(i, j, k, 3) - w(i-1, j, k, ivy)*si(i-1&
&       , j, k, 3) + w(i, j+1, k, ivy)*sj(i, j, k, 3) - w(i, j-1, k, ivy&
&       )*sj(i, j-1, k, 3) + w(i, j, k+1, ivy)*sk(i, j, k, 3) - w(i, j, &
&       k-1, ivy)*sk(i, j, k-1, 3)
! and for the gradient of w.
      wwx = w(i+1, j, k, ivz)*si(i, j, k, 1) - w(i-1, j, k, ivz)*si(i-1&
&       , j, k, 1) + w(i, j+1, k, ivz)*sj(i, j, k, 1) - w(i, j-1, k, ivz&
&       )*sj(i, j-1, k, 1) + w(i, j, k+1, ivz)*sk(i, j, k, 1) - w(i, j, &
&       k-1, ivz)*sk(i, j, k-1, 1)
      wwy = w(i+1, j, k, ivz)*si(i, j, k, 2) - w(i-1, j, k, ivz)*si(i-1&
&       , j, k, 2) + w(i, j+1, k, ivz)*sj(i, j, k, 2) - w(i, j-1, k, ivz&
&       )*sj(i, j-1, k, 2) + w(i, j, k+1, ivz)*sk(i, j, k, 2) - w(i, j, &
&       k-1, ivz)*sk(i, j, k-1, 2)
! compute the three components of the vorticity vector.
! substract the part coming from the rotating frame.
      fact = half/vol(i, j, k)
      vortx = fact*(wwy-vvz) - two*omegax
      vorty = fact*(uuz-wwx) - two*omegay
      vortz = fact*(vvx-uuy) - two*omegaz
! compute the magnitude squared of the vorticity.
      scratch(i, j, k, ivort) = vortx**2 + vorty**2 + vortz**2
    end do
  end subroutine prodwmag2
 
!  differentiation of sanuknowneddyratio in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: sanuknowneddyratio
!   with respect to varying inputs: nulam
  subroutine sanuknowneddyratio_b(eddyratio, nulam, nulamd, &
&   sanuknowneddyratiod)
!
!       sanuknowneddyratio computes the spalart-allmaras transport
!       variable nu for the given eddy viscosity ratio.
!
    use constants
    use paramturb
    implicit none
!
!      function type.
!
    real(kind=realtype) :: sanuknowneddyratio
    real(kind=realtype) :: sanuknowneddyratiod
!
!      function arguments.
!
    real(kind=realtype), intent(in) :: eddyratio, nulam
    real(kind=realtype) :: nulamd
!
!      local variables.
!
    real(kind=realtype) :: cv13, chi, chi2, chi3, chi4, f, df, dchi
    intrinsic abs
    real(kind=realtype) :: abs0
    integer :: ad_count
    integer :: i
    integer :: branch
! take care of the exceptional cases.
    if (eddyratio .le. zero) then
      nulamd = 0.0_8
    else
! set the value of cv1^3, which is the constant appearing in the
! sa function fv1 to compute the eddy viscosity
      cv13 = rsacv1**3
! determine the value of chi, which is given by the quartic
! polynomial chi^4 - ratio*(chi^3 + cv1^3) = 0.
! first determine the start value, depending on the eddyratio.
      if (eddyratio .lt. 1.e-4_realtype) then
        call pushcontrol2b(2)
        chi = 0.5_realtype
      else if (eddyratio .lt. 1.0_realtype) then
        call pushcontrol2b(1)
        chi = 5.0_realtype
      else if (eddyratio .lt. 10.0_realtype) then
        call pushcontrol2b(0)
        chi = 10.0_realtype
      else
        call pushcontrol2b(0)
        chi = eddyratio
      end if
      ad_count = 1
! the actual newton algorithm.
      do 
! compute the function value and the derivative.
        chi2 = chi*chi
        chi3 = chi*chi2
        chi4 = chi*chi3
        f = chi4 - eddyratio*(chi3+cv13)
        df = four*chi3 - three*eddyratio*chi2
! compute the negative update and the new value of chi.
        dchi = f/df
        chi = chi - dchi
        if (dchi/chi .ge. 0.) then
          abs0 = dchi/chi
        else
          abs0 = -(dchi/chi)
        end if
! condition to exit the loop.
        if (abs0 .le. thresholdreal) then
          exit
        else
          ad_count = ad_count + 1
        end if
      end do
      call pushinteger4(ad_count)
      nulamd = chi*sanuknowneddyratiod
      call popinteger4(ad_count)
      i = ad_count + 1
      call popcontrol2b(branch)
    end if
  end subroutine sanuknowneddyratio_b
 
  function sanuknowneddyratio(eddyratio, nulam)
!
!       sanuknowneddyratio computes the spalart-allmaras transport
!       variable nu for the given eddy viscosity ratio.
!
    use constants
    use paramturb
    implicit none
!
!      function type.
!
    real(kind=realtype) :: sanuknowneddyratio
!
!      function arguments.
!
    real(kind=realtype), intent(in) :: eddyratio, nulam
!
!      local variables.
!
    real(kind=realtype) :: cv13, chi, chi2, chi3, chi4, f, df, dchi
    intrinsic abs
    real(kind=realtype) :: abs0
! take care of the exceptional cases.
    if (eddyratio .le. zero) then
      sanuknowneddyratio = zero
      return
    else
! set the value of cv1^3, which is the constant appearing in the
! sa function fv1 to compute the eddy viscosity
      cv13 = rsacv1**3
! determine the value of chi, which is given by the quartic
! polynomial chi^4 - ratio*(chi^3 + cv1^3) = 0.
! first determine the start value, depending on the eddyratio.
      if (eddyratio .lt. 1.e-4_realtype) then
        chi = 0.5_realtype
      else if (eddyratio .lt. 1.0_realtype) then
        chi = 5.0_realtype
      else if (eddyratio .lt. 10.0_realtype) then
        chi = 10.0_realtype
      else
        chi = eddyratio
      end if
! the actual newton algorithm.
      do 
! compute the function value and the derivative.
        chi2 = chi*chi
        chi3 = chi*chi2
        chi4 = chi*chi3
        f = chi4 - eddyratio*(chi3+cv13)
        df = four*chi3 - three*eddyratio*chi2
! compute the negative update and the new value of chi.
        dchi = f/df
        chi = chi - dchi
        if (dchi/chi .ge. 0.) then
          abs0 = dchi/chi
        else
          abs0 = -(dchi/chi)
        end if
! condition to exit the loop.
        if (abs0 .le. thresholdreal) then
! chi is the ratio of the spalart allmaras transport variable and
! the laminar viscosity. so multiply chi with the laminar viscosity
! to obtain the correct value.
          sanuknowneddyratio = nulam*chi
          goto 100
        end if
      end do
    end if
 100 continue
  end function sanuknowneddyratio
 
  subroutine unsteadyturbterm(madv, nadv, offset, qq)
!
!       unsteadyturbterm discretizes the time derivative of the
!       turbulence transport equations and add it to the residual.
!       as the time derivative is the same for all turbulence models,
!       this generic routine can be used; both the discretization of
!       the time derivative and its contribution to the central
!       jacobian are computed by this routine.
!       only nadv equations are treated, while the actual system has
!       size madv. the reason is that some equations for some
!       turbulence equations do not have a time derivative, e.g. the
!       f equation in the v2-f model. the argument offset indicates
!       the offset in the w vector where this subsystem starts. as a
!       consequence it is assumed that the indices of the current
!       subsystem are contiguous, e.g. if a 2*2 system is solved the
!       last index in w is offset+1 and offset+2 respectively.
!
    use blockpointers
    use flowvarrefstate
    use inputphysics
    use inputtimespectral
    use inputunsteady
    use iteration
    use section
    use turbmod
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: madv, nadv, offset
    real(kind=realtype), dimension(2:il, 2:jl, 2:kl, madv, madv), &
&   intent(inout) :: qq
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii, jj, nn
    real(kind=realtype) :: oneoverdt, tmp
! determine the equation mode.
    select case  (equationmode) 
    case (steady) 
! steady computation. no time derivative present.
      return
!===============================================================
    case (unsteady) 
! the time deritvative term depends on the integration
! scheme used.
      select case  (timeintegrationscheme) 
      case (bdf) 
! backward difference formula is used as time
! integration scheme.
! store the inverse of the physical nondimensional
! time step a bit easier.
        oneoverdt = timeref/deltat
! loop over the number of turbulent transport equations.
nadvloopunsteady:do ii=1,nadv
! store the index of the current turbulent variable in jj.
          jj = ii + offset
! loop over the owned cells of this block to compute the
! time derivative.
          do k=2,kl
            do j=2,jl
              do i=2,il
! initialize tmp to the value of the current
! level multiplied by the corresponding coefficient
! in the time integration scheme.
                tmp = coeftime(0)*w(i, j, k, jj)
! loop over the old time levels and add the
! corresponding contribution to tmp.
                do nn=1,noldlevels
                  tmp = tmp + coeftime(nn)*wold(nn, i, j, k, jj)
                end do
! update the residual. note that in the turbulent
! routines the residual is defined with an opposite
! sign compared to the residual of the flow equations.
! therefore the time derivative must be substracted
! from dvt.
                scratch(i, j, k, idvt+ii-1) = scratch(i, j, k, idvt+ii-1&
&                 ) - oneoverdt*tmp
! update the central jacobian.
                qq(i, j, k, ii, ii) = qq(i, j, k, ii, ii) + coeftime(0)*&
&                 oneoverdt
              end do
            end do
          end do
        end do nadvloopunsteady
      case (explicitrk) 
!===========================================================
! explicit time integration scheme. the time derivative
! is handled differently.
        return
      end select
!===============================================================
    case (timespectral) 
! time spectral method.
! loop over the number of turbulent transport equations.
nadvloopspectral:do ii=1,nadv
! store the index of the current turbulent variable in jj.
        jj = ii + offset
! the time derivative has been computed earlier in
! unsteadyturbspectral and stored in entry jj of scratch.
! substract this value for all owned cells. it must be
! substracted, because in the turbulent routines the
! residual is defined with an opposite sign compared to
! the residual of the flow equations.
! also add a term to the diagonal matrix, which corresponds
! to to the contribution of the highest frequency. this is
! equivalent to an explicit treatment of the time derivative
! and may need to be changed.
        tmp = ntimeintervalsspectral*pi*timeref/sections(sectionid)%&
&         timeperiod
        do k=2,kl
          do j=2,jl
            do i=2,il
              scratch(i, j, k, idvt+ii-1) = scratch(i, j, k, idvt+ii-1) &
&               - dw(i, j, k, jj)
              qq(i, j, k, ii, ii) = qq(i, j, k, ii, ii) + tmp
            end do
          end do
        end do
      end do nadvloopspectral
    end select
  end subroutine unsteadyturbterm
 
!  differentiation of computeeddyviscosity in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *rev *w *rlv
!   with respect to varying inputs: *rev *w *rlv
!   rw status of diff variables: *rev:in-out *w:incr *rlv:incr
!   plus diff mem management of: rev:in w:in rlv:in
  subroutine computeeddyviscosity_b(includehalos)
!
!       computeeddyviscosity computes the eddy viscosity in the
!       owned cell centers of the given block. it is assumed that the
!       pointes already point to the correct block before entering
!       this subroutine.
!
    use constants
    use flowvarrefstate
    use inputphysics
    use iteration
    use blockpointers
    implicit none
! input parameter
    logical, intent(in) :: includehalos
!
!      local variables.
!
    logical :: returnimmediately
    integer(kind=inttype) :: ibeg, iend, jbeg, jend, kbeg, kend
    integer :: branch
! check if an immediate return can be made.
    if (eddymodel) then
      if (currentlevel .le. groundlevel) then
        call pushcontrol1b(0)
        returnimmediately = .false.
      else
        call pushcontrol1b(0)
        returnimmediately = .true.
      end if
    else
      call pushcontrol1b(1)
      returnimmediately = .true.
    end if
    if (.not.returnimmediately) then
! determine the turbulence model and call the appropriate
! routine to compute the eddy viscosity.
      if (includehalos) then
        ibeg = 1
        iend = ie
        jbeg = 1
        jend = je
        kbeg = 1
        kend = ke
      else
        ibeg = 2
        iend = il
        jbeg = 2
        jend = jl
        kbeg = 2
        kend = kl
      end if
      select case  (turbmodel) 
      case (spalartallmaras, spalartallmarasedwards) 
        call saeddyviscosity_b(ibeg, iend, jbeg, jend, kbeg, kend)
      end select
    end if
    call popcontrol1b(branch)
  end subroutine computeeddyviscosity_b
 
  subroutine computeeddyviscosity(includehalos)
!
!       computeeddyviscosity computes the eddy viscosity in the
!       owned cell centers of the given block. it is assumed that the
!       pointes already point to the correct block before entering
!       this subroutine.
!
    use constants
    use flowvarrefstate
    use inputphysics
    use iteration
    use blockpointers
    implicit none
! input parameter
    logical, intent(in) :: includehalos
!
!      local variables.
!
    logical :: returnimmediately
    integer(kind=inttype) :: ibeg, iend, jbeg, jend, kbeg, kend
! check if an immediate return can be made.
    if (eddymodel) then
      if (currentlevel .le. groundlevel) then
        returnimmediately = .false.
      else
        returnimmediately = .true.
      end if
    else
      returnimmediately = .true.
    end if
    if (returnimmediately) then
      return
    else
! determine the turbulence model and call the appropriate
! routine to compute the eddy viscosity.
      if (includehalos) then
        ibeg = 1
        iend = ie
        jbeg = 1
        jend = je
        kbeg = 1
        kend = ke
      else
        ibeg = 2
        iend = il
        jbeg = 2
        jend = jl
        kbeg = 2
        kend = kl
      end if
      select case  (turbmodel) 
      case (spalartallmaras, spalartallmarasedwards) 
        call saeddyviscosity(ibeg, iend, jbeg, jend, kbeg, kend)
      end select
    end if
  end subroutine computeeddyviscosity
 
!  differentiation of saeddyviscosity in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *rev *w *rlv
!   with respect to varying inputs: *rev *w *rlv
!   plus diff mem management of: rev:in w:in rlv:in
  subroutine saeddyviscosity_b(ibeg, iend, jbeg, jend, kbeg, kend)
!
!       saeddyviscosity computes the eddy-viscosity according to the
!       spalart-allmaras model for the block given in blockpointers.
!       this routine for both the original version as well as the
!       modified version according to edwards.
!
    use constants
    use blockpointers
    use constants
    use paramturb
    implicit none
! input variables
    integer(kind=inttype) :: ibeg, iend, jbeg, jend, kbeg, kend
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii, isize, jsize, ksize
    real(kind=realtype) :: chi, chi3, fv1, rnusa, cv13
    real(kind=realtype) :: chid, chi3d, fv1d, rnusad
    intrinsic mod
    real(kind=realtype) :: tempd
! store the cv1^3; cv1 is a constant of the spalart-allmaras model.
    cv13 = rsacv1**3
! loop over the cells of this block and compute the eddy viscosity.
! do not include halo's.
    isize = iend - ibeg + 1
    jsize = jend - jbeg + 1
    ksize = kend - kbeg + 1
!$bwd-of ii-loop 
    do ii=0,isize*jsize*ksize-1
      i = mod(ii, isize) + ibeg
      j = mod(ii/isize, jsize) + jbeg
      k = ii/(isize*jsize) + kbeg
      rnusa = w(i, j, k, itu1)*w(i, j, k, irho)
      chi = rnusa/rlv(i, j, k)
      chi3 = chi**3
      fv1 = chi3/(chi3+cv13)
      fv1d = rnusa*revd(i, j, k)
      tempd = fv1d/(cv13+chi3)
      chi3d = (1.0-chi3/(cv13+chi3))*tempd
      chid = 3*chi**2*chi3d
      tempd = chid/rlv(i, j, k)
      rnusad = fv1*revd(i, j, k) + tempd
      revd(i, j, k) = 0.0_8
      rlvd(i, j, k) = rlvd(i, j, k) - rnusa*tempd/rlv(i, j, k)
      wd(i, j, k, itu1) = wd(i, j, k, itu1) + w(i, j, k, irho)*rnusad
      wd(i, j, k, irho) = wd(i, j, k, irho) + w(i, j, k, itu1)*rnusad
    end do
  end subroutine saeddyviscosity_b
 
  subroutine saeddyviscosity(ibeg, iend, jbeg, jend, kbeg, kend)
!
!       saeddyviscosity computes the eddy-viscosity according to the
!       spalart-allmaras model for the block given in blockpointers.
!       this routine for both the original version as well as the
!       modified version according to edwards.
!
    use constants
    use blockpointers
    use constants
    use paramturb
    implicit none
! input variables
    integer(kind=inttype) :: ibeg, iend, jbeg, jend, kbeg, kend
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii, isize, jsize, ksize
    real(kind=realtype) :: chi, chi3, fv1, rnusa, cv13
    intrinsic mod
! store the cv1^3; cv1 is a constant of the spalart-allmaras model.
    cv13 = rsacv1**3
! loop over the cells of this block and compute the eddy viscosity.
! do not include halo's.
    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
      rnusa = w(i, j, k, itu1)*w(i, j, k, irho)
      chi = rnusa/rlv(i, j, k)
      chi3 = chi**3
      fv1 = chi3/(chi3+cv13)
      rev(i, j, k) = fv1*rnusa
    end do
  end subroutine saeddyviscosity
 
  subroutine kweddyviscosity(ibeg, iend, jbeg, jend, kbeg, kend)
!
!       kweddyviscosity computes the eddy viscosity according to the
!       k-omega models (both the original wilcox as well as the
!       modified version) for the block given in blockpointers.
!
    use constants
    use blockpointers
    implicit none
! input variables
    integer(kind=inttype) :: ibeg, iend, jbeg, jend, kbeg, kend
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii, isize, jsize, ksize
    intrinsic mod
    intrinsic abs
    real(kind=realtype) :: x1
! loop over the cells of this block and compute the eddy viscosity.
! do not include halo's.
    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
      x1 = w(i, j, k, irho)*w(i, j, k, itu1)/w(i, j, k, itu2)
      if (x1 .ge. 0.) then
        rev(i, j, k) = x1
      else
        rev(i, j, k) = -x1
      end if
    end do
  end subroutine kweddyviscosity
 
  subroutine ssteddyviscosity(ibeg, iend, jbeg, jend, kbeg, kend)
!
!       ssteddyviscosity computes the eddy viscosity according to
!       menter's sst variant of the k-omega turbulence model for the
!       block given in blockpointers.
!
    use constants
    use blockpointers
    use paramturb
    use turbmod
    implicit none
! input variables
    integer(kind=inttype) :: ibeg, iend, jbeg, jend, kbeg, kend
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii, isize, jsize, ksize
    real(kind=realtype) :: t1, t2, arg2, f2, vortmag
    intrinsic mod
    intrinsic sqrt
    intrinsic max
    intrinsic tanh
    real(kind=realtype) :: max1
! compute the vorticity squared in the cell centers. the reason
! for computing the vorticity squared is that a routine exists
! for it; for the actual eddy viscosity computation the vorticity
! itself is needed.
    call prodwmag2()
! loop over the cells of this block and compute the eddy viscosity.
! do not include halo's.
    isize = iend - ibeg + 1
    jsize = jend - jbeg + 1
    ksize = kend - kbeg + 1
!$ad ii-loop
    do ii=0,isize*jsize*ksize-1
      i = mod(ii, isize) + ibeg
      j = mod(ii/isize, jsize) + jbeg
      k = ii/(isize*jsize) + kbeg
! compute the value of the function f2, which occurs in the
! eddy-viscosity computation.
      t1 = two*sqrt(w(i, j, k, itu1))/(0.09_realtype*w(i, j, k, itu2)*&
&       d2wall(i, j, k))
      t2 = 500.0_realtype*rlv(i, j, k)/(w(i, j, k, irho)*w(i, j, k, itu2&
&       )*d2wall(i, j, k)**2)
      if (t1 .lt. t2) then
        arg2 = t2
      else
        arg2 = t1
      end if
      f2 = tanh(arg2**2)
! and compute the eddy viscosity.
      vortmag = sqrt(scratch(i, j, k, iprod))
      if (rssta1*w(i, j, k, itu2) .lt. f2*vortmag) then
        max1 = f2*vortmag
      else
        max1 = rssta1*w(i, j, k, itu2)
      end if
      rev(i, j, k) = w(i, j, k, irho)*rssta1*w(i, j, k, itu1)/max1
    end do
  end subroutine ssteddyviscosity
 
!  differentiation of turbadvection in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *sfacei *sfacej *sfacek *w
!                *scratch *vol *si *sj *sk
!   with respect to varying inputs: *sfacei *sfacej *sfacek *w
!                *scratch *vol *si *sj *sk
!   rw status of diff variables: *sfacei:incr *sfacej:incr *sfacek:incr
!                *w:incr *scratch:in-out *vol:incr *si:incr *sj:incr
!                *sk:incr
!   plus diff mem management of: sfacei:in sfacej:in sfacek:in
!                w:in scratch:in vol:in si:in sj:in sk:in
  subroutine turbadvection_b(madv, nadv, offset, qq)
!
!       turbadvection discretizes the advection part of the turbulent
!       transport equations. as the advection part is the same for all
!       models, this generic routine can be used. both the
!       discretization and the central jacobian are computed in this
!       subroutine. the former can either be 1st or 2nd order
!       accurate; the latter is always based on the 1st order upwind
!       discretization. when the discretization must be second order
!       accurate, the fully upwind (kappa = -1) scheme in combination
!       with the minmod limiter is used.
!       only nadv equations are treated, while the actual system has
!       size madv. the reason is that some equations for some
!       turbulence equations do not have an advection part, e.g. the
!       f equation in the v2-f model. the argument offset indicates
!       the offset in the w vector where this subsystem starts. as a
!       consequence it is assumed that the indices of the current
!       subsystem are contiguous, e.g. if a 2*2 system is solved the
!       last index in w is offset+1 and offset+2 respectively.
!
    use constants
    use blockpointers, only : nx, ny, nz, il, jl, kl, vol, vold, &
&   sfacei, sfaceid, sfacej, sfacejd, sfacek, sfacekd, w, wd, si, sid, &
&   sj, sjd, sk, skd, addgridvelocities, bmti1, bmti2, bmtj1, bmtj2, &
&   bmtk1, bmtk2, scratch, scratchd
    use inputdiscretization, only : orderturb
    use iteration, only : groundlevel
    use turbmod, only : secondord
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: nadv, madv, offset
    real(kind=realtype), dimension(2:il, 2:jl, 2:kl, madv, madv), &
&   intent(inout) :: qq
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii, jj, kk, iii
    real(kind=realtype) :: qs, voli, xa, ya, za
    real(kind=realtype) :: qsd, volid, xad, yad, zad
    real(kind=realtype) :: uu, dwt, dwtm1, dwtp1, dwti, dwtj, dwtk
    real(kind=realtype) :: uud, dwtd, dwtm1d, dwtp1d, dwtid, dwtjd, &
&   dwtkd
    real(kind=realtype), dimension(madv) :: impl
    intrinsic mod
    intrinsic abs
    real(kind=realtype) :: abs0
    real(kind=realtype) :: abs1
    real(kind=realtype) :: abs2
    real(kind=realtype) :: abs3
    real(kind=realtype) :: abs4
    real(kind=realtype) :: abs5
    real(kind=realtype) :: abs6
    real(kind=realtype) :: abs7
    real(kind=realtype) :: abs8
    real(kind=realtype) :: abs9
    real(kind=realtype) :: abs10
    real(kind=realtype) :: abs11
    real(kind=realtype) :: abs12
    real(kind=realtype) :: abs13
    real(kind=realtype) :: abs14
    real(kind=realtype) :: abs15
    real(kind=realtype) :: abs16
    real(kind=realtype) :: abs17
    real(kind=realtype) :: abs18
    real(kind=realtype) :: abs19
    real(kind=realtype) :: abs20
    real(kind=realtype) :: abs21
    real(kind=realtype) :: abs22
    real(kind=realtype) :: abs23
    integer :: branch
! determine whether or not a second order discretization for the
! advective terms must be used.
    secondord = .false.
    if (groundlevel .eq. 1_inttype .and. orderturb .eq. secondorder) &
&     secondord = .true.
    qsd = 0.0_8
    qs = zero
    qsd = 0.0_8
!$bwd-of ii-loop 
    do iii=0,nx*ny*nz-1
      i = mod(iii, nx) + 2
      j = mod(iii/nx, ny) + 2
      k = iii/(nx*ny) + 2
! compute the grid velocity if present.
! it is taken as the average of i and i-1,
      voli = half/vol(i, j, k)
      if (addgridvelocities) then
        qs = (sfacei(i, j, k)+sfacei(i-1, j, k))*voli
        call pushcontrol1b(0)
      else
        call pushcontrol1b(1)
      end if
! compute the normal velocity, where the normal direction
! is taken as the average of faces i and i-1.
      xa = (si(i, j, k, 1)+si(i-1, j, k, 1))*voli
      ya = (si(i, j, k, 2)+si(i-1, j, k, 2))*voli
      za = (si(i, j, k, 3)+si(i-1, j, k, 3))*voli
      uu = xa*w(i, j, k, ivx) + ya*w(i, j, k, ivy) + za*w(i, j, k, ivz) &
&       - qs
! determine the situation we are having here, i.e. positive
! or negative normal velocity.
      if (uu .gt. zero) then
        uud = 0.0_8
!$bwd-of ii-loop 
        do 100 ii=1,nadv
! set the value of jj such that it corresponds to the
! turbulent entry in w.
          jj = ii + offset
! check whether a first or a second order discretization
! must be used.
          if (secondord) then
! second order; store the three differences for the
! discretization of the derivative in i-direction.
            dwtm1 = w(i-1, j, k, jj) - w(i-2, j, k, jj)
            dwt = w(i, j, k, jj) - w(i-1, j, k, jj)
            dwtp1 = w(i+1, j, k, jj) - w(i, j, k, jj)
! construct the derivative in this cell center. this is
! the first order upwind derivative with two nonlinear
! corrections.
            dwti = dwt
            if (dwt*dwtp1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs8 = dwt
              else
                abs8 = -dwt
              end if
              if (dwtp1 .ge. 0.) then
                abs20 = dwtp1
              else
                abs20 = -dwtp1
              end if
              if (abs8 .lt. abs20) then
                dwti = dwti + half*dwt
                call pushcontrol2b(0)
              else
                dwti = dwti + half*dwtp1
                call pushcontrol2b(1)
              end if
            else
              call pushcontrol2b(2)
            end if
            if (dwt*dwtm1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs9 = dwt
              else
                abs9 = -dwt
              end if
              if (dwtm1 .ge. 0.) then
                abs21 = dwtm1
              else
                abs21 = -dwtm1
              end if
              if (abs9 .lt. abs21) then
                dwti = dwti - half*dwt
                call pushcontrol2b(0)
              else
                dwti = dwti - half*dwtm1
                call pushcontrol2b(1)
              end if
            else
              call pushcontrol2b(2)
            end if
          else
! 1st order upwind scheme.
            dwti = w(i, j, k, jj) - w(i-1, j, k, jj)
            call pushcontrol2b(3)
          end if
          uud = uud - dwti*scratchd(i, j, k, idvt+ii-1)
          dwtid = -(uu*scratchd(i, j, k, idvt+ii-1))
          call popcontrol2b(branch)
          if (branch .lt. 2) then
            if (branch .eq. 0) then
              dwtd = -(half*dwtid)
              dwtm1d = 0.0_8
            else
              dwtm1d = -(half*dwtid)
              dwtd = 0.0_8
            end if
          else if (branch .eq. 2) then
            dwtd = 0.0_8
            dwtm1d = 0.0_8
          else
            wd(i, j, k, jj) = wd(i, j, k, jj) + dwtid
            wd(i-1, j, k, jj) = wd(i-1, j, k, jj) - dwtid
            goto 100
          end if
          call popcontrol2b(branch)
          if (branch .eq. 0) then
            dwtd = dwtd + half*dwtid
            dwtp1d = 0.0_8
          else if (branch .eq. 1) then
            dwtp1d = half*dwtid
          else
            dwtp1d = 0.0_8
          end if
          dwtd = dwtd + dwtid
          wd(i+1, j, k, jj) = wd(i+1, j, k, jj) + dwtp1d
          wd(i, j, k, jj) = wd(i, j, k, jj) + dwtd - dwtp1d
          wd(i-1, j, k, jj) = wd(i-1, j, k, jj) + dwtm1d - dwtd
          wd(i-2, j, k, jj) = wd(i-2, j, k, jj) - dwtm1d
 100    continue
      else
        uud = 0.0_8
!$bwd-of ii-loop 
        do 110 ii=1,nadv
! set the value of jj such that it corresponds to the
! turbulent entry in w.
          jj = ii + offset
! check whether a first or a second order discretization
! must be used.
          if (secondord) then
! second order; store the three differences for the
! discretization of the derivative in i-direction.
            dwtm1 = w(i, j, k, jj) - w(i-1, j, k, jj)
            dwt = w(i+1, j, k, jj) - w(i, j, k, jj)
            dwtp1 = w(i+2, j, k, jj) - w(i+1, j, k, jj)
! construct the derivative in this cell center. this is
! the first order upwind derivative with two nonlinear
! corrections.
            dwti = dwt
            if (dwt*dwtp1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs10 = dwt
              else
                abs10 = -dwt
              end if
              if (dwtp1 .ge. 0.) then
                abs22 = dwtp1
              else
                abs22 = -dwtp1
              end if
              if (abs10 .lt. abs22) then
                dwti = dwti - half*dwt
                call pushcontrol2b(0)
              else
                dwti = dwti - half*dwtp1
                call pushcontrol2b(1)
              end if
            else
              call pushcontrol2b(2)
            end if
            if (dwt*dwtm1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs11 = dwt
              else
                abs11 = -dwt
              end if
              if (dwtm1 .ge. 0.) then
                abs23 = dwtm1
              else
                abs23 = -dwtm1
              end if
              if (abs11 .lt. abs23) then
                dwti = dwti + half*dwt
                call pushcontrol2b(0)
              else
                dwti = dwti + half*dwtm1
                call pushcontrol2b(1)
              end if
            else
              call pushcontrol2b(2)
            end if
          else
! 1st order upwind scheme.
            dwti = w(i+1, j, k, jj) - w(i, j, k, jj)
            call pushcontrol2b(3)
          end if
          uud = uud - dwti*scratchd(i, j, k, idvt+ii-1)
          dwtid = -(uu*scratchd(i, j, k, idvt+ii-1))
          call popcontrol2b(branch)
          if (branch .lt. 2) then
            if (branch .eq. 0) then
              dwtd = half*dwtid
              dwtm1d = 0.0_8
            else
              dwtm1d = half*dwtid
              dwtd = 0.0_8
            end if
          else if (branch .eq. 2) then
            dwtd = 0.0_8
            dwtm1d = 0.0_8
          else
            wd(i+1, j, k, jj) = wd(i+1, j, k, jj) + dwtid
            wd(i, j, k, jj) = wd(i, j, k, jj) - dwtid
            goto 110
          end if
          call popcontrol2b(branch)
          if (branch .eq. 0) then
            dwtd = dwtd - half*dwtid
            dwtp1d = 0.0_8
          else if (branch .eq. 1) then
            dwtp1d = -(half*dwtid)
          else
            dwtp1d = 0.0_8
          end if
          dwtd = dwtd + dwtid
          wd(i+2, j, k, jj) = wd(i+2, j, k, jj) + dwtp1d
          wd(i+1, j, k, jj) = wd(i+1, j, k, jj) + dwtd - dwtp1d
          wd(i, j, k, jj) = wd(i, j, k, jj) + dwtm1d - dwtd
          wd(i-1, j, k, jj) = wd(i-1, j, k, jj) - dwtm1d
 110    continue
      end if
      xad = w(i, j, k, ivx)*uud
      wd(i, j, k, ivx) = wd(i, j, k, ivx) + xa*uud
      yad = w(i, j, k, ivy)*uud
      wd(i, j, k, ivy) = wd(i, j, k, ivy) + ya*uud
      zad = w(i, j, k, ivz)*uud
      wd(i, j, k, ivz) = wd(i, j, k, ivz) + za*uud
      qsd = qsd - uud
      sid(i, j, k, 3) = sid(i, j, k, 3) + voli*zad
      sid(i-1, j, k, 3) = sid(i-1, j, k, 3) + voli*zad
      volid = (si(i, j, k, 3)+si(i-1, j, k, 3))*zad + (si(i, j, k, 2)+si&
&       (i-1, j, k, 2))*yad + (si(i, j, k, 1)+si(i-1, j, k, 1))*xad
      sid(i, j, k, 2) = sid(i, j, k, 2) + voli*yad
      sid(i-1, j, k, 2) = sid(i-1, j, k, 2) + voli*yad
      sid(i, j, k, 1) = sid(i, j, k, 1) + voli*xad
      sid(i-1, j, k, 1) = sid(i-1, j, k, 1) + voli*xad
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        sfaceid(i, j, k) = sfaceid(i, j, k) + voli*qsd
        sfaceid(i-1, j, k) = sfaceid(i-1, j, k) + voli*qsd
        volid = volid + (sfacei(i, j, k)+sfacei(i-1, j, k))*qsd
        qsd = 0.0_8
      end if
      vold(i, j, k) = vold(i, j, k) - half*volid/vol(i, j, k)**2
    end do
    qsd = 0.0_8
    qs = zero
    qsd = 0.0_8
!$bwd-of ii-loop 
    do iii=0,nx*ny*nz-1
      i = mod(iii, nx) + 2
      j = mod(iii/nx, ny) + 2
      k = iii/(nx*ny) + 2
! compute the grid velocity if present.
! it is taken as the average of j and j-1,
      voli = half/vol(i, j, k)
      if (addgridvelocities) then
        qs = (sfacej(i, j, k)+sfacej(i, j-1, k))*voli
        call pushcontrol1b(0)
      else
        call pushcontrol1b(1)
      end if
! compute the normal velocity, where the normal direction
! is taken as the average of faces j and j-1.
      xa = (sj(i, j, k, 1)+sj(i, j-1, k, 1))*voli
      ya = (sj(i, j, k, 2)+sj(i, j-1, k, 2))*voli
      za = (sj(i, j, k, 3)+sj(i, j-1, k, 3))*voli
      uu = xa*w(i, j, k, ivx) + ya*w(i, j, k, ivy) + za*w(i, j, k, ivz) &
&       - qs
! determine the situation we are having here, i.e. positive
! or negative normal velocity.
      if (uu .gt. zero) then
        uud = 0.0_8
!$bwd-of ii-loop 
        do 120 ii=1,nadv
! set the value of jj such that it corresponds to the
! turbulent entry in w.
          jj = ii + offset
! check whether a first or a second order discretization
! must be used.
          if (secondord) then
! second order; store the three differences for the
! discretization of the derivative in j-direction.
            dwtm1 = w(i, j-1, k, jj) - w(i, j-2, k, jj)
            dwt = w(i, j, k, jj) - w(i, j-1, k, jj)
            dwtp1 = w(i, j+1, k, jj) - w(i, j, k, jj)
! construct the derivative in this cell center. this is
! the first order upwind derivative with two nonlinear
! corrections.
            dwtj = dwt
            if (dwt*dwtp1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs4 = dwt
              else
                abs4 = -dwt
              end if
              if (dwtp1 .ge. 0.) then
                abs16 = dwtp1
              else
                abs16 = -dwtp1
              end if
              if (abs4 .lt. abs16) then
                dwtj = dwtj + half*dwt
                call pushcontrol2b(0)
              else
                dwtj = dwtj + half*dwtp1
                call pushcontrol2b(1)
              end if
            else
              call pushcontrol2b(2)
            end if
            if (dwt*dwtm1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs5 = dwt
              else
                abs5 = -dwt
              end if
              if (dwtm1 .ge. 0.) then
                abs17 = dwtm1
              else
                abs17 = -dwtm1
              end if
              if (abs5 .lt. abs17) then
                dwtj = dwtj - half*dwt
                call pushcontrol2b(0)
              else
                dwtj = dwtj - half*dwtm1
                call pushcontrol2b(1)
              end if
            else
              call pushcontrol2b(2)
            end if
          else
! 1st order upwind scheme.
            dwtj = w(i, j, k, jj) - w(i, j-1, k, jj)
            call pushcontrol2b(3)
          end if
          uud = uud - dwtj*scratchd(i, j, k, idvt+ii-1)
          dwtjd = -(uu*scratchd(i, j, k, idvt+ii-1))
          call popcontrol2b(branch)
          if (branch .lt. 2) then
            if (branch .eq. 0) then
              dwtd = -(half*dwtjd)
              dwtm1d = 0.0_8
            else
              dwtm1d = -(half*dwtjd)
              dwtd = 0.0_8
            end if
          else if (branch .eq. 2) then
            dwtd = 0.0_8
            dwtm1d = 0.0_8
          else
            wd(i, j, k, jj) = wd(i, j, k, jj) + dwtjd
            wd(i, j-1, k, jj) = wd(i, j-1, k, jj) - dwtjd
            goto 120
          end if
          call popcontrol2b(branch)
          if (branch .eq. 0) then
            dwtd = dwtd + half*dwtjd
            dwtp1d = 0.0_8
          else if (branch .eq. 1) then
            dwtp1d = half*dwtjd
          else
            dwtp1d = 0.0_8
          end if
          dwtd = dwtd + dwtjd
          wd(i, j+1, k, jj) = wd(i, j+1, k, jj) + dwtp1d
          wd(i, j, k, jj) = wd(i, j, k, jj) + dwtd - dwtp1d
          wd(i, j-1, k, jj) = wd(i, j-1, k, jj) + dwtm1d - dwtd
          wd(i, j-2, k, jj) = wd(i, j-2, k, jj) - dwtm1d
 120    continue
      else
        uud = 0.0_8
!$bwd-of ii-loop 
        do 130 ii=1,nadv
! set the value of jj such that it corresponds to the
! turbulent entry in w.
          jj = ii + offset
! check whether a first or a second order discretization
! must be used.
          if (secondord) then
! store the three differences for the discretization of
! the derivative in j-direction.
            dwtm1 = w(i, j, k, jj) - w(i, j-1, k, jj)
            dwt = w(i, j+1, k, jj) - w(i, j, k, jj)
            dwtp1 = w(i, j+2, k, jj) - w(i, j+1, k, jj)
! construct the derivative in this cell center. this is
! the first order upwind derivative with two nonlinear
! corrections.
            dwtj = dwt
            if (dwt*dwtp1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs6 = dwt
              else
                abs6 = -dwt
              end if
              if (dwtp1 .ge. 0.) then
                abs18 = dwtp1
              else
                abs18 = -dwtp1
              end if
              if (abs6 .lt. abs18) then
                dwtj = dwtj - half*dwt
                call pushcontrol2b(0)
              else
                dwtj = dwtj - half*dwtp1
                call pushcontrol2b(1)
              end if
            else
              call pushcontrol2b(2)
            end if
            if (dwt*dwtm1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs7 = dwt
              else
                abs7 = -dwt
              end if
              if (dwtm1 .ge. 0.) then
                abs19 = dwtm1
              else
                abs19 = -dwtm1
              end if
              if (abs7 .lt. abs19) then
                dwtj = dwtj + half*dwt
                call pushcontrol2b(0)
              else
                dwtj = dwtj + half*dwtm1
                call pushcontrol2b(1)
              end if
            else
              call pushcontrol2b(2)
            end if
          else
! 1st order upwind scheme.
            dwtj = w(i, j+1, k, jj) - w(i, j, k, jj)
            call pushcontrol2b(3)
          end if
          uud = uud - dwtj*scratchd(i, j, k, idvt+ii-1)
          dwtjd = -(uu*scratchd(i, j, k, idvt+ii-1))
          call popcontrol2b(branch)
          if (branch .lt. 2) then
            if (branch .eq. 0) then
              dwtd = half*dwtjd
              dwtm1d = 0.0_8
            else
              dwtm1d = half*dwtjd
              dwtd = 0.0_8
            end if
          else if (branch .eq. 2) then
            dwtd = 0.0_8
            dwtm1d = 0.0_8
          else
            wd(i, j+1, k, jj) = wd(i, j+1, k, jj) + dwtjd
            wd(i, j, k, jj) = wd(i, j, k, jj) - dwtjd
            goto 130
          end if
          call popcontrol2b(branch)
          if (branch .eq. 0) then
            dwtd = dwtd - half*dwtjd
            dwtp1d = 0.0_8
          else if (branch .eq. 1) then
            dwtp1d = -(half*dwtjd)
          else
            dwtp1d = 0.0_8
          end if
          dwtd = dwtd + dwtjd
          wd(i, j+2, k, jj) = wd(i, j+2, k, jj) + dwtp1d
          wd(i, j+1, k, jj) = wd(i, j+1, k, jj) + dwtd - dwtp1d
          wd(i, j, k, jj) = wd(i, j, k, jj) + dwtm1d - dwtd
          wd(i, j-1, k, jj) = wd(i, j-1, k, jj) - dwtm1d
 130    continue
      end if
      xad = w(i, j, k, ivx)*uud
      wd(i, j, k, ivx) = wd(i, j, k, ivx) + xa*uud
      yad = w(i, j, k, ivy)*uud
      wd(i, j, k, ivy) = wd(i, j, k, ivy) + ya*uud
      zad = w(i, j, k, ivz)*uud
      wd(i, j, k, ivz) = wd(i, j, k, ivz) + za*uud
      qsd = qsd - uud
      sjd(i, j, k, 3) = sjd(i, j, k, 3) + voli*zad
      sjd(i, j-1, k, 3) = sjd(i, j-1, k, 3) + voli*zad
      volid = (sj(i, j, k, 3)+sj(i, j-1, k, 3))*zad + (sj(i, j, k, 2)+sj&
&       (i, j-1, k, 2))*yad + (sj(i, j, k, 1)+sj(i, j-1, k, 1))*xad
      sjd(i, j, k, 2) = sjd(i, j, k, 2) + voli*yad
      sjd(i, j-1, k, 2) = sjd(i, j-1, k, 2) + voli*yad
      sjd(i, j, k, 1) = sjd(i, j, k, 1) + voli*xad
      sjd(i, j-1, k, 1) = sjd(i, j-1, k, 1) + voli*xad
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        sfacejd(i, j, k) = sfacejd(i, j, k) + voli*qsd
        sfacejd(i, j-1, k) = sfacejd(i, j-1, k) + voli*qsd
        volid = volid + (sfacej(i, j, k)+sfacej(i, j-1, k))*qsd
        qsd = 0.0_8
      end if
      vold(i, j, k) = vold(i, j, k) - half*volid/vol(i, j, k)**2
    end do
    qsd = 0.0_8
! initialize the grid velocity to zero. this value will be used
! if the block is not moving.
    qs = zero
    qsd = 0.0_8
!$bwd-of ii-loop 
    do iii=0,nx*ny*nz-1
      i = mod(iii, nx) + 2
      j = mod(iii/nx, ny) + 2
      k = iii/(nx*ny) + 2
! compute the grid velocity if present.
! it is taken as the average of k and k-1,
      voli = half/vol(i, j, k)
      if (addgridvelocities) then
        qs = (sfacek(i, j, k)+sfacek(i, j, k-1))*voli
        call pushcontrol1b(0)
      else
        call pushcontrol1b(1)
      end if
! compute the normal velocity, where the normal direction
! is taken as the average of faces k and k-1.
      xa = (sk(i, j, k, 1)+sk(i, j, k-1, 1))*voli
      ya = (sk(i, j, k, 2)+sk(i, j, k-1, 2))*voli
      za = (sk(i, j, k, 3)+sk(i, j, k-1, 3))*voli
      uu = xa*w(i, j, k, ivx) + ya*w(i, j, k, ivy) + za*w(i, j, k, ivz) &
&       - qs
! this term has unit: velocity/length
! determine the situation we are having here, i.e. positive
! or negative normal velocity.
      if (uu .gt. zero) then
        uud = 0.0_8
!$bwd-of ii-loop 
        do 140 ii=1,nadv
! set the value of jj such that it corresponds to the
! turbulent entry in w.
          jj = ii + offset
! check whether a first or a second order discretization
! must be used.
          if (secondord) then
! second order; store the three differences for the
! discretization of the derivative in k-direction.
            dwtm1 = w(i, j, k-1, jj) - w(i, j, k-2, jj)
            dwt = w(i, j, k, jj) - w(i, j, k-1, jj)
            dwtp1 = w(i, j, k+1, jj) - w(i, j, k, jj)
! construct the derivative in this cell center. this
! is the first order upwind derivative with two
! nonlinear corrections.
            dwtk = dwt
            if (dwt*dwtp1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs0 = dwt
              else
                abs0 = -dwt
              end if
              if (dwtp1 .ge. 0.) then
                abs12 = dwtp1
              else
                abs12 = -dwtp1
              end if
              if (abs0 .lt. abs12) then
                dwtk = dwtk + half*dwt
                call pushcontrol2b(0)
              else
                dwtk = dwtk + half*dwtp1
                call pushcontrol2b(1)
              end if
            else
              call pushcontrol2b(2)
            end if
            if (dwt*dwtm1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs1 = dwt
              else
                abs1 = -dwt
              end if
              if (dwtm1 .ge. 0.) then
                abs13 = dwtm1
              else
                abs13 = -dwtm1
              end if
              if (abs1 .lt. abs13) then
                dwtk = dwtk - half*dwt
                call pushcontrol2b(0)
              else
                dwtk = dwtk - half*dwtm1
                call pushcontrol2b(1)
              end if
            else
              call pushcontrol2b(2)
            end if
          else
! 1st order upwind scheme.
            dwtk = w(i, j, k, jj) - w(i, j, k-1, jj)
            call pushcontrol2b(3)
          end if
          uud = uud - dwtk*scratchd(i, j, k, idvt+ii-1)
          dwtkd = -(uu*scratchd(i, j, k, idvt+ii-1))
          call popcontrol2b(branch)
          if (branch .lt. 2) then
            if (branch .eq. 0) then
              dwtd = -(half*dwtkd)
              dwtm1d = 0.0_8
            else
              dwtm1d = -(half*dwtkd)
              dwtd = 0.0_8
            end if
          else if (branch .eq. 2) then
            dwtd = 0.0_8
            dwtm1d = 0.0_8
          else
            wd(i, j, k, jj) = wd(i, j, k, jj) + dwtkd
            wd(i, j, k-1, jj) = wd(i, j, k-1, jj) - dwtkd
            goto 140
          end if
          call popcontrol2b(branch)
          if (branch .eq. 0) then
            dwtd = dwtd + half*dwtkd
            dwtp1d = 0.0_8
          else if (branch .eq. 1) then
            dwtp1d = half*dwtkd
          else
            dwtp1d = 0.0_8
          end if
          dwtd = dwtd + dwtkd
          wd(i, j, k+1, jj) = wd(i, j, k+1, jj) + dwtp1d
          wd(i, j, k, jj) = wd(i, j, k, jj) + dwtd - dwtp1d
          wd(i, j, k-1, jj) = wd(i, j, k-1, jj) + dwtm1d - dwtd
          wd(i, j, k-2, jj) = wd(i, j, k-2, jj) - dwtm1d
 140    continue
      else
        uud = 0.0_8
!$bwd-of ii-loop 
        do 150 ii=1,nadv
! set the value of jj such that it corresponds to the
! turbulent entry in w.
          jj = ii + offset
! check whether a first or a second order discretization
! must be used.
          if (secondord) then
! store the three differences for the discretization of
! the derivative in k-direction.
            dwtm1 = w(i, j, k, jj) - w(i, j, k-1, jj)
            dwt = w(i, j, k+1, jj) - w(i, j, k, jj)
            dwtp1 = w(i, j, k+2, jj) - w(i, j, k+1, jj)
! construct the derivative in this cell center. this is
! the first order upwind derivative with two nonlinear
! corrections.
            dwtk = dwt
            if (dwt*dwtp1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs2 = dwt
              else
                abs2 = -dwt
              end if
              if (dwtp1 .ge. 0.) then
                abs14 = dwtp1
              else
                abs14 = -dwtp1
              end if
              if (abs2 .lt. abs14) then
                dwtk = dwtk - half*dwt
                call pushcontrol2b(0)
              else
                dwtk = dwtk - half*dwtp1
                call pushcontrol2b(1)
              end if
            else
              call pushcontrol2b(2)
            end if
            if (dwt*dwtm1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs3 = dwt
              else
                abs3 = -dwt
              end if
              if (dwtm1 .ge. 0.) then
                abs15 = dwtm1
              else
                abs15 = -dwtm1
              end if
              if (abs3 .lt. abs15) then
                dwtk = dwtk + half*dwt
                call pushcontrol2b(0)
              else
                dwtk = dwtk + half*dwtm1
                call pushcontrol2b(1)
              end if
            else
              call pushcontrol2b(2)
            end if
          else
! 1st order upwind scheme.
            dwtk = w(i, j, k+1, jj) - w(i, j, k, jj)
            call pushcontrol2b(3)
          end if
          uud = uud - dwtk*scratchd(i, j, k, idvt+ii-1)
          dwtkd = -(uu*scratchd(i, j, k, idvt+ii-1))
          call popcontrol2b(branch)
          if (branch .lt. 2) then
            if (branch .eq. 0) then
              dwtd = half*dwtkd
              dwtm1d = 0.0_8
            else
              dwtm1d = half*dwtkd
              dwtd = 0.0_8
            end if
          else if (branch .eq. 2) then
            dwtd = 0.0_8
            dwtm1d = 0.0_8
          else
            wd(i, j, k+1, jj) = wd(i, j, k+1, jj) + dwtkd
            wd(i, j, k, jj) = wd(i, j, k, jj) - dwtkd
            goto 150
          end if
          call popcontrol2b(branch)
          if (branch .eq. 0) then
            dwtd = dwtd - half*dwtkd
            dwtp1d = 0.0_8
          else if (branch .eq. 1) then
            dwtp1d = -(half*dwtkd)
          else
            dwtp1d = 0.0_8
          end if
          dwtd = dwtd + dwtkd
          wd(i, j, k+2, jj) = wd(i, j, k+2, jj) + dwtp1d
          wd(i, j, k+1, jj) = wd(i, j, k+1, jj) + dwtd - dwtp1d
          wd(i, j, k, jj) = wd(i, j, k, jj) + dwtm1d - dwtd
          wd(i, j, k-1, jj) = wd(i, j, k-1, jj) - dwtm1d
 150    continue
      end if
      xad = w(i, j, k, ivx)*uud
      wd(i, j, k, ivx) = wd(i, j, k, ivx) + xa*uud
      yad = w(i, j, k, ivy)*uud
      wd(i, j, k, ivy) = wd(i, j, k, ivy) + ya*uud
      zad = w(i, j, k, ivz)*uud
      wd(i, j, k, ivz) = wd(i, j, k, ivz) + za*uud
      qsd = qsd - uud
      skd(i, j, k, 3) = skd(i, j, k, 3) + voli*zad
      skd(i, j, k-1, 3) = skd(i, j, k-1, 3) + voli*zad
      volid = (sk(i, j, k, 3)+sk(i, j, k-1, 3))*zad + (sk(i, j, k, 2)+sk&
&       (i, j, k-1, 2))*yad + (sk(i, j, k, 1)+sk(i, j, k-1, 1))*xad
      skd(i, j, k, 2) = skd(i, j, k, 2) + voli*yad
      skd(i, j, k-1, 2) = skd(i, j, k-1, 2) + voli*yad
      skd(i, j, k, 1) = skd(i, j, k, 1) + voli*xad
      skd(i, j, k-1, 1) = skd(i, j, k-1, 1) + voli*xad
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        sfacekd(i, j, k) = sfacekd(i, j, k) + voli*qsd
        sfacekd(i, j, k-1) = sfacekd(i, j, k-1) + voli*qsd
        volid = volid + (sfacek(i, j, k)+sfacek(i, j, k-1))*qsd
        qsd = 0.0_8
      end if
      vold(i, j, k) = vold(i, j, k) - half*volid/vol(i, j, k)**2
    end do
  end subroutine turbadvection_b
 
  subroutine turbadvection(madv, nadv, offset, qq)
!
!       turbadvection discretizes the advection part of the turbulent
!       transport equations. as the advection part is the same for all
!       models, this generic routine can be used. both the
!       discretization and the central jacobian are computed in this
!       subroutine. the former can either be 1st or 2nd order
!       accurate; the latter is always based on the 1st order upwind
!       discretization. when the discretization must be second order
!       accurate, the fully upwind (kappa = -1) scheme in combination
!       with the minmod limiter is used.
!       only nadv equations are treated, while the actual system has
!       size madv. the reason is that some equations for some
!       turbulence equations do not have an advection part, e.g. the
!       f equation in the v2-f model. the argument offset indicates
!       the offset in the w vector where this subsystem starts. as a
!       consequence it is assumed that the indices of the current
!       subsystem are contiguous, e.g. if a 2*2 system is solved the
!       last index in w is offset+1 and offset+2 respectively.
!
    use constants
    use blockpointers, only : nx, ny, nz, il, jl, kl, vol, sfacei, &
&   sfacej, sfacek, w, si, sj, sk, addgridvelocities, bmti1, bmti2, &
&   bmtj1, bmtj2, bmtk1, bmtk2, scratch
    use inputdiscretization, only : orderturb
    use iteration, only : groundlevel
    use turbmod, only : secondord
    implicit none
!
!      subroutine arguments.
!
    integer(kind=inttype), intent(in) :: nadv, madv, offset
    real(kind=realtype), dimension(2:il, 2:jl, 2:kl, madv, madv), &
&   intent(inout) :: qq
!
!      local variables.
!
    integer(kind=inttype) :: i, j, k, ii, jj, kk, iii
    real(kind=realtype) :: qs, voli, xa, ya, za
    real(kind=realtype) :: uu, dwt, dwtm1, dwtp1, dwti, dwtj, dwtk
    real(kind=realtype), dimension(madv) :: impl
    intrinsic mod
    intrinsic abs
    real(kind=realtype) :: abs0
    real(kind=realtype) :: abs1
    real(kind=realtype) :: abs2
    real(kind=realtype) :: abs3
    real(kind=realtype) :: abs4
    real(kind=realtype) :: abs5
    real(kind=realtype) :: abs6
    real(kind=realtype) :: abs7
    real(kind=realtype) :: abs8
    real(kind=realtype) :: abs9
    real(kind=realtype) :: abs10
    real(kind=realtype) :: abs11
    real(kind=realtype) :: abs12
    real(kind=realtype) :: abs13
    real(kind=realtype) :: abs14
    real(kind=realtype) :: abs15
    real(kind=realtype) :: abs16
    real(kind=realtype) :: abs17
    real(kind=realtype) :: abs18
    real(kind=realtype) :: abs19
    real(kind=realtype) :: abs20
    real(kind=realtype) :: abs21
    real(kind=realtype) :: abs22
    real(kind=realtype) :: abs23
! determine whether or not a second order discretization for the
! advective terms must be used.
    secondord = .false.
    if (groundlevel .eq. 1_inttype .and. orderturb .eq. secondorder) &
&     secondord = .true.
!$ad checkpoint-start
! initialize the grid velocity to zero. this value will be used
! if the block is not moving.
    qs = zero
!$ad ii-loop
!
!       upwind discretization of the convective term in k (zeta)
!       direction. either the 1st order upwind or the second order
!       fully upwind interpolation scheme, kappa = -1, is used in
!       combination with the minmod limiter.
!       the possible grid velocity must be taken into account.
!
    do iii=0,nx*ny*nz-1
      i = mod(iii, nx) + 2
      j = mod(iii/nx, ny) + 2
      k = iii/(nx*ny) + 2
! compute the grid velocity if present.
! it is taken as the average of k and k-1,
      voli = half/vol(i, j, k)
      if (addgridvelocities) qs = (sfacek(i, j, k)+sfacek(i, j, k-1))*&
&         voli
! compute the normal velocity, where the normal direction
! is taken as the average of faces k and k-1.
      xa = (sk(i, j, k, 1)+sk(i, j, k-1, 1))*voli
      ya = (sk(i, j, k, 2)+sk(i, j, k-1, 2))*voli
      za = (sk(i, j, k, 3)+sk(i, j, k-1, 3))*voli
      uu = xa*w(i, j, k, ivx) + ya*w(i, j, k, ivy) + za*w(i, j, k, ivz) &
&       - qs
! this term has unit: velocity/length
! determine the situation we are having here, i.e. positive
! or negative normal velocity.
      if (uu .gt. zero) then
!$ad ii-loop
! velocity has a component in positive k-direction.
! loop over the number of advection equations.
        do ii=1,nadv
! set the value of jj such that it corresponds to the
! turbulent entry in w.
          jj = ii + offset
! check whether a first or a second order discretization
! must be used.
          if (secondord) then
! second order; store the three differences for the
! discretization of the derivative in k-direction.
            dwtm1 = w(i, j, k-1, jj) - w(i, j, k-2, jj)
            dwt = w(i, j, k, jj) - w(i, j, k-1, jj)
            dwtp1 = w(i, j, k+1, jj) - w(i, j, k, jj)
! construct the derivative in this cell center. this
! is the first order upwind derivative with two
! nonlinear corrections.
            dwtk = dwt
            if (dwt*dwtp1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs0 = dwt
              else
                abs0 = -dwt
              end if
              if (dwtp1 .ge. 0.) then
                abs12 = dwtp1
              else
                abs12 = -dwtp1
              end if
              if (abs0 .lt. abs12) then
                dwtk = dwtk + half*dwt
              else
                dwtk = dwtk + half*dwtp1
              end if
            end if
            if (dwt*dwtm1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs1 = dwt
              else
                abs1 = -dwt
              end if
              if (dwtm1 .ge. 0.) then
                abs13 = dwtm1
              else
                abs13 = -dwtm1
              end if
              if (abs1 .lt. abs13) then
                dwtk = dwtk - half*dwt
              else
                dwtk = dwtk - half*dwtm1
              end if
            end if
          else
! 1st order upwind scheme.
            dwtk = w(i, j, k, jj) - w(i, j, k-1, jj)
          end if
! update the residual. the convective term must be
! substracted, because it appears on the other side of
! the equation as the source and viscous terms.
! uu*dwtk = (v.dot.face_normal)*delta(nutilde)/delta(x)
          scratch(i, j, k, idvt+ii-1) = scratch(i, j, k, idvt+ii-1) - uu&
&           *dwtk
        end do
      else
!$ad ii-loop
! velocity has a component in negative k-direction.
! loop over the number of advection equations
        do ii=1,nadv
! set the value of jj such that it corresponds to the
! turbulent entry in w.
          jj = ii + offset
! check whether a first or a second order discretization
! must be used.
          if (secondord) then
! store the three differences for the discretization of
! the derivative in k-direction.
            dwtm1 = w(i, j, k, jj) - w(i, j, k-1, jj)
            dwt = w(i, j, k+1, jj) - w(i, j, k, jj)
            dwtp1 = w(i, j, k+2, jj) - w(i, j, k+1, jj)
! construct the derivative in this cell center. this is
! the first order upwind derivative with two nonlinear
! corrections.
            dwtk = dwt
            if (dwt*dwtp1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs2 = dwt
              else
                abs2 = -dwt
              end if
              if (dwtp1 .ge. 0.) then
                abs14 = dwtp1
              else
                abs14 = -dwtp1
              end if
              if (abs2 .lt. abs14) then
                dwtk = dwtk - half*dwt
              else
                dwtk = dwtk - half*dwtp1
              end if
            end if
            if (dwt*dwtm1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs3 = dwt
              else
                abs3 = -dwt
              end if
              if (dwtm1 .ge. 0.) then
                abs15 = dwtm1
              else
                abs15 = -dwtm1
              end if
              if (abs3 .lt. abs15) then
                dwtk = dwtk + half*dwt
              else
                dwtk = dwtk + half*dwtm1
              end if
            end if
          else
! 1st order upwind scheme.
            dwtk = w(i, j, k+1, jj) - w(i, j, k, jj)
          end if
! update the residual. the convective term must be
! substracted, because it appears on the other side
! of the equation as the source and viscous terms.
          scratch(i, j, k, idvt+ii-1) = scratch(i, j, k, idvt+ii-1) - uu&
&           *dwtk
! update the central jacobian. first the term which is
! always present, i.e. -uu.
        end do
      end if
    end do
!$ad checkpoint-end
!
!       upwind discretization of the convective term in j (eta)
!       direction. either the 1st order upwind or the second order
!       fully upwind interpolation scheme, kappa = -1, is used in
!       combination with the minmod limiter.
!       the possible grid velocity must be taken into account.
!
    continue
!$ad checkpoint-start
    qs = zero
!$ad ii-loop
    do iii=0,nx*ny*nz-1
      i = mod(iii, nx) + 2
      j = mod(iii/nx, ny) + 2
      k = iii/(nx*ny) + 2
! compute the grid velocity if present.
! it is taken as the average of j and j-1,
      voli = half/vol(i, j, k)
      if (addgridvelocities) qs = (sfacej(i, j, k)+sfacej(i, j-1, k))*&
&         voli
! compute the normal velocity, where the normal direction
! is taken as the average of faces j and j-1.
      xa = (sj(i, j, k, 1)+sj(i, j-1, k, 1))*voli
      ya = (sj(i, j, k, 2)+sj(i, j-1, k, 2))*voli
      za = (sj(i, j, k, 3)+sj(i, j-1, k, 3))*voli
      uu = xa*w(i, j, k, ivx) + ya*w(i, j, k, ivy) + za*w(i, j, k, ivz) &
&       - qs
! determine the situation we are having here, i.e. positive
! or negative normal velocity.
      if (uu .gt. zero) then
!$ad ii-loop
! velocity has a component in positive j-direction.
! loop over the number of advection equations.
        do ii=1,nadv
! set the value of jj such that it corresponds to the
! turbulent entry in w.
          jj = ii + offset
! check whether a first or a second order discretization
! must be used.
          if (secondord) then
! second order; store the three differences for the
! discretization of the derivative in j-direction.
            dwtm1 = w(i, j-1, k, jj) - w(i, j-2, k, jj)
            dwt = w(i, j, k, jj) - w(i, j-1, k, jj)
            dwtp1 = w(i, j+1, k, jj) - w(i, j, k, jj)
! construct the derivative in this cell center. this is
! the first order upwind derivative with two nonlinear
! corrections.
            dwtj = dwt
            if (dwt*dwtp1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs4 = dwt
              else
                abs4 = -dwt
              end if
              if (dwtp1 .ge. 0.) then
                abs16 = dwtp1
              else
                abs16 = -dwtp1
              end if
              if (abs4 .lt. abs16) then
                dwtj = dwtj + half*dwt
              else
                dwtj = dwtj + half*dwtp1
              end if
            end if
            if (dwt*dwtm1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs5 = dwt
              else
                abs5 = -dwt
              end if
              if (dwtm1 .ge. 0.) then
                abs17 = dwtm1
              else
                abs17 = -dwtm1
              end if
              if (abs5 .lt. abs17) then
                dwtj = dwtj - half*dwt
              else
                dwtj = dwtj - half*dwtm1
              end if
            end if
          else
! 1st order upwind scheme.
            dwtj = w(i, j, k, jj) - w(i, j-1, k, jj)
          end if
! update the residual. the convective term must be
! substracted, because it appears on the other side of
! the equation as the source and viscous terms.
          scratch(i, j, k, idvt+ii-1) = scratch(i, j, k, idvt+ii-1) - uu&
&           *dwtj
! update the central jacobian. first the term which is
! always present, i.e. uu.
        end do
      else
!$ad ii-loop
! velocity has a component in negative j-direction.
! loop over the number of advection equations.
        do ii=1,nadv
! set the value of jj such that it corresponds to the
! turbulent entry in w.
          jj = ii + offset
! check whether a first or a second order discretization
! must be used.
          if (secondord) then
! store the three differences for the discretization of
! the derivative in j-direction.
            dwtm1 = w(i, j, k, jj) - w(i, j-1, k, jj)
            dwt = w(i, j+1, k, jj) - w(i, j, k, jj)
            dwtp1 = w(i, j+2, k, jj) - w(i, j+1, k, jj)
! construct the derivative in this cell center. this is
! the first order upwind derivative with two nonlinear
! corrections.
            dwtj = dwt
            if (dwt*dwtp1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs6 = dwt
              else
                abs6 = -dwt
              end if
              if (dwtp1 .ge. 0.) then
                abs18 = dwtp1
              else
                abs18 = -dwtp1
              end if
              if (abs6 .lt. abs18) then
                dwtj = dwtj - half*dwt
              else
                dwtj = dwtj - half*dwtp1
              end if
            end if
            if (dwt*dwtm1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs7 = dwt
              else
                abs7 = -dwt
              end if
              if (dwtm1 .ge. 0.) then
                abs19 = dwtm1
              else
                abs19 = -dwtm1
              end if
              if (abs7 .lt. abs19) then
                dwtj = dwtj + half*dwt
              else
                dwtj = dwtj + half*dwtm1
              end if
            end if
          else
! 1st order upwind scheme.
            dwtj = w(i, j+1, k, jj) - w(i, j, k, jj)
          end if
! update the residual. the convective term must be
! substracted, because it appears on the other side
! of the equation as the source and viscous terms.
          scratch(i, j, k, idvt+ii-1) = scratch(i, j, k, idvt+ii-1) - uu&
&           *dwtj
! update the central jacobian. first the term which is
! always present, i.e. -uu.
        end do
      end if
    end do
!$ad checkpoint-end
!
!       upwind discretization of the convective term in i (xi)
!       direction. either the 1st order upwind or the second order
!       fully upwind interpolation scheme, kappa = -1, is used in
!       combination with the minmod limiter.
!       the possible grid velocity must be taken into account.
!
    continue
!$ad checkpoint-start
    qs = zero
!$ad ii-loop
    do iii=0,nx*ny*nz-1
      i = mod(iii, nx) + 2
      j = mod(iii/nx, ny) + 2
      k = iii/(nx*ny) + 2
! compute the grid velocity if present.
! it is taken as the average of i and i-1,
      voli = half/vol(i, j, k)
      if (addgridvelocities) qs = (sfacei(i, j, k)+sfacei(i-1, j, k))*&
&         voli
! compute the normal velocity, where the normal direction
! is taken as the average of faces i and i-1.
      xa = (si(i, j, k, 1)+si(i-1, j, k, 1))*voli
      ya = (si(i, j, k, 2)+si(i-1, j, k, 2))*voli
      za = (si(i, j, k, 3)+si(i-1, j, k, 3))*voli
      uu = xa*w(i, j, k, ivx) + ya*w(i, j, k, ivy) + za*w(i, j, k, ivz) &
&       - qs
! determine the situation we are having here, i.e. positive
! or negative normal velocity.
      if (uu .gt. zero) then
!$ad ii-loop
! velocity has a component in positive i-direction.
! loop over the number of advection equations.
        do ii=1,nadv
! set the value of jj such that it corresponds to the
! turbulent entry in w.
          jj = ii + offset
! check whether a first or a second order discretization
! must be used.
          if (secondord) then
! second order; store the three differences for the
! discretization of the derivative in i-direction.
            dwtm1 = w(i-1, j, k, jj) - w(i-2, j, k, jj)
            dwt = w(i, j, k, jj) - w(i-1, j, k, jj)
            dwtp1 = w(i+1, j, k, jj) - w(i, j, k, jj)
! construct the derivative in this cell center. this is
! the first order upwind derivative with two nonlinear
! corrections.
            dwti = dwt
            if (dwt*dwtp1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs8 = dwt
              else
                abs8 = -dwt
              end if
              if (dwtp1 .ge. 0.) then
                abs20 = dwtp1
              else
                abs20 = -dwtp1
              end if
              if (abs8 .lt. abs20) then
                dwti = dwti + half*dwt
              else
                dwti = dwti + half*dwtp1
              end if
            end if
            if (dwt*dwtm1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs9 = dwt
              else
                abs9 = -dwt
              end if
              if (dwtm1 .ge. 0.) then
                abs21 = dwtm1
              else
                abs21 = -dwtm1
              end if
              if (abs9 .lt. abs21) then
                dwti = dwti - half*dwt
              else
                dwti = dwti - half*dwtm1
              end if
            end if
          else
! 1st order upwind scheme.
            dwti = w(i, j, k, jj) - w(i-1, j, k, jj)
          end if
! update the residual. the convective term must be
! substracted, because it appears on the other side of
! the equation as the source and viscous terms.
          scratch(i, j, k, idvt+ii-1) = scratch(i, j, k, idvt+ii-1) - uu&
&           *dwti
! update the central jacobian. first the term which is
! always present, i.e. uu.
        end do
      else
!$ad ii-loop
! velocity has a component in negative i-direction.
! loop over the number of advection equations.
        do ii=1,nadv
! set the value of jj such that it corresponds to the
! turbulent entry in w.
          jj = ii + offset
! check whether a first or a second order discretization
! must be used.
          if (secondord) then
! second order; store the three differences for the
! discretization of the derivative in i-direction.
            dwtm1 = w(i, j, k, jj) - w(i-1, j, k, jj)
            dwt = w(i+1, j, k, jj) - w(i, j, k, jj)
            dwtp1 = w(i+2, j, k, jj) - w(i+1, j, k, jj)
! construct the derivative in this cell center. this is
! the first order upwind derivative with two nonlinear
! corrections.
            dwti = dwt
            if (dwt*dwtp1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs10 = dwt
              else
                abs10 = -dwt
              end if
              if (dwtp1 .ge. 0.) then
                abs22 = dwtp1
              else
                abs22 = -dwtp1
              end if
              if (abs10 .lt. abs22) then
                dwti = dwti - half*dwt
              else
                dwti = dwti - half*dwtp1
              end if
            end if
            if (dwt*dwtm1 .gt. zero) then
              if (dwt .ge. 0.) then
                abs11 = dwt
              else
                abs11 = -dwt
              end if
              if (dwtm1 .ge. 0.) then
                abs23 = dwtm1
              else
                abs23 = -dwtm1
              end if
              if (abs11 .lt. abs23) then
                dwti = dwti + half*dwt
              else
                dwti = dwti + half*dwtm1
              end if
            end if
          else
! 1st order upwind scheme.
            dwti = w(i+1, j, k, jj) - w(i, j, k, jj)
          end if
! update the residual. the convective term must be
! substracted, because it appears on the other side
! of the equation as the source and viscous terms.
          scratch(i, j, k, idvt+ii-1) = scratch(i, j, k, idvt+ii-1) - uu&
&           *dwti
! update the central jacobian. first the term which is
! always present, i.e. -uu.
        end do
      end if
    end do
!$ad checkpoint-end
 
  end subroutine turbadvection
! ----------------------------------------------------------------------
!                                                                      |
!                    no tapenade routine below this line               |
!                                                                      |
! ----------------------------------------------------------------------
 
end module turbutils_b