sa_b.f90

Go to the documentation of this file.

!        generated by tapenade     (inria, ecuador team)
!  tapenade 3.16 (develop) - 22 aug 2023 15:51
!
! this module contains the source code related to the sa turbulence
! model. it is slightly more modularized than the original which makes
! performing reverse mode ad simplier.
module sa_b
  use constants
  implicit none
  real(kind=realtype) :: cv13, kar2inv, cw36, cb3inv
  real(kind=realtype), dimension(:, :, :), allocatable :: qq
  real(kind=realtype), dimension(:, :, :), pointer :: ddw, ww, ddvt
  real(kind=realtype), dimension(:, :), pointer :: rrlv
  real(kind=realtype), dimension(:, :), pointer :: dd2wall
 
contains
!  differentiation of sasource in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: timeref *w *rlv *scratch *vol
!                *si *sj *sk
!   with respect to varying inputs: timeref *w *rlv *scratch *vol
!                *d2wall *si *sj *sk
!   rw status of diff variables: timeref:incr *w:incr *rlv:incr
!                *scratch:in-out *vol:incr *d2wall:out *si:incr
!                *sj:incr *sk:incr
!   plus diff mem management of: w:in rlv:in scratch:in vol:in
!                d2wall:in si:in sj:in sk:in
  subroutine sasource_b()
!
!  source terms.
!  determine the source term and its derivative w.r.t. nutilde
!  for all internal cells of the block.
!  remember that the sa field variable nutilde = w(i,j,k,itu1)
    use blockpointers
    use constants
    use paramturb
    use section
    use inputphysics
    use inputdiscretization, only : approxsa
    use flowvarrefstate
    implicit none
! local parameters
    real(kind=realtype), parameter :: f23=two*third
! local variables.
    integer(kind=inttype) :: i, j, k, nn, ii
    real(kind=realtype) :: fv1, fv2, ft2
    real(kind=realtype) :: fv1d, fv2d, ft2d
    real(kind=realtype) :: ss, sst, nu, dist2inv, chi, chi2, chi3
    real(kind=realtype) :: ssd, sstd, nud, dist2invd, chid, chi2d, chi3d
    real(kind=realtype) :: rr, gg, gg6, termfw, fwsa, term1, term2
    real(kind=realtype) :: rrd, ggd, gg6d, termfwd, fwsad, term1d, &
&   term2d
    real(kind=realtype) :: dfv1, dfv2, dft2, drr, dgg, dfw
    real(kind=realtype) :: uux, uuy, uuz, vvx, vvy, vvz, wwx, wwy, wwz
    real(kind=realtype) :: uuxd, uuyd, uuzd, vvxd, vvyd, vvzd, wwxd, &
&   wwyd, wwzd
    real(kind=realtype) :: div2, fact, sxx, syy, szz, sxy, sxz, syz
    real(kind=realtype) :: div2d, factd, sxxd, syyd, szzd, sxyd, sxzd, &
&   syzd
    real(kind=realtype) :: vortx, vorty, vortz
    real(kind=realtype) :: vortxd, vortyd, vortzd
    real(kind=realtype) :: omegax, omegay, omegaz
    real(kind=realtype) :: omegaxd, omegayd, omegazd
    real(kind=realtype) :: strainmag2, strainprod, vortprod
    real(kind=realtype) :: strainmag2d, strainprodd, vortprodd
    real(kind=realtype), parameter :: xminn=1.e-10_realtype
    intrinsic mod
    intrinsic sqrt
    intrinsic exp
    intrinsic min
    intrinsic max
    real(kind=realtype) :: y1
    real(kind=realtype) :: y1d
    real(kind=realtype) :: min1
    real(kind=realtype) :: min1d
    real(kind=realtype) :: temp
    real(kind=realtype) :: tempd
    real(kind=realtype) :: temp0
    real(kind=realtype) :: tempd0
    integer :: branch
! set model constants
    cv13 = rsacv1**3
    kar2inv = one/rsak**2
    cw36 = rsacw3**6
! 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)
! create switches to production term depending on the variable that
! should be used
    if (turbprod .eq. katolaunder) then
      stop
    else
      if (associated(d2walld)) d2walld = 0.0_8
      omegaxd = 0.0_8
      omegayd = 0.0_8
      omegazd = 0.0_8
      strainmag2d = 0.0_8
      ssd = 0.0_8
!$bwd-of ii-loop 
      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)
        if (turbprod .eq. strain) then
          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
! compute strain production term
          strainmag2 = two*(sxy**2+sxz**2+syz**2) + sxx**2 + syy**2 + &
&           szz**2
          strainprod = two*strainmag2 - div2
          ss = sqrt(strainprod)
          call pushcontrol2b(0)
        else if (turbprod .eq. vorticity) then
! compute the three components of the vorticity vector.
! substract the part coming from the rotating frame.
          vortx = two*fact*(wwy-vvz) - two*omegax
          vorty = two*fact*(uuz-wwx) - two*omegay
          vortz = two*fact*(vvx-uuy) - two*omegaz
! compute the vorticity production term
          vortprod = vortx**2 + vorty**2 + vortz**2
! first take the square root of the production term to
! obtain the correct production term for spalart-allmaras.
! we do this to avoid if statements.
          ss = sqrt(vortprod)
          call pushcontrol2b(1)
        else
          call pushcontrol2b(2)
        end if
! compute the laminar kinematic viscosity, the inverse of
! wall distance squared, the ratio chi (ratio of nutilde
! and nu) and the functions fv1 and fv2. the latter corrects
! the production term near a viscous wall.
        nu = rlv(i, j, k)/w(i, j, k, irho)
        dist2inv = one/d2wall(i, j, k)**2
        chi = w(i, j, k, itu1)/nu
        chi2 = chi*chi
        chi3 = chi*chi2
        fv1 = chi3/(chi3+cv13)
        fv2 = one - chi/(one+chi*fv1)
! the function ft2, which is designed to keep a laminar
! solution laminar. when running in fully turbulent mode
! this function should be set to 0.0.
        if (useft2sa) then
          ft2 = rsact3*exp(-(rsact4*chi2))
          call pushcontrol1b(0)
        else
          ft2 = zero
          call pushcontrol1b(1)
        end if
! correct the production term to account for the influence
! of the wall.
        sst = ss + w(i, j, k, itu1)*fv2*kar2inv*dist2inv
! add rotation term (userotationsa defined in inputparams.f90)
        if (userotationsa) then
          y1 = sqrt(two*strainmag2)
          if (zero .gt. y1) then
            min1 = y1
            call pushcontrol1b(0)
          else
            min1 = zero
            call pushcontrol1b(1)
          end if
          sst = sst + rsacrot*min1
          call pushcontrol1b(1)
        else
          call pushcontrol1b(0)
        end if
        if (sst .lt. xminn) then
          sst = xminn
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
          sst = sst
        end if
! compute the function fw. the argument rr is cut off at 10
! to avoid numerical problems. this is ok, because the
! asymptotical value of fw is then already reached.
        rr = w(i, j, k, itu1)*kar2inv*dist2inv/sst
        if (rr .gt. 10.0_realtype) then
          rr = 10.0_realtype
          call pushcontrol1b(0)
        else
          call pushcontrol1b(1)
          rr = rr
        end if
        gg = rr + rsacw2*(rr**6-rr)
        gg6 = gg**6
        termfw = ((one+cw36)/(gg6+cw36))**sixth
        fwsa = gg*termfw
! compute the source term; some terms are saved for the
! linearization. the source term is stored in dvt.
        if (approxsa) then
          call pushcontrol1b(0)
          term1 = zero
        else
          term1 = rsacb1*(one-ft2)*ss
          call pushcontrol1b(1)
        end if
        term2 = dist2inv*(kar2inv*rsacb1*((one-ft2)*fv2+ft2)-rsacw1*fwsa&
&         )
        temp = w(i, j, k, itu1)
        tempd0 = w(i, j, k, itu1)*scratchd(i, j, k, idvt)
        wd(i, j, k, itu1) = wd(i, j, k, itu1) + (term1+term2*temp)*&
&         scratchd(i, j, k, idvt) + term2*tempd0
        scratchd(i, j, k, idvt) = 0.0_8
        term1d = tempd0
        term2d = temp*tempd0
        dist2invd = (kar2inv*rsacb1*((one-ft2)*fv2+ft2)-rsacw1*fwsa)*&
