turbUtils_d.f90

Go to the documentation of this file.

!        generated by tapenade     (inria, ecuador team)
!  tapenade 3.16 (develop) - 22 aug 2023 15:51
!
module turbutils_d
  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 sqrt
    real(kind=realtype) :: arg1
    real(kind=realtype) :: result1
! 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)
! 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 k=2,kl
      do j=2,jl
        do i=2,il
! 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.
          arg1 = sijsij*oijoij
          result1 = sqrt(arg1)
          scratch(i, j, k, iprod) = two*result1
        end do
      end do
    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
! 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 k=2,kl
      do j=2,jl
        do i=2,il
! 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 do
    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
! 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)
! 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 k=2,kl
      do j=2,jl
        do i=2,il
! 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 do
    end do
  end subroutine prodwmag2
 
!  differentiation of sanuknowneddyratio in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: sanuknowneddyratio
!   with respect to varying inputs: nulam
  real(kind=realtype) function sanuknowneddyratio_d(eddyratio, nulam, &
&   nulamd, sanuknowneddyratio)
!
!       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
    real(kind=realtype), intent(in) :: nulamd
!
!      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
      sanuknowneddyratio_d = 0.0_8
      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) exit
      end do
! 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_d = chi*nulamd
      sanuknowneddyratio = nulam*chi
    end if
  end function sanuknowneddyratio_d
 
  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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *rev