&         term2d
        tempd0 = kar2inv*rsacb1*dist2inv*term2d
        fwsad = -(rsacw1*dist2inv*term2d)
        ft2d = (1.0-fv2)*tempd0
        fv2d = (one-ft2)*tempd0
        call popcontrol1b(branch)
        if (branch .ne. 0) then
          ft2d = ft2d - ss*rsacb1*term1d
          ssd = ssd + (one-ft2)*rsacb1*term1d
        end if
        termfwd = gg*fwsad
        temp0 = (one+cw36)/(cw36+gg6)
        if (temp0 .le. 0.0_8 .and. (sixth .eq. 0.0_8 .or. sixth .ne. int&
&           (sixth))) then
          gg6d = 0.0_8
        else
          gg6d = -(temp0*sixth*temp0**(sixth-1)*termfwd/(cw36+gg6))
        end if
        ggd = termfw*fwsad + 6*gg**5*gg6d
        rrd = (6*rr**5*rsacw2-rsacw2+1.0)*ggd
        call popcontrol1b(branch)
        if (branch .eq. 0) rrd = 0.0_8
        wd(i, j, k, itu1) = wd(i, j, k, itu1) + dist2inv*kar2inv*rrd/sst
        tempd0 = w(i, j, k, itu1)*kar2inv*rrd/sst
        dist2invd = dist2invd + tempd0
        sstd = -(dist2inv*tempd0/sst)
        call popcontrol1b(branch)
        if (branch .eq. 0) sstd = 0.0_8
        call popcontrol1b(branch)
        if (branch .ne. 0) then
          min1d = rsacrot*sstd
          call popcontrol1b(branch)
          if (branch .eq. 0) then
            y1d = min1d
          else
            y1d = 0.0_8
          end if
          if (.not.two*strainmag2 .eq. 0.0_8) strainmag2d = strainmag2d &
&             + two*y1d/(2.0*sqrt(two*strainmag2))
        end if
        ssd = ssd + sstd
        wd(i, j, k, itu1) = wd(i, j, k, itu1) + fv2*dist2inv*kar2inv*&
&         sstd
        tempd0 = w(i, j, k, itu1)*kar2inv*sstd
        fv2d = fv2d + dist2inv*tempd0
        dist2invd = dist2invd + fv2*tempd0
        call popcontrol1b(branch)
        if (branch .eq. 0) then
          chi2d = -(rsact4*exp(-(rsact4*chi2))*rsact3*ft2d)
        else
          chi2d = 0.0_8
        end if
        tempd = -(fv2d/(one+chi*fv1))
        chid = tempd
        tempd0 = -(chi*tempd/(one+chi*fv1))
        fv1d = chi*tempd0
        tempd = fv1d/(cv13+chi3)
        chi3d = (1.0-chi3/(cv13+chi3))*tempd
        chi2d = chi2d + chi*chi3d
        chid = chid + fv1*tempd0 + chi2*chi3d + 2*chi*chi2d
        wd(i, j, k, itu1) = wd(i, j, k, itu1) + chid/nu
        nud = -(w(i, j, k, itu1)*chid/nu**2)
        temp = d2wall(i, j, k)
        d2walld(i, j, k) = d2walld(i, j, k) - 2*one*dist2invd/temp**3
        temp = w(i, j, k, irho)
        rlvd(i, j, k) = rlvd(i, j, k) + nud/temp
        wd(i, j, k, irho) = wd(i, j, k, irho) - rlv(i, j, k)*nud/temp**2
        call popcontrol2b(branch)
        if (branch .eq. 0) then
          if (strainprod .eq. 0.0_8) then
            strainprodd = 0.0_8
          else
            strainprodd = ssd/(2.0*sqrt(strainprod))
          end if
          strainmag2d = strainmag2d + two*strainprodd
          div2d = -strainprodd
          tempd = two*strainmag2d
          sxyd = 2*sxy*tempd
          sxzd = 2*sxz*tempd
          syzd = 2*syz*tempd
          tempd = 2*(sxx+syy+szz)*f23*div2d
          sxxd = 2*sxx*strainmag2d + tempd
          syyd = 2*syy*strainmag2d + tempd
          szzd = 2*szz*strainmag2d + tempd
          factd = (vvz+wwy)*syzd + (uuz+wwx)*sxzd + (uuy+vvx)*sxyd + wwz&
&           *two*szzd + vvy*two*syyd + uux*two*sxxd
          vvzd = fact*syzd
          wwyd = fact*syzd
          uuzd = fact*sxzd
          wwxd = fact*sxzd
          uuyd = fact*sxyd
          vvxd = fact*sxyd
          wwzd = fact*two*szzd
          vvyd = fact*two*syyd
          uuxd = fact*two*sxxd
          strainmag2d = 0.0_8
          ssd = 0.0_8
        else
          if (branch .eq. 1) then
            if (vortprod .eq. 0.0_8) then
              vortprodd = 0.0_8
            else
              vortprodd = ssd/(2.0*sqrt(vortprod))
            end if
            vortxd = 2*vortx*vortprodd
            vortyd = 2*vorty*vortprodd
            vortzd = 2*vortz*vortprodd
            tempd = two*vortzd
            omegazd = omegazd - two*vortzd
            factd = (vvx-uuy)*tempd
            vvxd = fact*tempd
            uuyd = -(fact*tempd)
            tempd = two*vortyd
            omegayd = omegayd - two*vortyd
            factd = factd + (uuz-wwx)*tempd
            uuzd = fact*tempd
            wwxd = -(fact*tempd)
            tempd = two*vortxd
            omegaxd = omegaxd - two*vortxd
            factd = factd + (wwy-vvz)*tempd
            wwyd = fact*tempd
            vvzd = -(fact*tempd)
            ssd = 0.0_8
          else
            wwxd = 0.0_8
            wwyd = 0.0_8
            vvxd = 0.0_8
            vvzd = 0.0_8
            factd = 0.0_8
            uuyd = 0.0_8
            uuzd = 0.0_8
          end if
          wwzd = 0.0_8
          vvyd = 0.0_8
          uuxd = 0.0_8
        end if
        wd(i, j, k-1, ivz) = wd(i, j, k-1, ivz) - sk(i, j, k-1, 3)*wwzd &
&         - sk(i, j, k-1, 2)*wwyd - sk(i, j, k-1, 1)*wwxd
        wd(i, j-1, k, ivz) = wd(i, j-1, k, ivz) - sj(i, j-1, k, 3)*wwzd &
&         - sj(i, j-1, k, 2)*wwyd - sj(i, j-1, k, 1)*wwxd
        wd(i, j, k+1, ivz) = wd(i, j, k+1, ivz) + sk(i, j, k, 3)*wwzd + &
&         sk(i, j, k, 2)*wwyd + sk(i, j, k, 1)*wwxd
        wd(i, j+1, k, ivz) = wd(i, j+1, k, ivz) + sj(i, j, k, 3)*wwzd + &
&         sj(i, j, k, 2)*wwyd + sj(i, j, k, 1)*wwxd
        wd(i-1, j, k, ivz) = wd(i-1, j, k, ivz) - si(i-1, j, k, 3)*wwzd &
&         - si(i-1, j, k, 2)*wwyd - si(i-1, j, k, 1)*wwxd
        vold(i, j, k) = vold(i, j, k) - fourth*factd/vol(i, j, k)**2
        wd(i+1, j, k, ivz) = wd(i+1, j, k, ivz) + si(i, j, k, 3)*wwzd + &
&         si(i, j, k, 2)*wwyd + si(i, j, k, 1)*wwxd
        sid(i, j, k, 3) = sid(i, j, k, 3) + w(i+1, j, k, ivz)*wwzd
        sid(i-1, j, k, 3) = sid(i-1, j, k, 3) - w(i-1, j, k, ivz)*wwzd
        sjd(i, j, k, 3) = sjd(i, j, k, 3) + w(i, j+1, k, ivz)*wwzd
        skd(i, j, k, 3) = skd(i, j, k, 3) + w(i, j, k+1, ivz)*wwzd
        sjd(i, j-1, k, 3) = sjd(i, j-1, k, 3) - w(i, j-1, k, ivz)*wwzd
        skd(i, j, k-1, 3) = skd(i, j, k-1, 3) - w(i, j, k-1, ivz)*wwzd
        sid(i, j, k, 2) = sid(i, j, k, 2) + w(i+1, j, k, ivz)*wwyd
        sid(i-1, j, k, 2) = sid(i-1, j, k, 2) - w(i-1, j, k, ivz)*wwyd
        sjd(i, j, k, 2) = sjd(i, j, k, 2) + w(i, j+1, k, ivz)*wwyd
        skd(i, j, k, 2) = skd(i, j, k, 2) + w(i, j, k+1, ivz)*wwyd
        sjd(i, j-1, k, 2) = sjd(i, j-1, k, 2) - w(i, j-1, k, ivz)*wwyd
        skd(i, j, k-1, 2) = skd(i, j, k-1, 2) - w(i, j, k-1, ivz)*wwyd
        sid(i, j, k, 1) = sid(i, j, k, 1) + w(i+1, j, k, ivz)*wwxd
        sid(i-1, j, k, 1) = sid(i-1, j, k, 1) - w(i-1, j, k, ivz)*wwxd
        sjd(i, j, k, 1) = sjd(i, j, k, 1) + w(i, j+1, k, ivz)*wwxd
        skd(i, j, k, 1) = skd(i, j, k, 1) + w(i, j, k+1, ivz)*wwxd
        sjd(i, j-1, k, 1) = sjd(i, j-1, k, 1) - w(i, j-1, k, ivz)*wwxd
        skd(i, j, k-1, 1) = skd(i, j, k-1, 1) - w(i, j, k-1, ivz)*wwxd
        sid(i, j, k, 3) = sid(i, j, k, 3) + w(i+1, j, k, ivy)*vvzd
        sid(i-1, j, k, 3) = sid(i-1, j, k, 3) - w(i-1, j, k, ivy)*vvzd
        sjd(i, j, k, 3) = sjd(i, j, k, 3) + w(i, j+1, k, ivy)*vvzd
        skd(i, j, k, 3) = skd(i, j, k, 3) + w(i, j, k+1, ivy)*vvzd
        sjd(i, j-1, k, 3) = sjd(i, j-1, k, 3) - w(i, j-1, k, ivy)*vvzd
        wd(i, j, k-1, ivy) = wd(i, j, k-1, ivy) - sk(i, j, k-1, 3)*vvzd &
&         - sk(i, j, k-1, 2)*vvyd - sk(i, j, k-1, 1)*vvxd
        wd(i, j-1, k, ivy) = wd(i, j-1, k, ivy) - sj(i, j-1, k, 3)*vvzd &
&         - sj(i, j-1, k, 2)*vvyd - sj(i, j-1, k, 1)*vvxd
        wd(i, j, k+1, ivy) = wd(i, j, k+1, ivy) + sk(i, j, k, 3)*vvzd + &
&         sk(i, j, k, 2)*vvyd + sk(i, j, k, 1)*vvxd
        wd(i, j+1, k, ivy) = wd(i, j+1, k, ivy) + sj(i, j, k, 3)*vvzd + &
&         sj(i, j, k, 2)*vvyd + sj(i, j, k, 1)*vvxd
        wd(i-1, j, k, ivy) = wd(i-1, j, k, ivy) - si(i-1, j, k, 3)*vvzd &
&         - si(i-1, j, k, 2)*vvyd - si(i-1, j, k, 1)*vvxd
        wd(i+1, j, k, ivy) = wd(i+1, j, k, ivy) + si(i, j, k, 3)*vvzd + &
&         si(i, j, k, 2)*vvyd + si(i, j, k, 1)*vvxd
        skd(i, j, k-1, 3) = skd(i, j, k-1, 3) - w(i, j, k-1, ivy)*vvzd
        sid(i, j, k, 2) = sid(i, j, k, 2) + w(i+1, j, k, ivy)*vvyd
        sid(i-1, j, k, 2) = sid(i-1, j, k, 2) - w(i-1, j, k, ivy)*vvyd
        sjd(i, j, k, 2) = sjd(i, j, k, 2) + w(i, j+1, k, ivy)*vvyd
        skd(i, j, k, 2) = skd(i, j, k, 2) + w(i, j, k+1, ivy)*vvyd
        sjd(i, j-1, k, 2) = sjd(i, j-1, k, 2) - w(i, j-1, k, ivy)*vvyd
        skd(i, j, k-1, 2) = skd(i, j, k-1, 2) - w(i, j, k-1, ivy)*vvyd
        sid(i, j, k, 1) = sid(i, j, k, 1) + w(i+1, j, k, ivy)*vvxd
        sid(i-1, j, k, 1) = sid(i-1, j, k, 1) - w(i-1, j, k, ivy)*vvxd
        sjd(i, j, k, 1) = sjd(i, j, k, 1) + w(i, j+1, k, ivy)*vvxd
        skd(i, j, k, 1) = skd(i, j, k, 1) + w(i, j, k+1, ivy)*vvxd
        sjd(i, j-1, k, 1) = sjd(i, j-1, k, 1) - w(i, j-1, k, ivy)*vvxd
        skd(i, j, k-1, 1) = skd(i, j, k-1, 1) - w(i, j, k-1, ivy)*vvxd
        sid(i, j, k, 3) = sid(i, j, k, 3) + w(i+1, j, k, ivx)*uuzd
        sid(i-1, j, k, 3) = sid(i-1, j, k, 3) - w(i-1, j, k, ivx)*uuzd
        sjd(i, j, k, 3) = sjd(i, j, k, 3) + w(i, j+1, k, ivx)*uuzd
        skd(i, j, k, 3) = skd(i, j, k, 3) + w(i, j, k+1, ivx)*uuzd
        sjd(i, j-1, k, 3) = sjd(i, j-1, k, 3) - w(i, j-1, k, ivx)*uuzd
        wd(i, j, k-1, ivx) = wd(i, j, k-1, ivx) - sk(i, j, k-1, 3)*uuzd &
&         - sk(i, j, k-1, 2)*uuyd - sk(i, j, k-1, 1)*uuxd
        wd(i, j-1, k, ivx) = wd(i, j-1, k, ivx) - sj(i, j-1, k, 3)*uuzd &
&         - sj(i, j-1, k, 2)*uuyd - sj(i, j-1, k, 1)*uuxd
        wd(i, j, k+1, ivx) = wd(i, j, k+1, ivx) + sk(i, j, k, 3)*uuzd + &
&         sk(i, j, k, 2)*uuyd + sk(i, j, k, 1)*uuxd
        wd(i, j+1, k, ivx) = wd(i, j+1, k, ivx) + sj(i, j, k, 3)*uuzd + &
&         sj(i, j, k, 2)*uuyd + sj(i, j, k, 1)*uuxd
        wd(i-1, j, k, ivx) = wd(i-1, j, k, ivx) - si(i-1, j, k, 3)*uuzd &
&         - si(i-1, j, k, 2)*uuyd - si(i-1, j, k, 1)*uuxd
        wd(i+1, j, k, ivx) = wd(i+1, j, k, ivx) + si(i, j, k, 3)*uuzd + &
&         si(i, j, k, 2)*uuyd + si(i, j, k, 1)*uuxd
        skd(i, j, k-1, 3) = skd(i, j, k-1, 3) - w(i, j, k-1, ivx)*uuzd
        sid(i, j, k, 2) = sid(i, j, k, 2) + w(i+1, j, k, ivx)*uuyd
        sid(i-1, j, k, 2) = sid(i-1, j, k, 2) - w(i-1, j, k, ivx)*uuyd
        sjd(i, j, k, 2) = sjd(i, j, k, 2) + w(i, j+1, k, ivx)*uuyd
        skd(i, j, k, 2) = skd(i, j, k, 2) + w(i, j, k+1, ivx)*uuyd
        sjd(i, j-1, k, 2) = sjd(i, j-1, k, 2) - w(i, j-1, k, ivx)*uuyd
        skd(i, j, k-1, 2) = skd(i, j, k-1, 2) - w(i, j, k-1, ivx)*uuyd
        sid(i, j, k, 1) = sid(i, j, k, 1) + w(i+1, j, k, ivx)*uuxd
        sid(i-1, j, k, 1) = sid(i-1, j, k, 1) - w(i-1, j, k, ivx)*uuxd
        sjd(i, j, k, 1) = sjd(i, j, k, 1) + w(i, j+1, k, ivx)*uuxd
        skd(i, j, k, 1) = skd(i, j, k, 1) + w(i, j, k+1, ivx)*uuxd
        sjd(i, j-1, k, 1) = sjd(i, j-1, k, 1) - w(i, j-1, k, ivx)*uuxd
        skd(i, j, k-1, 1) = skd(i, j, k-1, 1) - w(i, j, k-1, ivx)*uuxd
      end do
      timerefd = timerefd + sections(sectionid)%rotrate(3)*omegazd + &
&       sections(sectionid)%rotrate(2)*omegayd + sections(sectionid)%&
&       rotrate(1)*omegaxd
    end if
  end subroutine sasource_b
 
  subroutine sasource()
!
!  source terms.
!  determine the source term and its derivative w.r.t. nutilde
!  for all internal cells of the block.
!  remember that the sa field variable nutilde = w(i,j,k,itu1)
    use blockpointers
    use constants
    use paramturb
    use section
    use inputphysics
    use inputdiscretization, only : approxsa
    use flowvarrefstate
    implicit none
! local parameters
    real(kind=realtype), parameter :: f23=two*third