!   with respect to varying inputs: *w *rlv
!   rw status of diff variables: *rev:out *w:in *rlv:in
!   plus diff mem management of: rev:in w:in rlv:in
  subroutine computeeddyviscosity_d(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
      if (associated(revd)) revd = 0.0_8
      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_d(ibeg, iend, jbeg, jend, kbeg, kend)
      case default
        if (associated(revd)) revd = 0.0_8
      end select
    end if
  end subroutine computeeddyviscosity_d
 
  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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *rev
!   with respect to varying inputs: *w *rlv
!   plus diff mem management of: rev:in w:in rlv:in
  subroutine saeddyviscosity_d(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
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
! store the cv1^3; cv1 is a constant of the spalart-allmaras model.
    cv13 = rsacv1**3
    if (associated(revd)) revd = 0.0_8
! loop over the cells of this block and compute the eddy viscosity.
! do not include halo's.
    do k=kbeg,kend
      do j=jbeg,jend
        do i=ibeg,iend
          temp = w(i, j, k, irho)
          temp0 = w(i, j, k, itu1)
          rnusad = temp*wd(i, j, k, itu1) + temp0*wd(i, j, k, irho)
          rnusa = temp0*temp
          temp0 = rnusa/rlv(i, j, k)
          chid = (rnusad-temp0*rlvd(i, j, k))/rlv(i, j, k)
          chi = temp0
          chi3d = 3*chi**2*chid
          chi3 = chi**3
          temp0 = chi3/(cv13+chi3)
          fv1d = (1.0-temp0)*chi3d/(cv13+chi3)
          fv1 = temp0
          revd(i, j, k) = rnusa*fv1d + fv1*rnusad
          rev(i, j, k) = fv1*rnusa
        end do
      end do
    end do
  end subroutine saeddyviscosity_d
 
  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
! 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.
    do k=kbeg,kend
      do j=jbeg,jend
        do i=ibeg,iend
          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 do
    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 abs
    real(kind=realtype) :: x1
! loop over the cells of this block and compute the eddy viscosity.
! do not include halo's.
    do k=kbeg,kend
      do j=jbeg,jend
        do i=ibeg,iend
          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 do
    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 sqrt
    intrinsic max
    intrinsic tanh
    real(kind=realtype) :: max1
    real(kind=realtype) :: result1
    real(kind=realtype) :: arg1
! 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.
    do k=kbeg,kend
      do j=jbeg,jend
        do i=ibeg,iend
! compute the value of the function f2, which occurs in the
! eddy-viscosity computation.
          result1 = sqrt(w(i, j, k, itu1))
          t1 = two*result1/(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
          arg1 = arg2**2
          f2 = tanh(arg1)
! 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 do
    end do
  end subroutine ssteddyviscosity
 
!  differentiation of turbadvection in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *scratch
!   with respect to varying inputs: *sfacei *sfacej *sfacek *w
!                *scratch *vol *si *sj *sk
!   rw status of diff variables: *sfacei:in *sfacej:in *sfacek:in
!                *w:in *scratch:in-out *vol:in *si:in *sj:in *sk:in
!   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_d(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 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
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype) :: temp1
! 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.
! initialize the grid velocity to zero. this value will be used
! if the block is not moving.
    qs = zero
    qsd = 0.0_8
!
!       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 k=2,kl
      do j=2,jl
        do i=2,il
! compute the grid velocity if present.
! it is taken as the average of k and k-1,
          temp = half/vol(i, j, k)
          volid = -(temp*vold(i, j, k)/vol(i, j, k))
          voli = temp
          if (addgridvelocities) then
            temp = sfacek(i, j, k) + sfacek(i, j, k-1)
            qsd = voli*(sfacekd(i, j, k)+sfacekd(i, j, k-1)) + temp*&
&             volid
            qs = temp*voli
          end if
! compute the normal velocity, where the normal direction
! is taken as the average of faces k and k-1.
          temp = sk(i, j, k, 1) + sk(i, j, k-1, 1)
          xad = voli*(skd(i, j, k, 1)+skd(i, j, k-1, 1)) + temp*volid
          xa = temp*voli
          temp = sk(i, j, k, 2) + sk(i, j, k-1, 2)
          yad = voli*(skd(i, j, k, 2)+skd(i, j, k-1, 2)) + temp*volid
          ya = temp*voli
          temp = sk(i, j, k, 3) + sk(i, j, k-1, 3)
          zad = voli*(skd(i, j, k, 3)+skd(i, j, k-1, 3)) + temp*volid
          za = temp*voli
          temp = w(i, j, k, ivx)
          temp0 = w(i, j, k, ivy)
          temp1 = w(i, j, k, ivz)
          uud = temp*xad + xa*wd(i, j, k, ivx) + temp0*yad + ya*wd(i, j&
&           , k, ivy) + temp1*zad + za*wd(i, j, k, ivz) - qsd
          uu = xa*temp + ya*temp0 + za*temp1 - 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
! 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.
                dwtm1d = wd(i, j, k-1, jj) - wd(i, j, k-2, jj)
                dwtm1 = w(i, j, k-1, jj) - w(i, j, k-2, jj)
                dwtd = wd(i, j, k, jj) - wd(i, j, k-1, jj)
                dwt = w(i, j, k, jj) - w(i, j, k-1, jj)
                dwtp1d = wd(i, j, k+1, jj) - wd(i, j, k, 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.
                dwtkd = dwtd
                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
                    dwtkd = dwtkd + half*dwtd
                    dwtk = dwtk + half*dwt
                  else
                    dwtkd = dwtkd + half*dwtp1d
                    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
                    dwtkd = dwtkd - half*dwtd
                    dwtk = dwtk - half*dwt
                  else
                    dwtkd = dwtkd - half*dwtm1d
                    dwtk = dwtk - half*dwtm1
                  end if
                end if
              else
! 1st order upwind scheme.
                dwtkd = wd(i, j, k, jj) - wd(i, j, k-1, jj)
                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)
              scratchd(i, j, k, idvt+ii-1) = scratchd(i, j, k, idvt+ii-1&
&               ) - dwtk*uud - uu*dwtkd
              scratch(i, j, k, idvt+ii-1) = scratch(i, j, k, idvt+ii-1) &
&               - uu*dwtk
            end do
          else
! 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.
                dwtm1d = wd(i, j, k, jj) - wd(i, j, k-1, jj)
                dwtm1 = w(i, j, k, jj) - w(i, j, k-1, jj)
                dwtd = wd(i, j, k+1, jj) - wd(i, j, k, jj)
                dwt = w(i, j, k+1, jj) - w(i, j, k, jj)
                dwtp1d = wd(i, j, k+2, jj) - wd(i, j, k+1, 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.
                dwtkd = dwtd
                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
                    dwtkd = dwtkd - half*dwtd
                    dwtk = dwtk - half*dwt
                  else
                    dwtkd = dwtkd - half*dwtp1d
                    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
                    dwtkd = dwtkd + half*dwtd
                    dwtk = dwtk + half*dwt
                  else
                    dwtkd = dwtkd + half*dwtm1d
                    dwtk = dwtk + half*dwtm1
                  end if
                end if
              else
! 1st order upwind scheme.
                dwtkd = wd(i, j, k+1, jj) - wd(i, j, k, jj)
                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.
              scratchd(i, j, k, idvt+ii-1) = scratchd(i, j, k, idvt+ii-1&
&               ) - dwtk*uud - uu*dwtkd
              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
      end do
    end do
!
!       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.
!
    qs = zero
    qsd = 0.0_8
    do k=2,kl
      do j=2,jl
        do i=2,il
! compute the grid velocity if present.
! it is taken as the average of j and j-1,
          temp1 = half/vol(i, j, k)
          volid = -(temp1*vold(i, j, k)/vol(i, j, k))
          voli = temp1
          if (addgridvelocities) then
            temp1 = sfacej(i, j, k) + sfacej(i, j-1, k)
            qsd = voli*(sfacejd(i, j, k)+sfacejd(i, j-1, k)) + temp1*&
&             volid
            qs = temp1*voli
          end if
! compute the normal velocity, where the normal direction
! is taken as the average of faces j and j-1.
          temp1 = sj(i, j, k, 1) + sj(i, j-1, k, 1)
          xad = voli*(sjd(i, j, k, 1)+sjd(i, j-1, k, 1)) + temp1*volid
          xa = temp1*voli
          temp1 = sj(i, j, k, 2) + sj(i, j-1, k, 2)
          yad = voli*(sjd(i, j, k, 2)+sjd(i, j-1, k, 2)) + temp1*volid
          ya = temp1*voli
          temp1 = sj(i, j, k, 3) + sj(i, j-1, k, 3)
          zad = voli*(sjd(i, j, k, 3)+sjd(i, j-1, k, 3)) + temp1*volid
          za = temp1*voli
          temp1 = w(i, j, k, ivx)
          temp0 = w(i, j, k, ivy)
          temp = w(i, j, k, ivz)
          uud = temp1*xad + xa*wd(i, j, k, ivx) + temp0*yad + ya*wd(i, j&
&           , k, ivy) + temp*zad + za*wd(i, j, k, ivz) - qsd
          uu = xa*temp1 + ya*temp0 + za*temp - qs
! determine the situation we are having here, i.e. positive
! or negative normal velocity.
          if (uu .gt. zero) then
! 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.
                dwtm1d = wd(i, j-1, k, jj) - wd(i, j-2, k, jj)
                dwtm1 = w(i, j-1, k, jj) - w(i, j-2, k, jj)
                dwtd = wd(i, j, k, jj) - wd(i, j-1, k, jj)
                dwt = w(i, j, k, jj) - w(i, j-1, k, jj)
                dwtp1d = wd(i, j+1, k, jj) - wd(i, j, 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.
                dwtjd = dwtd
                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
                    dwtjd = dwtjd + half*dwtd
                    dwtj = dwtj + half*dwt
                  else
                    dwtjd = dwtjd + half*dwtp1d
                    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
                    dwtjd = dwtjd - half*dwtd
                    dwtj = dwtj - half*dwt
                  else
                    dwtjd = dwtjd - half*dwtm1d
                    dwtj = dwtj - half*dwtm1
                  end if
                end if
              else
! 1st order upwind scheme.
                dwtjd = wd(i, j, k, jj) - wd(i, j-1, k, jj)
                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.
              scratchd(i, j, k, idvt+ii-1) = scratchd(i, j, k, idvt+ii-1&
&               ) - dwtj*uud - uu*dwtjd
              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
! 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.
                dwtm1d = wd(i, j, k, jj) - wd(i, j-1, k, jj)
                dwtm1 = w(i, j, k, jj) - w(i, j-1, k, jj)
                dwtd = wd(i, j+1, k, jj) - wd(i, j, k, jj)
                dwt = w(i, j+1, k, jj) - w(i, j, k, jj)
                dwtp1d = wd(i, j+2, k, jj) - wd(i, j+1, 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.
                dwtjd = dwtd
                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
                    dwtjd = dwtjd - half*dwtd
                    dwtj = dwtj - half*dwt
                  else
                    dwtjd = dwtjd - half*dwtp1d
                    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
                    dwtjd = dwtjd + half*dwtd
                    dwtj = dwtj + half*dwt
                  else
                    dwtjd = dwtjd + half*dwtm1d
                    dwtj = dwtj + half*dwtm1
                  end if
                end if
              else
! 1st order upwind scheme.
                dwtjd = wd(i, j+1, k, jj) - wd(i, j, k, jj)
                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.
              scratchd(i, j, k, idvt+ii-1) = scratchd(i, j, k, idvt+ii-1&
&               ) - dwtj*uud - uu*dwtjd
              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
      end do
    end do
!
!       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.
!
    qs = zero
    qsd = 0.0_8
    do k=2,kl
      do j=2,jl
        do i=2,il
! compute the grid velocity if present.
! it is taken as the average of i and i-1,
          temp1 = half/vol(i, j, k)
          volid = -(temp1*vold(i, j, k)/vol(i, j, k))
          voli = temp1
          if (addgridvelocities) then
            temp1 = sfacei(i, j, k) + sfacei(i-1, j, k)
            qsd = voli*(sfaceid(i, j, k)+sfaceid(i-1, j, k)) + temp1*&
&             volid
            qs = temp1*voli
          end if
! compute the normal velocity, where the normal direction
! is taken as the average of faces i and i-1.
          temp1 = si(i, j, k, 1) + si(i-1, j, k, 1)
          xad = voli*(sid(i, j, k, 1)+sid(i-1, j, k, 1)) + temp1*volid
          xa = temp1*voli
          temp1 = si(i, j, k, 2) + si(i-1, j, k, 2)
          yad = voli*(sid(i, j, k, 2)+sid(i-1, j, k, 2)) + temp1*volid
          ya = temp1*voli
          temp1 = si(i, j, k, 3) + si(i-1, j, k, 3)
          zad = voli*(sid(i, j, k, 3)+sid(i-1, j, k, 3)) + temp1*volid
          za = temp1*voli
          temp1 = w(i, j, k, ivx)
          temp0 = w(i, j, k, ivy)
          temp = w(i, j, k, ivz)
          uud = temp1*xad + xa*wd(i, j, k, ivx) + temp0*yad + ya*wd(i, j&
&           , k, ivy) + temp*zad + za*wd(i, j, k, ivz) - qsd
          uu = xa*temp1 + ya*temp0 + za*temp - qs
! determine the situation we are having here, i.e. positive
! or negative normal velocity.
          if (uu .gt. zero) then
! 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.
                dwtm1d = wd(i-1, j, k, jj) - wd(i-2, j, k, jj)
                dwtm1 = w(i-1, j, k, jj) - w(i-2, j, k, jj)
                dwtd = wd(i, j, k, jj) - wd(i-1, j, k, jj)
                dwt = w(i, j, k, jj) - w(i-1, j, k, jj)
                dwtp1d = wd(i+1, j, k, jj) - wd(i, 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.
                dwtid = dwtd
                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
                    dwtid = dwtid + half*dwtd
                    dwti = dwti + half*dwt
                  else
                    dwtid = dwtid + half*dwtp1d
                    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
                    dwtid = dwtid - half*dwtd
                    dwti = dwti - half*dwt
                  else
                    dwtid = dwtid - half*dwtm1d
                    dwti = dwti - half*dwtm1
                  end if
                end if
              else
! 1st order upwind scheme.
                dwtid = wd(i, j, k, jj) - wd(i-1, j, k, jj)
                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.
              scratchd(i, j, k, idvt+ii-1) = scratchd(i, j, k, idvt+ii-1&
&               ) - dwti*uud - uu*dwtid
              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
! 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.
                dwtm1d = wd(i, j, k, jj) - wd(i-1, j, k, jj)
                dwtm1 = w(i, j, k, jj) - w(i-1, j, k, jj)
                dwtd = wd(i+1, j, k, jj) - wd(i, j, k, jj)
                dwt = w(i+1, j, k, jj) - w(i, j, k, jj)
                dwtp1d = wd(i+2, j, k, jj) - wd(i+1, 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.
                dwtid = dwtd
                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
                    dwtid = dwtid - half*dwtd
                    dwti = dwti - half*dwt
                  else
                    dwtid = dwtid - half*dwtp1d
                    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
                    dwtid = dwtid + half*dwtd
                    dwti = dwti + half*dwt
                  else
                    dwtid = dwtid + half*dwtm1d
                    dwti = dwti + half*dwtm1
                  end if
                end if
              else
! 1st order upwind scheme.
                dwtid = wd(i+1, j, k, jj) - wd(i, j, k, jj)
                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.
              scratchd(i, j, k, idvt+ii-1) = scratchd(i, j, k, idvt+ii-1&
&               ) - dwti*uud - uu*dwtid
              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
      end do
    end do
  end subroutine turbadvection_d
 
  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 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
!
!       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 k=2,kl
      do j=2,jl
        do i=2,il
! 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
      end do
    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
    do k=2,kl
      do j=2,jl
        do i=2,il
! 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
      end do
    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
    do k=2,kl
      do j=2,jl
        do i=2,il
! 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
      end do
    end do
!$ad checkpoint-end
 
  end subroutine turbadvection
! ----------------------------------------------------------------------
!                                                                      |
!                    no tapenade routine below this line               |
!                                                                      |
! ----------------------------------------------------------------------
 
end module turbutils_d