! local variables.
    integer(kind=inttype) :: i, j, k, nn, ii
    real(kind=realtype) :: fv1, fv2, ft2
    real(kind=realtype) :: ss, sst, nu, dist2inv, chi, chi2, chi3
    real(kind=realtype) :: rr, gg, gg6, termfw, fwsa, term1, term2
    real(kind=realtype) :: dfv1, dfv2, dft2, drr, dgg, dfw
    real(kind=realtype) :: uux, uuy, uuz, vvx, vvy, vvz, wwx, wwy, wwz
    real(kind=realtype) :: div2, fact, sxx, syy, szz, sxy, sxz, syz
    real(kind=realtype) :: vortx, vorty, vortz
    real(kind=realtype) :: omegax, omegay, omegaz
    real(kind=realtype) :: strainmag2, strainprod, vortprod
    real(kind=realtype), parameter :: xminn=1.e-10_realtype
    intrinsic mod
    intrinsic sqrt
    intrinsic exp
    intrinsic min
    intrinsic max
    real(kind=realtype) :: y1
    real(kind=realtype) :: min1
! set model constants
    cv13 = rsacv1**3
    kar2inv = one/rsak**2
    cw36 = rsacw3**6
    cb3inv = one/rsacb3
! 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)
! create switches to production term depending on the variable that
! should be used
    if (turbprod .eq. katolaunder) then
      print*, 'katolaunder production term not supported for sa'
      stop
    else
!$ad ii-loop
      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)
        if (turbprod .eq. strain) then
          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
! compute strain production term
          strainmag2 = two*(sxy**2+sxz**2+syz**2) + sxx**2 + syy**2 + &
&           szz**2
          strainprod = two*strainmag2 - div2
          ss = sqrt(strainprod)
        else if (turbprod .eq. vorticity) then
! compute the three components of the vorticity vector.
! substract the part coming from the rotating frame.
          vortx = two*fact*(wwy-vvz) - two*omegax
          vorty = two*fact*(uuz-wwx) - two*omegay
          vortz = two*fact*(vvx-uuy) - two*omegaz
! compute the vorticity production term
          vortprod = vortx**2 + vorty**2 + vortz**2
! first take the square root of the production term to
! obtain the correct production term for spalart-allmaras.
! we do this to avoid if statements.
          ss = sqrt(vortprod)
        end if
! compute the laminar kinematic viscosity, the inverse of
! wall distance squared, the ratio chi (ratio of nutilde
! and nu) and the functions fv1 and fv2. the latter corrects
! the production term near a viscous wall.
        nu = rlv(i, j, k)/w(i, j, k, irho)
        dist2inv = one/d2wall(i, j, k)**2
        chi = w(i, j, k, itu1)/nu
        chi2 = chi*chi
        chi3 = chi*chi2
        fv1 = chi3/(chi3+cv13)
        fv2 = one - chi/(one+chi*fv1)
! the function ft2, which is designed to keep a laminar
! solution laminar. when running in fully turbulent mode
! this function should be set to 0.0.
        if (useft2sa) then
          ft2 = rsact3*exp(-(rsact4*chi2))
        else
          ft2 = zero
        end if
! correct the production term to account for the influence
! of the wall.
        sst = ss + w(i, j, k, itu1)*fv2*kar2inv*dist2inv
! add rotation term (userotationsa defined in inputparams.f90)
        if (userotationsa) then
          y1 = sqrt(two*strainmag2)
          if (zero .gt. y1) then
            min1 = y1
          else
            min1 = zero
          end if
          sst = sst + rsacrot*min1
        end if
        if (sst .lt. xminn) then
          sst = xminn
        else
          sst = sst
        end if
! compute the function fw. the argument rr is cut off at 10
! to avoid numerical problems. this is ok, because the
! asymptotical value of fw is then already reached.
        rr = w(i, j, k, itu1)*kar2inv*dist2inv/sst
        if (rr .gt. 10.0_realtype) then
          rr = 10.0_realtype
        else
          rr = rr
        end if
        gg = rr + rsacw2*(rr**6-rr)
        gg6 = gg**6
        termfw = ((one+cw36)/(gg6+cw36))**sixth
        fwsa = gg*termfw
! compute the source term; some terms are saved for the
! linearization. the source term is stored in dvt.
        if (approxsa) then
          term1 = zero
        else
          term1 = rsacb1*(one-ft2)*ss
        end if
        term2 = dist2inv*(kar2inv*rsacb1*((one-ft2)*fv2+ft2)-rsacw1*fwsa&
&         )
        scratch(i, j, k, idvt) = (term1+term2*w(i, j, k, itu1))*w(i, j, &
&         k, itu1)
      end do
    end if
  end subroutine sasource
 
!  differentiation of saviscous in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *w *rlv *scratch *vol *si *sj
!                *sk
!   with respect to varying inputs: *w *rlv *scratch *vol *si *sj
!                *sk
!   rw status of diff variables: *w:incr *rlv:incr *scratch:in-out
!                *vol:incr *si:incr *sj:incr *sk:incr
!   plus diff mem management of: w:in rlv:in scratch:in vol:in
!                si:in sj:in sk:in
  subroutine saviscous_b()
!
!  viscous term.
!  determine the viscous contribution to the residual
!  for all internal cells of the block.
    use blockpointers
    use paramturb
    implicit none
! local variables.
    integer(kind=inttype) :: i, j, k, nn, ii
    real(kind=realtype) :: nu
    real(kind=realtype) :: nud
    real(kind=realtype) :: fv1, fv2, ft2
    real(kind=realtype) :: voli, volmi, volpi, xm, ym, zm, xp, yp, zp
    real(kind=realtype) :: volid, volmid, volpid, xmd, ymd, zmd, xpd, &
&   ypd, zpd
    real(kind=realtype) :: xa, ya, za, ttm, ttp, cnud, cam, cap
    real(kind=realtype) :: xad, yad, zad, ttmd, ttpd, cnudd, camd, capd
    real(kind=realtype) :: nutm, nutp, num, nup, cdm, cdp
    real(kind=realtype) :: nutmd, nutpd, numd, nupd, cdmd, cdpd
    real(kind=realtype) :: c1m, c1p, c10, b1, c1, d1, qs
    real(kind=realtype) :: c1md, c1pd, c10d
    intrinsic mod
    intrinsic max
    real(kind=realtype) :: temp
    real(kind=realtype) :: tempd
    real(kind=realtype) :: temp0
    real(kind=realtype) :: tempd0
    integer :: branch
! set model constants
    cb3inv = one/rsacb3
!$bwd-of ii-loop 
    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 metrics in xi-direction, i.e. along the
! line i = constant.
      voli = one/vol(i, j, k)
      volmi = two/(vol(i, j, k)+vol(i-1, j, k))
      volpi = two/(vol(i, j, k)+vol(i+1, j, k))
      xm = si(i-1, j, k, 1)*volmi
      ym = si(i-1, j, k, 2)*volmi
      zm = si(i-1, j, k, 3)*volmi
      xp = si(i, j, k, 1)*volpi
      yp = si(i, j, k, 2)*volpi
      zp = si(i, j, k, 3)*volpi
      xa = half*(si(i, j, k, 1)+si(i-1, j, k, 1))*voli
      ya = half*(si(i, j, k, 2)+si(i-1, j, k, 2))*voli
      za = half*(si(i, j, k, 3)+si(i-1, j, k, 3))*voli
      ttm = xm*xa + ym*ya + zm*za
      ttp = xp*xa + yp*ya + zp*za
! computation of the viscous terms in xi-direction; note
! that cross-derivatives are neglected, i.e. the mesh is
! assumed to be orthogonal.
! furthermore, the grad(nu)**2 has been rewritten as
! div(nu grad(nu)) - nu div(grad nu) to enhance stability.
! the second derivative in xi-direction is constructed as
! the central difference of the first order derivatives, i.e.
! d^2/dxi^2 = d/dxi (d/dxi i+1/2 - d/dxi i-1/2).
! in this way the metric can be taken into account.
! compute the diffusion coefficients multiplying the nodes
! i+1, i and i-1 in the second derivative. make sure that
! these coefficients are nonnegative.
      cnud = -(rsacb2*w(i, j, k, itu1)*cb3inv)
      cam = ttm*cnud
      cap = ttp*cnud
      nutm = half*(w(i-1, j, k, itu1)+w(i, j, k, itu1))
      nutp = half*(w(i+1, j, k, itu1)+w(i, j, k, itu1))
      nu = rlv(i, j, k)/w(i, j, k, irho)
      num = half*(rlv(i-1, j, k)/w(i-1, j, k, irho)+nu)
      nup = half*(rlv(i+1, j, k)/w(i+1, j, k, irho)+nu)
      cdm = (num+(one+rsacb2)*nutm)*ttm*cb3inv
      cdp = (nup+(one+rsacb2)*nutp)*ttp*cb3inv
      if (cdm + cam .lt. zero) then
        c1m = zero
        call pushcontrol1b(0)
      else
        c1m = cdm + cam
        call pushcontrol1b(1)
      end if
      if (cdp + cap .lt. zero) then
        c1p = zero
        call pushcontrol1b(0)
      else
        c1p = cdp + cap
        call pushcontrol1b(1)
      end if
      c10 = c1m + c1p
! update the residual for this cell and store the possible
! coefficients for the matrix in b1, c1 and d1.
      c10d = -(w(i, j, k, itu1)*scratchd(i, j, k, idvt))
      c1md = w(i-1, j, k, itu1)*scratchd(i, j, k, idvt) + c10d
      wd(i-1, j, k, itu1) = wd(i-1, j, k, itu1) + c1m*scratchd(i, j, k, &
&       idvt)
      c1pd = w(i+1, j, k, itu1)*scratchd(i, j, k, idvt) + c10d
      wd(i+1, j, k, itu1) = wd(i+1, j, k, itu1) + c1p*scratchd(i, j, k, &
&       idvt)
      wd(i, j, k, itu1) = wd(i, j, k, itu1) - c10*scratchd(i, j, k, idvt&
&       )
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        cdpd = 0.0_8
        capd = 0.0_8
      else
        cdpd = c1pd
        capd = c1pd
      end if
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        cdmd = 0.0_8
        camd = 0.0_8
      else
        cdmd = c1md
        camd = c1md
      end if
      cnudd = ttp*capd + ttm*camd
      tempd0 = cb3inv*cdpd
      nupd = ttp*tempd0
      nutpd = (one+rsacb2)*ttp*tempd0
      ttpd = (nup+(one+rsacb2)*nutp)*tempd0 + cnud*capd
      tempd0 = cb3inv*cdmd
      numd = ttm*tempd0
      nutmd = (one+rsacb2)*ttm*tempd0
      ttmd = (num+(one+rsacb2)*nutm)*tempd0 + cnud*camd
      temp0 = w(i+1, j, k, irho)
      tempd = half*nupd/temp0
      nud = half*nupd + half*numd
      rlvd(i+1, j, k) = rlvd(i+1, j, k) + tempd
      wd(i+1, j, k, irho) = wd(i+1, j, k, irho) - rlv(i+1, j, k)*tempd/&
&       temp0
      temp0 = w(i-1, j, k, irho)
      tempd = half*numd/temp0
      rlvd(i-1, j, k) = rlvd(i-1, j, k) + tempd
      wd(i-1, j, k, irho) = wd(i-1, j, k, irho) - rlv(i-1, j, k)*tempd/&
&       temp0
      temp0 = w(i, j, k, irho)
      rlvd(i, j, k) = rlvd(i, j, k) + nud/temp0
      wd(i, j, k, irho) = wd(i, j, k, irho) - rlv(i, j, k)*nud/temp0**2
      wd(i+1, j, k, itu1) = wd(i+1, j, k, itu1) + half*nutpd
      wd(i, j, k, itu1) = wd(i, j, k, itu1) + half*nutpd + half*nutmd - &
&       rsacb2*cb3inv*cnudd
      wd(i-1, j, k, itu1) = wd(i-1, j, k, itu1) + half*nutmd
      xpd = xa*ttpd
      xad = xp*ttpd + xm*ttmd
      ypd = ya*ttpd
      yad = yp*ttpd + ym*ttmd
      zpd = za*ttpd
      zad = zp*ttpd + zm*ttmd
      xmd = xa*ttmd
      ymd = ya*ttmd
      zmd = za*ttmd
      tempd0 = half*zad
      sid(i, j, k, 3) = sid(i, j, k, 3) + voli*tempd0
      sid(i-1, j, k, 3) = sid(i-1, j, k, 3) + voli*tempd0
      volid = (si(i, j, k, 3)+si(i-1, j, k, 3))*tempd0
      tempd0 = half*yad
      sid(i, j, k, 2) = sid(i, j, k, 2) + voli*tempd0
      sid(i-1, j, k, 2) = sid(i-1, j, k, 2) + voli*tempd0
      volid = volid + (si(i, j, k, 2)+si(i-1, j, k, 2))*tempd0
      tempd0 = half*xad
      sid(i, j, k, 1) = sid(i, j, k, 1) + voli*tempd0
      sid(i-1, j, k, 1) = sid(i-1, j, k, 1) + voli*tempd0
      volid = volid + (si(i, j, k, 1)+si(i-1, j, k, 1))*tempd0
      sid(i, j, k, 3) = sid(i, j, k, 3) + volpi*zpd
      volpid = si(i, j, k, 3)*zpd + si(i, j, k, 2)*ypd + si(i, j, k, 1)*&
&       xpd
      sid(i, j, k, 2) = sid(i, j, k, 2) + volpi*ypd
      sid(i, j, k, 1) = sid(i, j, k, 1) + volpi*xpd
      sid(i-1, j, k, 3) = sid(i-1, j, k, 3) + volmi*zmd
      volmid = si(i-1, j, k, 3)*zmd + si(i-1, j, k, 2)*ymd + si(i-1, j, &
&       k, 1)*xmd
      sid(i-1, j, k, 2) = sid(i-1, j, k, 2) + volmi*ymd
      sid(i-1, j, k, 1) = sid(i-1, j, k, 1) + volmi*xmd
      temp0 = vol(i, j, k) + vol(i+1, j, k)
      tempd0 = -(two*volpid/temp0**2)
      vold(i, j, k) = vold(i, j, k) + tempd0
      vold(i+1, j, k) = vold(i+1, j, k) + tempd0
      temp0 = vol(i, j, k) + vol(i-1, j, k)
      tempd0 = -(two*volmid/temp0**2)
      vold(i-1, j, k) = vold(i-1, j, k) + tempd0
      vold(i, j, k) = vold(i, j, k) + tempd0 - one*volid/vol(i, j, k)**2
    end do
!$bwd-of ii-loop 
    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 metrics in eta-direction, i.e. along the
! line j = constant.
      voli = one/vol(i, j, k)
      volmi = two/(vol(i, j, k)+vol(i, j-1, k))
      volpi = two/(vol(i, j, k)+vol(i, j+1, k))
      xm = sj(i, j-1, k, 1)*volmi
      ym = sj(i, j-1, k, 2)*volmi
      zm = sj(i, j-1, k, 3)*volmi
      xp = sj(i, j, k, 1)*volpi
      yp = sj(i, j, k, 2)*volpi
      zp = sj(i, j, k, 3)*volpi
      xa = half*(sj(i, j, k, 1)+sj(i, j-1, k, 1))*voli
      ya = half*(sj(i, j, k, 2)+sj(i, j-1, k, 2))*voli
      za = half*(sj(i, j, k, 3)+sj(i, j-1, k, 3))*voli
      ttm = xm*xa + ym*ya + zm*za
      ttp = xp*xa + yp*ya + zp*za
! computation of the viscous terms in eta-direction; note
! that cross-derivatives are neglected, i.e. the mesh is
! assumed to be orthogonal.
! furthermore, the grad(nu)**2 has been rewritten as
! div(nu grad(nu)) - nu div(grad nu) to enhance stability.
! the second derivative in eta-direction is constructed as
! the central difference of the first order derivatives, i.e.
! d^2/deta^2 = d/deta (d/deta j+1/2 - d/deta j-1/2).
! in this way the metric can be taken into account.
! compute the diffusion coefficients multiplying the nodes
! j+1, j and j-1 in the second derivative. make sure that
! these coefficients are nonnegative.
      cnud = -(rsacb2*w(i, j, k, itu1)*cb3inv)
      cam = ttm*cnud
      cap = ttp*cnud
      nutm = half*(w(i, j-1, k, itu1)+w(i, j, k, itu1))
      nutp = half*(w(i, j+1, k, itu1)+w(i, j, k, itu1))
      nu = rlv(i, j, k)/w(i, j, k, irho)
      num = half*(rlv(i, j-1, k)/w(i, j-1, k, irho)+nu)
      nup = half*(rlv(i, j+1, k)/w(i, j+1, k, irho)+nu)
      cdm = (num+(one+rsacb2)*nutm)*ttm*cb3inv
      cdp = (nup+(one+rsacb2)*nutp)*ttp*cb3inv
      if (cdm + cam .lt. zero) then
        c1m = zero
        call pushcontrol1b(0)
      else
        c1m = cdm + cam
        call pushcontrol1b(1)
      end if
      if (cdp + cap .lt. zero) then
        c1p = zero
        call pushcontrol1b(0)
      else
        c1p = cdp + cap
        call pushcontrol1b(1)
      end if
      c10 = c1m + c1p
! update the residual for this cell and store the possible
! coefficients for the matrix in b1, c1 and d1.
      c10d = -(w(i, j, k, itu1)*scratchd(i, j, k, idvt))
      c1md = w(i, j-1, k, itu1)*scratchd(i, j, k, idvt) + c10d
      wd(i, j-1, k, itu1) = wd(i, j-1, k, itu1) + c1m*scratchd(i, j, k, &
&       idvt)
      c1pd = w(i, j+1, k, itu1)*scratchd(i, j, k, idvt) + c10d
      wd(i, j+1, k, itu1) = wd(i, j+1, k, itu1) + c1p*scratchd(i, j, k, &
&       idvt)
      wd(i, j, k, itu1) = wd(i, j, k, itu1) - c10*scratchd(i, j, k, idvt&
&       )
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        cdpd = 0.0_8
        capd = 0.0_8
      else
        cdpd = c1pd
        capd = c1pd
      end if
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        cdmd = 0.0_8
        camd = 0.0_8
      else
        cdmd = c1md
        camd = c1md
      end if
      cnudd = ttp*capd + ttm*camd
      tempd0 = cb3inv*cdpd
      nupd = ttp*tempd0
      nutpd = (one+rsacb2)*ttp*tempd0
      ttpd = (nup+(one+rsacb2)*nutp)*tempd0 + cnud*capd
      tempd0 = cb3inv*cdmd
      numd = ttm*tempd0
      nutmd = (one+rsacb2)*ttm*tempd0
      ttmd = (num+(one+rsacb2)*nutm)*tempd0 + cnud*camd
      temp0 = w(i, j+1, k, irho)
      tempd = half*nupd/temp0
      nud = half*nupd + half*numd
      rlvd(i, j+1, k) = rlvd(i, j+1, k) + tempd
      wd(i, j+1, k, irho) = wd(i, j+1, k, irho) - rlv(i, j+1, k)*tempd/&
&       temp0
      temp0 = w(i, j-1, k, irho)
      tempd = half*numd/temp0
      rlvd(i, j-1, k) = rlvd(i, j-1, k) + tempd
      wd(i, j-1, k, irho) = wd(i, j-1, k, irho) - rlv(i, j-1, k)*tempd/&
&       temp0
      temp0 = w(i, j, k, irho)
      rlvd(i, j, k) = rlvd(i, j, k) + nud/temp0
      wd(i, j, k, irho) = wd(i, j, k, irho) - rlv(i, j, k)*nud/temp0**2
      wd(i, j+1, k, itu1) = wd(i, j+1, k, itu1) + half*nutpd
      wd(i, j, k, itu1) = wd(i, j, k, itu1) + half*nutpd + half*nutmd - &
&       rsacb2*cb3inv*cnudd
      wd(i, j-1, k, itu1) = wd(i, j-1, k, itu1) + half*nutmd
      xpd = xa*ttpd
      xad = xp*ttpd + xm*ttmd
      ypd = ya*ttpd
      yad = yp*ttpd + ym*ttmd
      zpd = za*ttpd
      zad = zp*ttpd + zm*ttmd
      xmd = xa*ttmd
      ymd = ya*ttmd
      zmd = za*ttmd
      tempd0 = half*zad
      sjd(i, j, k, 3) = sjd(i, j, k, 3) + voli*tempd0
      sjd(i, j-1, k, 3) = sjd(i, j-1, k, 3) + voli*tempd0
      volid = (sj(i, j, k, 3)+sj(i, j-1, k, 3))*tempd0
      tempd0 = half*yad
      sjd(i, j, k, 2) = sjd(i, j, k, 2) + voli*tempd0
      sjd(i, j-1, k, 2) = sjd(i, j-1, k, 2) + voli*tempd0
      volid = volid + (sj(i, j, k, 2)+sj(i, j-1, k, 2))*tempd0
      tempd0 = half*xad
      sjd(i, j, k, 1) = sjd(i, j, k, 1) + voli*tempd0
      sjd(i, j-1, k, 1) = sjd(i, j-1, k, 1) + voli*tempd0
      volid = volid + (sj(i, j, k, 1)+sj(i, j-1, k, 1))*tempd0
      sjd(i, j, k, 3) = sjd(i, j, k, 3) + volpi*zpd
      volpid = sj(i, j, k, 3)*zpd + sj(i, j, k, 2)*ypd + sj(i, j, k, 1)*&
&       xpd
      sjd(i, j, k, 2) = sjd(i, j, k, 2) + volpi*ypd
      sjd(i, j, k, 1) = sjd(i, j, k, 1) + volpi*xpd
      sjd(i, j-1, k, 3) = sjd(i, j-1, k, 3) + volmi*zmd
      volmid = sj(i, j-1, k, 3)*zmd + sj(i, j-1, k, 2)*ymd + sj(i, j-1, &
&       k, 1)*xmd
      sjd(i, j-1, k, 2) = sjd(i, j-1, k, 2) + volmi*ymd
      sjd(i, j-1, k, 1) = sjd(i, j-1, k, 1) + volmi*xmd
      temp0 = vol(i, j, k) + vol(i, j+1, k)
      tempd0 = -(two*volpid/temp0**2)
      vold(i, j, k) = vold(i, j, k) + tempd0
      vold(i, j+1, k) = vold(i, j+1, k) + tempd0
      temp0 = vol(i, j, k) + vol(i, j-1, k)
      tempd0 = -(two*volmid/temp0**2)
      vold(i, j-1, k) = vold(i, j-1, k) + tempd0
      vold(i, j, k) = vold(i, j, k) + tempd0 - one*volid/vol(i, j, k)**2
    end do
!$bwd-of ii-loop 
    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 metrics in zeta-direction, i.e. along the
! line k = constant.
      voli = one/vol(i, j, k)
      volmi = two/(vol(i, j, k)+vol(i, j, k-1))
      volpi = two/(vol(i, j, k)+vol(i, j, k+1))
      xm = sk(i, j, k-1, 1)*volmi
      ym = sk(i, j, k-1, 2)*volmi
      zm = sk(i, j, k-1, 3)*volmi
      xp = sk(i, j, k, 1)*volpi
      yp = sk(i, j, k, 2)*volpi
      zp = sk(i, j, k, 3)*volpi
      xa = half*(sk(i, j, k, 1)+sk(i, j, k-1, 1))*voli
      ya = half*(sk(i, j, k, 2)+sk(i, j, k-1, 2))*voli
      za = half*(sk(i, j, k, 3)+sk(i, j, k-1, 3))*voli
      ttm = xm*xa + ym*ya + zm*za
      ttp = xp*xa + yp*ya + zp*za
! ttm and ttp ~ 1/deltax^2
! computation of the viscous terms in zeta-direction; note
! that cross-derivatives are neglected, i.e. the mesh is
! assumed to be orthogonal.
! furthermore, the grad(nu)**2 has been rewritten as
! div(nu grad(nu)) - nu div(grad nu) to enhance stability.
! the second derivative in zeta-direction is constructed as
! the central difference of the first order derivatives, i.e.
! d^2/dzeta^2 = d/dzeta (d/dzeta k+1/2 - d/dzeta k-1/2).
! in this way the metric can be taken into account.
! compute the diffusion coefficients multiplying the nodes
! k+1, k and k-1 in the second derivative. make sure that
! these coefficients are nonnegative.
      cnud = -(rsacb2*w(i, j, k, itu1)*cb3inv)
      cam = ttm*cnud
      cap = ttp*cnud
! compute nutilde at the faces
      nutm = half*(w(i, j, k-1, itu1)+w(i, j, k, itu1))
      nutp = half*(w(i, j, k+1, itu1)+w(i, j, k, itu1))
! compute nu at the faces
      nu = rlv(i, j, k)/w(i, j, k, irho)
      num = half*(rlv(i, j, k-1)/w(i, j, k-1, irho)+nu)
      nup = half*(rlv(i, j, k+1)/w(i, j, k+1, irho)+nu)
      cdm = (num+(one+rsacb2)*nutm)*ttm*cb3inv
      cdp = (nup+(one+rsacb2)*nutp)*ttp*cb3inv
      if (cdm + cam .lt. zero) then
        c1m = zero
        call pushcontrol1b(0)
      else
        c1m = cdm + cam
        call pushcontrol1b(1)
      end if
      if (cdp + cap .lt. zero) then
        c1p = zero
        call pushcontrol1b(0)
      else
        c1p = cdp + cap
        call pushcontrol1b(1)
      end if
      c10 = c1m + c1p
! update the residual for this cell and store the possible
! coefficients for the matrix in b1, c1 and d1.
      c10d = -(w(i, j, k, itu1)*scratchd(i, j, k, idvt))
      c1md = w(i, j, k-1, itu1)*scratchd(i, j, k, idvt) + c10d
      wd(i, j, k-1, itu1) = wd(i, j, k-1, itu1) + c1m*scratchd(i, j, k, &
&       idvt)
      c1pd = w(i, j, k+1, itu1)*scratchd(i, j, k, idvt) + c10d
      wd(i, j, k+1, itu1) = wd(i, j, k+1, itu1) + c1p*scratchd(i, j, k, &
&       idvt)
      wd(i, j, k, itu1) = wd(i, j, k, itu1) - c10*scratchd(i, j, k, idvt&
&       )
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        cdpd = 0.0_8
        capd = 0.0_8
      else
        cdpd = c1pd
        capd = c1pd
      end if
      call popcontrol1b(branch)
      if (branch .eq. 0) then
        cdmd = 0.0_8
        camd = 0.0_8
      else
        cdmd = c1md
        camd = c1md
      end if
      cnudd = ttp*capd + ttm*camd
      tempd0 = cb3inv*cdpd
      nupd = ttp*tempd0
      nutpd = (one+rsacb2)*ttp*tempd0
      ttpd = (nup+(one+rsacb2)*nutp)*tempd0 + cnud*capd
      tempd0 = cb3inv*cdmd
      numd = ttm*tempd0
      nutmd = (one+rsacb2)*ttm*tempd0
      ttmd = (num+(one+rsacb2)*nutm)*tempd0 + cnud*camd
      temp0 = w(i, j, k+1, irho)
      tempd = half*nupd/temp0
      nud = half*nupd + half*numd
      rlvd(i, j, k+1) = rlvd(i, j, k+1) + tempd
      wd(i, j, k+1, irho) = wd(i, j, k+1, irho) - rlv(i, j, k+1)*tempd/&
&       temp0
      temp = w(i, j, k-1, irho)
      tempd0 = half*numd/temp
      rlvd(i, j, k-1) = rlvd(i, j, k-1) + tempd0
      wd(i, j, k-1, irho) = wd(i, j, k-1, irho) - rlv(i, j, k-1)*tempd0/&
&       temp
      temp = w(i, j, k, irho)
      rlvd(i, j, k) = rlvd(i, j, k) + nud/temp
      wd(i, j, k, irho) = wd(i, j, k, irho) - rlv(i, j, k)*nud/temp**2
      wd(i, j, k+1, itu1) = wd(i, j, k+1, itu1) + half*nutpd
      wd(i, j, k, itu1) = wd(i, j, k, itu1) + half*nutpd + half*nutmd - &
&       rsacb2*cb3inv*cnudd
      wd(i, j, k-1, itu1) = wd(i, j, k-1, itu1) + half*nutmd
      xpd = xa*ttpd
      xad = xp*ttpd + xm*ttmd
      ypd = ya*ttpd
      yad = yp*ttpd + ym*ttmd
      zpd = za*ttpd
      zad = zp*ttpd + zm*ttmd
      xmd = xa*ttmd
      ymd = ya*ttmd
      zmd = za*ttmd
      tempd = half*zad
      skd(i, j, k, 3) = skd(i, j, k, 3) + voli*tempd
      skd(i, j, k-1, 3) = skd(i, j, k-1, 3) + voli*tempd
      volid = (sk(i, j, k, 3)+sk(i, j, k-1, 3))*tempd
      tempd = half*yad
      skd(i, j, k, 2) = skd(i, j, k, 2) + voli*tempd
      skd(i, j, k-1, 2) = skd(i, j, k-1, 2) + voli*tempd
      volid = volid + (sk(i, j, k, 2)+sk(i, j, k-1, 2))*tempd
      tempd = half*xad
      skd(i, j, k, 1) = skd(i, j, k, 1) + voli*tempd
      skd(i, j, k-1, 1) = skd(i, j, k-1, 1) + voli*tempd
      volid = volid + (sk(i, j, k, 1)+sk(i, j, k-1, 1))*tempd
      skd(i, j, k, 3) = skd(i, j, k, 3) + volpi*zpd
      volpid = sk(i, j, k, 3)*zpd + sk(i, j, k, 2)*ypd + sk(i, j, k, 1)*&
&       xpd
      skd(i, j, k, 2) = skd(i, j, k, 2) + volpi*ypd
      skd(i, j, k, 1) = skd(i, j, k, 1) + volpi*xpd
      skd(i, j, k-1, 3) = skd(i, j, k-1, 3) + volmi*zmd
      volmid = sk(i, j, k-1, 3)*zmd + sk(i, j, k-1, 2)*ymd + sk(i, j, k-&
&       1, 1)*xmd
      skd(i, j, k-1, 2) = skd(i, j, k-1, 2) + volmi*ymd
      skd(i, j, k-1, 1) = skd(i, j, k-1, 1) + volmi*xmd
      temp = vol(i, j, k) + vol(i, j, k+1)
      tempd = -(two*volpid/temp**2)
      vold(i, j, k) = vold(i, j, k) + tempd
      vold(i, j, k+1) = vold(i, j, k+1) + tempd
      temp = vol(i, j, k) + vol(i, j, k-1)
      tempd = -(two*volmid/temp**2)
      vold(i, j, k-1) = vold(i, j, k-1) + tempd
      vold(i, j, k) = vold(i, j, k) + tempd - one*volid/vol(i, j, k)**2
    end do
  end subroutine saviscous_b
 
  subroutine saviscous()
!
!  viscous term.
!  determine the viscous contribution to the residual
!  for all internal cells of the block.
    use blockpointers
    use paramturb
    implicit none
! local variables.
    integer(kind=inttype) :: i, j, k, nn, ii
    real(kind=realtype) :: nu
    real(kind=realtype) :: fv1, fv2, ft2
    real(kind=realtype) :: voli, volmi, volpi, xm, ym, zm, xp, yp, zp
    real(kind=realtype) :: xa, ya, za, ttm, ttp, cnud, cam, cap
    real(kind=realtype) :: nutm, nutp, num, nup, cdm, cdp
    real(kind=realtype) :: c1m, c1p, c10, b1, c1, d1, qs
    intrinsic mod
    intrinsic max
! set model constants
    cv13 = rsacv1**3
    kar2inv = one/rsak**2
    cw36 = rsacw3**6
    cb3inv = one/rsacb3
!$ad ii-loop
!
!       viscous terms in k-direction.
!
    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 metrics in zeta-direction, i.e. along the
! line k = constant.
      voli = one/vol(i, j, k)
      volmi = two/(vol(i, j, k)+vol(i, j, k-1))
      volpi = two/(vol(i, j, k)+vol(i, j, k+1))
      xm = sk(i, j, k-1, 1)*volmi
      ym = sk(i, j, k-1, 2)*volmi
      zm = sk(i, j, k-1, 3)*volmi
      xp = sk(i, j, k, 1)*volpi
      yp = sk(i, j, k, 2)*volpi
      zp = sk(i, j, k, 3)*volpi
      xa = half*(sk(i, j, k, 1)+sk(i, j, k-1, 1))*voli
      ya = half*(sk(i, j, k, 2)+sk(i, j, k-1, 2))*voli
      za = half*(sk(i, j, k, 3)+sk(i, j, k-1, 3))*voli
      ttm = xm*xa + ym*ya + zm*za
      ttp = xp*xa + yp*ya + zp*za
! ttm and ttp ~ 1/deltax^2
! computation of the viscous terms in zeta-direction; note
! that cross-derivatives are neglected, i.e. the mesh is
! assumed to be orthogonal.
! furthermore, the grad(nu)**2 has been rewritten as
! div(nu grad(nu)) - nu div(grad nu) to enhance stability.
! the second derivative in zeta-direction is constructed as
! the central difference of the first order derivatives, i.e.
! d^2/dzeta^2 = d/dzeta (d/dzeta k+1/2 - d/dzeta k-1/2).
! in this way the metric can be taken into account.
! compute the diffusion coefficients multiplying the nodes
! k+1, k and k-1 in the second derivative. make sure that
! these coefficients are nonnegative.
      cnud = -(rsacb2*w(i, j, k, itu1)*cb3inv)
      cam = ttm*cnud
      cap = ttp*cnud
! compute nutilde at the faces
      nutm = half*(w(i, j, k-1, itu1)+w(i, j, k, itu1))
      nutp = half*(w(i, j, k+1, itu1)+w(i, j, k, itu1))
! compute nu at the faces
      nu = rlv(i, j, k)/w(i, j, k, irho)
      num = half*(rlv(i, j, k-1)/w(i, j, k-1, irho)+nu)
      nup = half*(rlv(i, j, k+1)/w(i, j, k+1, irho)+nu)
      cdm = (num+(one+rsacb2)*nutm)*ttm*cb3inv
      cdp = (nup+(one+rsacb2)*nutp)*ttp*cb3inv
      if (cdm + cam .lt. zero) then
        c1m = zero
      else
        c1m = cdm + cam
      end if
      if (cdp + cap .lt. zero) then
        c1p = zero
      else
        c1p = cdp + cap
      end if
      c10 = c1m + c1p
! update the residual for this cell and store the possible
! coefficients for the matrix in b1, c1 and d1.
      scratch(i, j, k, idvt) = scratch(i, j, k, idvt) + c1m*w(i, j, k-1&
&       , itu1) - c10*w(i, j, k, itu1) + c1p*w(i, j, k+1, itu1)
    end do
!$ad ii-loop
!
!       viscous terms in j-direction.
!
    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 metrics in eta-direction, i.e. along the
! line j = constant.
      voli = one/vol(i, j, k)
      volmi = two/(vol(i, j, k)+vol(i, j-1, k))
      volpi = two/(vol(i, j, k)+vol(i, j+1, k))
      xm = sj(i, j-1, k, 1)*volmi
      ym = sj(i, j-1, k, 2)*volmi
      zm = sj(i, j-1, k, 3)*volmi
      xp = sj(i, j, k, 1)*volpi
      yp = sj(i, j, k, 2)*volpi
      zp = sj(i, j, k, 3)*volpi
      xa = half*(sj(i, j, k, 1)+sj(i, j-1, k, 1))*voli
      ya = half*(sj(i, j, k, 2)+sj(i, j-1, k, 2))*voli
      za = half*(sj(i, j, k, 3)+sj(i, j-1, k, 3))*voli
      ttm = xm*xa + ym*ya + zm*za
      ttp = xp*xa + yp*ya + zp*za
! computation of the viscous terms in eta-direction; note
! that cross-derivatives are neglected, i.e. the mesh is
! assumed to be orthogonal.
! furthermore, the grad(nu)**2 has been rewritten as
! div(nu grad(nu)) - nu div(grad nu) to enhance stability.
! the second derivative in eta-direction is constructed as
! the central difference of the first order derivatives, i.e.
! d^2/deta^2 = d/deta (d/deta j+1/2 - d/deta j-1/2).
! in this way the metric can be taken into account.
! compute the diffusion coefficients multiplying the nodes
! j+1, j and j-1 in the second derivative. make sure that
! these coefficients are nonnegative.
      cnud = -(rsacb2*w(i, j, k, itu1)*cb3inv)
      cam = ttm*cnud
      cap = ttp*cnud
      nutm = half*(w(i, j-1, k, itu1)+w(i, j, k, itu1))
      nutp = half*(w(i, j+1, k, itu1)+w(i, j, k, itu1))
      nu = rlv(i, j, k)/w(i, j, k, irho)
      num = half*(rlv(i, j-1, k)/w(i, j-1, k, irho)+nu)
      nup = half*(rlv(i, j+1, k)/w(i, j+1, k, irho)+nu)
      cdm = (num+(one+rsacb2)*nutm)*ttm*cb3inv
      cdp = (nup+(one+rsacb2)*nutp)*ttp*cb3inv
      if (cdm + cam .lt. zero) then
        c1m = zero
      else
        c1m = cdm + cam
      end if
      if (cdp + cap .lt. zero) then
        c1p = zero
      else
        c1p = cdp + cap
      end if
      c10 = c1m + c1p
! update the residual for this cell and store the possible
! coefficients for the matrix in b1, c1 and d1.
      scratch(i, j, k, idvt) = scratch(i, j, k, idvt) + c1m*w(i, j-1, k&
&       , itu1) - c10*w(i, j, k, itu1) + c1p*w(i, j+1, k, itu1)
    end do
!$ad ii-loop
!
!       viscous terms in i-direction.
!
    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 metrics in xi-direction, i.e. along the
! line i = constant.
      voli = one/vol(i, j, k)
      volmi = two/(vol(i, j, k)+vol(i-1, j, k))
      volpi = two/(vol(i, j, k)+vol(i+1, j, k))
      xm = si(i-1, j, k, 1)*volmi
      ym = si(i-1, j, k, 2)*volmi
      zm = si(i-1, j, k, 3)*volmi
      xp = si(i, j, k, 1)*volpi
      yp = si(i, j, k, 2)*volpi
      zp = si(i, j, k, 3)*volpi
      xa = half*(si(i, j, k, 1)+si(i-1, j, k, 1))*voli
      ya = half*(si(i, j, k, 2)+si(i-1, j, k, 2))*voli
      za = half*(si(i, j, k, 3)+si(i-1, j, k, 3))*voli
      ttm = xm*xa + ym*ya + zm*za
      ttp = xp*xa + yp*ya + zp*za
! computation of the viscous terms in xi-direction; note
! that cross-derivatives are neglected, i.e. the mesh is
! assumed to be orthogonal.
! furthermore, the grad(nu)**2 has been rewritten as
! div(nu grad(nu)) - nu div(grad nu) to enhance stability.
! the second derivative in xi-direction is constructed as
! the central difference of the first order derivatives, i.e.
! d^2/dxi^2 = d/dxi (d/dxi i+1/2 - d/dxi i-1/2).
! in this way the metric can be taken into account.
! compute the diffusion coefficients multiplying the nodes
! i+1, i and i-1 in the second derivative. make sure that
! these coefficients are nonnegative.
      cnud = -(rsacb2*w(i, j, k, itu1)*cb3inv)
      cam = ttm*cnud
      cap = ttp*cnud
      nutm = half*(w(i-1, j, k, itu1)+w(i, j, k, itu1))
      nutp = half*(w(i+1, j, k, itu1)+w(i, j, k, itu1))
      nu = rlv(i, j, k)/w(i, j, k, irho)
      num = half*(rlv(i-1, j, k)/w(i-1, j, k, irho)+nu)
      nup = half*(rlv(i+1, j, k)/w(i+1, j, k, irho)+nu)
      cdm = (num+(one+rsacb2)*nutm)*ttm*cb3inv
      cdp = (nup+(one+rsacb2)*nutp)*ttp*cb3inv
      if (cdm + cam .lt. zero) then
        c1m = zero
      else
        c1m = cdm + cam
      end if
      if (cdp + cap .lt. zero) then
        c1p = zero
      else
        c1p = cdp + cap
      end if
      c10 = c1m + c1p
! update the residual for this cell and store the possible
! coefficients for the matrix in b1, c1 and d1.
      scratch(i, j, k, idvt) = scratch(i, j, k, idvt) + c1m*w(i-1, j, k&
&       , itu1) - c10*w(i, j, k, itu1) + c1p*w(i+1, j, k, itu1)
    end do
  end subroutine saviscous
 
!  differentiation of saresscale in reverse (adjoint) mode (with options noisize i4 dr8 r8):
!   gradient     of useful results: *dw
!   with respect to varying inputs: *dw *scratch
!   rw status of diff variables: *dw:in-out *scratch:out
!   plus diff mem management of: dw:in scratch:in
  subroutine saresscale_b()
!
!  multiply the residual by the volume and store this in dw; this
! * is done for monitoring reasons only. the multiplication with the
! * volume is present to be consistent with the flow residuals; also
!  the negative value is taken, again to be consistent with the
! * flow equations. also multiply by iblank so that no updates occur
!  in holes or the overset boundary.
    use blockpointers
    implicit none
! local variables
    integer(kind=inttype) :: i, j, k, ii
    real(kind=realtype) :: rblank
    intrinsic mod
    intrinsic real
    intrinsic max
    real(realtype) :: x1
    if (associated(scratchd)) scratchd = 0.0_8
!$bwd-of ii-loop 
    do ii=0,nx*ny*nz-1
      i = mod(ii, nx) + 2
      j = mod(ii/nx, ny) + 2
      k = ii/(nx*ny) + 2
      x1 = real(iblank(i, j, k), realtype)
      if (x1 .lt. zero) then
        rblank = zero
      else
        rblank = x1
      end if
      scratchd(i, j, k, idvt) = scratchd(i, j, k, idvt) - volref(i, j, k&
&       )*rblank*dwd(i, j, k, itu1)
      dwd(i, j, k, itu1) = 0.0_8
    end do
  end subroutine saresscale_b
 
  subroutine saresscale()
!
!  multiply the residual by the volume and store this in dw; this
! * is done for monitoring reasons only. the multiplication with the
! * volume is present to be consistent with the flow residuals; also
!  the negative value is taken, again to be consistent with the
! * flow equations. also multiply by iblank so that no updates occur
!  in holes or the overset boundary.
    use blockpointers
    implicit none
! local variables
    integer(kind=inttype) :: i, j, k, ii
    real(kind=realtype) :: rblank
    intrinsic mod
    intrinsic real
    intrinsic max
    real(realtype) :: x1
!$ad ii-loop
    do ii=0,nx*ny*nz-1
      i = mod(ii, nx) + 2
      j = mod(ii/nx, ny) + 2
      k = ii/(nx*ny) + 2
      x1 = real(iblank(i, j, k), realtype)
      if (x1 .lt. zero) then
        rblank = zero
      else
        rblank = x1
      end if
      dw(i, j, k, itu1) = -(volref(i, j, k)*scratch(i, j, k, idvt)*&
&       rblank)
    end do
  end subroutine saresscale
 
end module sa_b