sa_d.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_d
  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 forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *scratch
!   with respect to varying inputs: timeref *w *rlv *vol *d2wall
!                *si *sj *sk
!   rw status of diff variables: timeref:in *w:in *rlv:in *scratch:out
!                *vol:in *d2wall:in *si:in *sj:in *sk:in
!   plus diff mem management of: w:in rlv:in scratch:in vol:in
!                d2wall:in si:in sj:in sk:in
  subroutine sasource_d()
!
!  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 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) :: temp0
    real(kind=realtype) :: temp1
    real(kind=realtype) :: temp2
    real(kind=realtype) :: temp3
    real(kind=realtype) :: temp4
    real(kind=realtype) :: temp5
    real(kind=realtype) :: temp6
    real(kind=realtype) :: temp7
    real(kind=realtype) :: temp8
    real(kind=realtype) :: temp9
    real(kind=realtype) :: temp10
! 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.
    omegaxd = sections(sectionid)%rotrate(1)*timerefd
    omegax = timeref*sections(sectionid)%rotrate(1)
    omegayd = sections(sectionid)%rotrate(2)*timerefd
    omegay = timeref*sections(sectionid)%rotrate(2)
    omegazd = sections(sectionid)%rotrate(3)*timerefd
    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
      if (associated(scratchd)) scratchd = 0.0_8
      strainmag2d = 0.0_8
      ssd = 0.0_8
      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.
            temp = si(i, j, k, 1)
            temp0 = w(i+1, j, k, ivx)
            temp1 = si(i-1, j, k, 1)
            temp2 = w(i-1, j, k, ivx)
            temp3 = sj(i, j, k, 1)
            temp4 = w(i, j+1, k, ivx)
            temp5 = sk(i, j, k, 1)
            temp6 = w(i, j, k+1, ivx)
            temp7 = sj(i, j-1, k, 1)
            temp8 = w(i, j-1, k, ivx)
            temp9 = sk(i, j, k-1, 1)
            temp10 = w(i, j, k-1, ivx)
            uuxd = temp*wd(i+1, j, k, ivx) + temp0*sid(i, j, k, 1) - &
&             temp1*wd(i-1, j, k, ivx) - temp2*sid(i-1, j, k, 1) + temp3&
&             *wd(i, j+1, k, ivx) + temp4*sjd(i, j, k, 1) + temp5*wd(i, &
&             j, k+1, ivx) + temp6*skd(i, j, k, 1) - temp7*wd(i, j-1, k&
&             , ivx) - temp8*sjd(i, j-1, k, 1) - temp9*wd(i, j, k-1, ivx&
&             ) - temp10*skd(i, j, k-1, 1)
            uux = temp0*temp - temp2*temp1 + temp4*temp3 + temp6*temp5 -&
&             temp8*temp7 - temp10*temp9
            temp10 = si(i, j, k, 2)
            temp9 = w(i+1, j, k, ivx)
            temp8 = si(i-1, j, k, 2)
            temp7 = w(i-1, j, k, ivx)
            temp6 = sj(i, j, k, 2)
            temp5 = w(i, j+1, k, ivx)
            temp4 = sk(i, j, k, 2)
            temp3 = w(i, j, k+1, ivx)
            temp2 = sj(i, j-1, k, 2)
            temp1 = w(i, j-1, k, ivx)
            temp0 = sk(i, j, k-1, 2)
            temp = w(i, j, k-1, ivx)
            uuyd = temp10*wd(i+1, j, k, ivx) + temp9*sid(i, j, k, 2) - &
&             temp8*wd(i-1, j, k, ivx) - temp7*sid(i-1, j, k, 2) + temp6&
&             *wd(i, j+1, k, ivx) + temp5*sjd(i, j, k, 2) + temp4*wd(i, &
&             j, k+1, ivx) + temp3*skd(i, j, k, 2) - temp2*wd(i, j-1, k&
&             , ivx) - temp1*sjd(i, j-1, k, 2) - temp0*wd(i, j, k-1, ivx&
&             ) - temp*skd(i, j, k-1, 2)
            uuy = temp9*temp10 - temp7*temp8 + temp5*temp6 + temp3*temp4&
&             - temp1*temp2 - temp*temp0
            temp10 = si(i, j, k, 3)
            temp9 = w(i+1, j, k, ivx)
            temp8 = si(i-1, j, k, 3)
            temp7 = w(i-1, j, k, ivx)
            temp6 = sj(i, j, k, 3)
            temp5 = w(i, j+1, k, ivx)
            temp4 = sk(i, j, k, 3)
            temp3 = w(i, j, k+1, ivx)
            temp2 = sj(i, j-1, k, 3)
            temp1 = w(i, j-1, k, ivx)
            temp0 = sk(i, j, k-1, 3)
            temp = w(i, j, k-1, ivx)
            uuzd = temp10*wd(i+1, j, k, ivx) + temp9*sid(i, j, k, 3) - &
&             temp8*wd(i-1, j, k, ivx) - temp7*sid(i-1, j, k, 3) + temp6&
&             *wd(i, j+1, k, ivx) + temp5*sjd(i, j, k, 3) + temp4*wd(i, &
&             j, k+1, ivx) + temp3*skd(i, j, k, 3) - temp2*wd(i, j-1, k&
&             , ivx) - temp1*sjd(i, j-1, k, 3) - temp0*wd(i, j, k-1, ivx&
&             ) - temp*skd(i, j, k-1, 3)
            uuz = temp9*temp10 - temp7*temp8 + temp5*temp6 + temp3*temp4&
&             - temp1*temp2 - temp*temp0
! idem for the gradient of v.
            temp10 = si(i, j, k, 1)
            temp9 = w(i+1, j, k, ivy)
            temp8 = si(i-1, j, k, 1)
            temp7 = w(i-1, j, k, ivy)
            temp6 = sj(i, j, k, 1)
            temp5 = w(i, j+1, k, ivy)
            temp4 = sk(i, j, k, 1)
            temp3 = w(i, j, k+1, ivy)
            temp2 = sj(i, j-1, k, 1)
            temp1 = w(i, j-1, k, ivy)
            temp0 = sk(i, j, k-1, 1)
            temp = w(i, j, k-1, ivy)
            vvxd = temp10*wd(i+1, j, k, ivy) + temp9*sid(i, j, k, 1) - &
&             temp8*wd(i-1, j, k, ivy) - temp7*sid(i-1, j, k, 1) + temp6&
&             *wd(i, j+1, k, ivy) + temp5*sjd(i, j, k, 1) + temp4*wd(i, &
&             j, k+1, ivy) + temp3*skd(i, j, k, 1) - temp2*wd(i, j-1, k&
&             , ivy) - temp1*sjd(i, j-1, k, 1) - temp0*wd(i, j, k-1, ivy&
&             ) - temp*skd(i, j, k-1, 1)
            vvx = temp9*temp10 - temp7*temp8 + temp5*temp6 + temp3*temp4&
&             - temp1*temp2 - temp*temp0
            temp10 = si(i, j, k, 2)
            temp9 = w(i+1, j, k, ivy)
            temp8 = si(i-1, j, k, 2)
            temp7 = w(i-1, j, k, ivy)
            temp6 = sj(i, j, k, 2)
            temp5 = w(i, j+1, k, ivy)
            temp4 = sk(i, j, k, 2)
            temp3 = w(i, j, k+1, ivy)
            temp2 = sj(i, j-1, k, 2)
            temp1 = w(i, j-1, k, ivy)
            temp0 = sk(i, j, k-1, 2)
            temp = w(i, j, k-1, ivy)
            vvyd = temp10*wd(i+1, j, k, ivy) + temp9*sid(i, j, k, 2) - &
&             temp8*wd(i-1, j, k, ivy) - temp7*sid(i-1, j, k, 2) + temp6&
&             *wd(i, j+1, k, ivy) + temp5*sjd(i, j, k, 2) + temp4*wd(i, &
&             j, k+1, ivy) + temp3*skd(i, j, k, 2) - temp2*wd(i, j-1, k&
&             , ivy) - temp1*sjd(i, j-1, k, 2) - temp0*wd(i, j, k-1, ivy&
&             ) - temp*skd(i, j, k-1, 2)
            vvy = temp9*temp10 - temp7*temp8 + temp5*temp6 + temp3*temp4&
&             - temp1*temp2 - temp*temp0
            temp10 = si(i, j, k, 3)
            temp9 = w(i+1, j, k, ivy)
            temp8 = si(i-1, j, k, 3)
            temp7 = w(i-1, j, k, ivy)
            temp6 = sj(i, j, k, 3)
            temp5 = w(i, j+1, k, ivy)
            temp4 = sk(i, j, k, 3)
            temp3 = w(i, j, k+1, ivy)
            temp2 = sj(i, j-1, k, 3)
            temp1 = w(i, j-1, k, ivy)
            temp0 = sk(i, j, k-1, 3)
            temp = w(i, j, k-1, ivy)
            vvzd = temp10*wd(i+1, j, k, ivy) + temp9*sid(i, j, k, 3) - &
&             temp8*wd(i-1, j, k, ivy) - temp7*sid(i-1, j, k, 3) + temp6&
&             *wd(i, j+1, k, ivy) + temp5*sjd(i, j, k, 3) + temp4*wd(i, &
&             j, k+1, ivy) + temp3*skd(i, j, k, 3) - temp2*wd(i, j-1, k&
&             , ivy) - temp1*sjd(i, j-1, k, 3) - temp0*wd(i, j, k-1, ivy&
&             ) - temp*skd(i, j, k-1, 3)
            vvz = temp9*temp10 - temp7*temp8 + temp5*temp6 + temp3*temp4&
&             - temp1*temp2 - temp*temp0
! and for the gradient of w.
            temp10 = si(i, j, k, 1)
            temp9 = w(i+1, j, k, ivz)
            temp8 = si(i-1, j, k, 1)
            temp7 = w(i-1, j, k, ivz)
            temp6 = sj(i, j, k, 1)
            temp5 = w(i, j+1, k, ivz)
            temp4 = sk(i, j, k, 1)
            temp3 = w(i, j, k+1, ivz)
            temp2 = sj(i, j-1, k, 1)
            temp1 = w(i, j-1, k, ivz)
            temp0 = sk(i, j, k-1, 1)
            temp = w(i, j, k-1, ivz)
            wwxd = temp10*wd(i+1, j, k, ivz) + temp9*sid(i, j, k, 1) - &
&             temp8*wd(i-1, j, k, ivz) - temp7*sid(i-1, j, k, 1) + temp6&
&             *wd(i, j+1, k, ivz) + temp5*sjd(i, j, k, 1) + temp4*wd(i, &
&             j, k+1, ivz) + temp3*skd(i, j, k, 1) - temp2*wd(i, j-1, k&
&             , ivz) - temp1*sjd(i, j-1, k, 1) - temp0*wd(i, j, k-1, ivz&
&             ) - temp*skd(i, j, k-1, 1)
            wwx = temp9*temp10 - temp7*temp8 + temp5*temp6 + temp3*temp4&
&             - temp1*temp2 - temp*temp0
            temp10 = si(i, j, k, 2)
            temp9 = w(i+1, j, k, ivz)
            temp8 = si(i-1, j, k, 2)
            temp7 = w(i-1, j, k, ivz)
            temp6 = sj(i, j, k, 2)
            temp5 = w(i, j+1, k, ivz)
            temp4 = sk(i, j, k, 2)
            temp3 = w(i, j, k+1, ivz)
            temp2 = sj(i, j-1, k, 2)
            temp1 = w(i, j-1, k, ivz)
            temp0 = sk(i, j, k-1, 2)
            temp = w(i, j, k-1, ivz)
            wwyd = temp10*wd(i+1, j, k, ivz) + temp9*sid(i, j, k, 2) - &
&             temp8*wd(i-1, j, k, ivz) - temp7*sid(i-1, j, k, 2) + temp6&
&             *wd(i, j+1, k, ivz) + temp5*sjd(i, j, k, 2) + temp4*wd(i, &
&             j, k+1, ivz) + temp3*skd(i, j, k, 2) - temp2*wd(i, j-1, k&
&             , ivz) - temp1*sjd(i, j-1, k, 2) - temp0*wd(i, j, k-1, ivz&
&             ) - temp*skd(i, j, k-1, 2)
            wwy = temp9*temp10 - temp7*temp8 + temp5*temp6 + temp3*temp4&
&             - temp1*temp2 - temp*temp0
            temp10 = si(i, j, k, 3)
            temp9 = w(i+1, j, k, ivz)
            temp8 = si(i-1, j, k, 3)
            temp7 = w(i-1, j, k, ivz)
            temp6 = sj(i, j, k, 3)
            temp5 = w(i, j+1, k, ivz)
            temp4 = sk(i, j, k, 3)
            temp3 = w(i, j, k+1, ivz)
            temp2 = sj(i, j-1, k, 3)
            temp1 = w(i, j-1, k, ivz)
            temp0 = sk(i, j, k-1, 3)
            temp = w(i, j, k-1, ivz)
            wwzd = temp10*wd(i+1, j, k, ivz) + temp9*sid(i, j, k, 3) - &
&             temp8*wd(i-1, j, k, ivz) - temp7*sid(i-1, j, k, 3) + temp6&
&             *wd(i, j+1, k, ivz) + temp5*sjd(i, j, k, 3) + temp4*wd(i, &
&             j, k+1, ivz) + temp3*skd(i, j, k, 3) - temp2*wd(i, j-1, k&
&             , ivz) - temp1*sjd(i, j-1, k, 3) - temp0*wd(i, j, k-1, ivz&
&             ) - temp*skd(i, j, k-1, 3)
            wwz = temp9*temp10 - temp7*temp8 + temp5*temp6 + temp3*temp4&
&             - temp1*temp2 - temp*temp0
! 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).
            temp10 = fourth/vol(i, j, k)
            factd = -(temp10*vold(i, j, k)/vol(i, j, k))
            fact = temp10
            if (turbprod .eq. strain) then
              sxxd = two*(uux*factd+fact*uuxd)
              sxx = two*fact*uux
              syyd = two*(vvy*factd+fact*vvyd)
              syy = two*fact*vvy
              szzd = two*(wwz*factd+fact*wwzd)
              szz = two*fact*wwz
              sxyd = (uuy+vvx)*factd + fact*(uuyd+vvxd)
              sxy = fact*(uuy+vvx)
              sxzd = (uuz+wwx)*factd + fact*(uuzd+wwxd)
              sxz = fact*(uuz+wwx)
              syzd = (vvz+wwy)*factd + fact*(vvzd+wwyd)
              syz = fact*(vvz+wwy)
! compute 2/3 * divergence of velocity squared
              div2d = f23*2*(sxx+syy+szz)*(sxxd+syyd+szzd)
              div2 = f23*(sxx+syy+szz)**2
! compute strain production term
              strainmag2d = two*(2*sxy*sxyd+2*sxz*sxzd+2*syz*syzd) + 2*&
&               sxx*sxxd + 2*syy*syyd + 2*szz*szzd
              strainmag2 = two*(sxy**2+sxz**2+syz**2) + sxx**2 + syy**2 &
&               + szz**2
              strainprodd = two*strainmag2d - div2d
              strainprod = two*strainmag2 - div2
              temp10 = sqrt(strainprod)
              if (strainprod .eq. 0.0_8) then
                ssd = 0.0_8
              else
                ssd = strainprodd/(2.0*temp10)
              end if
              ss = temp10
            else if (turbprod .eq. vorticity) then
! compute the three components of the vorticity vector.
! substract the part coming from the rotating frame.
              vortxd = two*((wwy-vvz)*factd+fact*(wwyd-vvzd)) - two*&
&               omegaxd
              vortx = two*fact*(wwy-vvz) - two*omegax
              vortyd = two*((uuz-wwx)*factd+fact*(uuzd-wwxd)) - two*&
&               omegayd
              vorty = two*fact*(uuz-wwx) - two*omegay
              vortzd = two*((vvx-uuy)*factd+fact*(vvxd-uuyd)) - two*&
&               omegazd
              vortz = two*fact*(vvx-uuy) - two*omegaz
! compute the vorticity production term
              vortprodd = 2*vortx*vortxd + 2*vorty*vortyd + 2*vortz*&
&               vortzd
              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.
              temp10 = sqrt(vortprod)
              if (vortprod .eq. 0.0_8) then
                ssd = 0.0_8
              else
                ssd = vortprodd/(2.0*temp10)
              end if
              ss = temp10
            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.
            temp10 = w(i, j, k, irho)
            temp9 = rlv(i, j, k)/temp10
            nud = (rlvd(i, j, k)-temp9*wd(i, j, k, irho))/temp10
            nu = temp9
            temp10 = one/(d2wall(i, j, k)*d2wall(i, j, k))
            dist2invd = -(temp10*2*d2walld(i, j, k)/d2wall(i, j, k))
            dist2inv = temp10
            temp10 = w(i, j, k, itu1)/nu
            chid = (wd(i, j, k, itu1)-temp10*nud)/nu
            chi = temp10
            chi2d = 2*chi*chid
            chi2 = chi*chi
            chi3d = chi2*chid + chi*chi2d
            chi3 = chi*chi2
            temp10 = chi3/(cv13+chi3)
            fv1d = (1.0-temp10)*chi3d/(cv13+chi3)
            fv1 = temp10
            temp10 = chi/(one+chi*fv1)
            fv2d = -((chid-temp10*(fv1*chid+chi*fv1d))/(one+chi*fv1))
            fv2 = one - temp10
! 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
              ft2d = -(rsact3*exp(-(rsact4*chi2))*rsact4*chi2d)
              ft2 = rsact3*exp(-(rsact4*chi2))
            else
              ft2 = zero
              ft2d = 0.0_8
            end if
! correct the production term to account for the influence
! of the wall.
            temp10 = w(i, j, k, itu1)
            sstd = ssd + kar2inv*(fv2*dist2inv*wd(i, j, k, itu1)+temp10*&
&             (dist2inv*fv2d+fv2*dist2invd))
            sst = ss + kar2inv*(temp10*(fv2*dist2inv))
! add rotation term (userotationsa defined in inputparams.f90)
            if (userotationsa) then
              temp10 = sqrt(two*strainmag2)
              if (two*strainmag2 .eq. 0.0_8) then
                y1d = 0.0_8
              else
                y1d = two*strainmag2d/(2.0*temp10)
              end if
              y1 = temp10
              if (zero .gt. y1) then
                min1d = y1d
                min1 = y1
              else
                min1 = zero
                min1d = 0.0_8
              end if
              sstd = sstd + rsacrot*min1d
              sst = sst + rsacrot*min1
            end if
            if (sst .lt. xminn) then
              sst = xminn
              sstd = 0.0_8
            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.
            temp10 = dist2inv/sst
            temp9 = w(i, j, k, itu1)
            rrd = kar2inv*(temp10*wd(i, j, k, itu1)+temp9*(dist2invd-&
&             temp10*sstd)/sst)
            rr = kar2inv*(temp9*temp10)
            if (rr .gt. 10.0_realtype) then
              rr = 10.0_realtype
              rrd = 0.0_8
            else
              rr = rr
            end if
            ggd = (rsacw2*(6*rr**5-1.0)+1.0)*rrd
            gg = rr + rsacw2*(rr**6-rr)
            gg6d = 6*gg**5*ggd
            gg6 = gg**6
            temp10 = (one+cw36)/(cw36+gg6)
            if (temp10 .le. 0.0_8 .and. (sixth .eq. 0.0_8 .or. sixth &
&               .ne. int(sixth))) then
              termfwd = 0.0_8
            else
              termfwd = -(sixth*temp10**(sixth-1)*temp10*gg6d/(cw36+gg6)&
&               )
            end if
            termfw = temp10**sixth
            fwsad = termfw*ggd + gg*termfwd
            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
              term1d = 0.0_8
            else
              term1d = rsacb1*((one-ft2)*ssd-ss*ft2d)
              term1 = rsacb1*(one-ft2)*ss
            end if
            temp10 = kar2inv*rsacb1*((one-ft2)*fv2+ft2) - rsacw1*fwsa
            term2d = temp10*dist2invd + dist2inv*(kar2inv*rsacb1*((one-&
&             ft2)*fv2d-(fv2-1.0)*ft2d)-rsacw1*fwsad)
            term2 = dist2inv*temp10
            temp10 = w(i, j, k, itu1)
            temp9 = w(i, j, k, itu1)
            scratchd(i, j, k, idvt) = temp10*(term1d+temp9*term2d+term2*&
&             wd(i, j, k, itu1)) + (term1+term2*temp9)*wd(i, j, k, itu1)
            scratch(i, j, k, idvt) = (term1+term2*temp9)*temp10
          end do
        end do
      end do
    end if
  end subroutine sasource_d
 
  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 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
      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)
            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 do
      end do
    end if
  end subroutine sasource
 
!  differentiation of saviscous in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *scratch
!   with respect to varying inputs: *w *rlv *scratch *vol *si *sj
!                *sk
!   rw status of diff variables: *w:in *rlv:in *scratch:in-out
!                *vol:in *si:in *sj:in *sk:in
!   plus diff mem management of: w:in rlv:in scratch:in vol:in
!                si:in sj:in sk:in
  subroutine saviscous_d()
!
!  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 max
    real(kind=realtype) :: temp
    real(kind=realtype) :: temp0
    real(kind=realtype) :: temp1
! set model constants
    cv13 = rsacv1**3
    kar2inv = one/rsak**2
    cw36 = rsacw3**6
    cb3inv = one/rsacb3
!
!       viscous terms in k-direction.
!
    do k=2,kl
      do j=2,jl
        do i=2,il
! compute the metrics in zeta-direction, i.e. along the
! line k = constant.
          temp = one/vol(i, j, k)
          volid = -(temp*vold(i, j, k)/vol(i, j, k))
          voli = temp
          temp = two/(vol(i, j, k)+vol(i, j, k-1))
          volmid = -(temp*(vold(i, j, k)+vold(i, j, k-1))/(vol(i, j, k)+&
&           vol(i, j, k-1)))
          volmi = temp
          temp = two/(vol(i, j, k)+vol(i, j, k+1))
          volpid = -(temp*(vold(i, j, k)+vold(i, j, k+1))/(vol(i, j, k)+&
&           vol(i, j, k+1)))
          volpi = temp
          temp = sk(i, j, k-1, 1)
          xmd = volmi*skd(i, j, k-1, 1) + temp*volmid
          xm = temp*volmi
          temp = sk(i, j, k-1, 2)
          ymd = volmi*skd(i, j, k-1, 2) + temp*volmid
          ym = temp*volmi
          temp = sk(i, j, k-1, 3)
          zmd = volmi*skd(i, j, k-1, 3) + temp*volmid
          zm = temp*volmi
          temp = sk(i, j, k, 1)
          xpd = volpi*skd(i, j, k, 1) + temp*volpid
          xp = temp*volpi
          temp = sk(i, j, k, 2)
          ypd = volpi*skd(i, j, k, 2) + temp*volpid
          yp = temp*volpi
          temp = sk(i, j, k, 3)
          zpd = volpi*skd(i, j, k, 3) + temp*volpid
          zp = temp*volpi
          temp = sk(i, j, k, 1) + sk(i, j, k-1, 1)
          xad = half*(voli*(skd(i, j, k, 1)+skd(i, j, k-1, 1))+temp*&
&           volid)
          xa = half*(temp*voli)
          temp = sk(i, j, k, 2) + sk(i, j, k-1, 2)
          yad = half*(voli*(skd(i, j, k, 2)+skd(i, j, k-1, 2))+temp*&
&           volid)
          ya = half*(temp*voli)
          temp = sk(i, j, k, 3) + sk(i, j, k-1, 3)
          zad = half*(voli*(skd(i, j, k, 3)+skd(i, j, k-1, 3))+temp*&
&           volid)
          za = half*(temp*voli)
          ttmd = xa*xmd + xm*xad + ya*ymd + ym*yad + za*zmd + zm*zad
          ttm = xm*xa + ym*ya + zm*za
          ttpd = xa*xpd + xp*xad + ya*ypd + yp*yad + za*zpd + zp*zad
          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.
          cnudd = -(rsacb2*cb3inv*wd(i, j, k, itu1))
          cnud = -(rsacb2*w(i, j, k, itu1)*cb3inv)
          camd = cnud*ttmd + ttm*cnudd
          cam = ttm*cnud
          capd = cnud*ttpd + ttp*cnudd
          cap = ttp*cnud
! compute nutilde at the faces
          nutmd = half*(wd(i, j, k-1, itu1)+wd(i, j, k, itu1))
          nutm = half*(w(i, j, k-1, itu1)+w(i, j, k, itu1))
          nutpd = half*(wd(i, j, k+1, itu1)+wd(i, j, k, itu1))
          nutp = half*(w(i, j, k+1, itu1)+w(i, j, k, itu1))
! compute nu at the faces
          temp = w(i, j, k, irho)
          temp0 = rlv(i, j, k)/temp
          nud = (rlvd(i, j, k)-temp0*wd(i, j, k, irho))/temp
          nu = temp0
          temp0 = w(i, j, k-1, irho)
          temp = rlv(i, j, k-1)/temp0
          numd = half*((rlvd(i, j, k-1)-temp*wd(i, j, k-1, irho))/temp0+&
&           nud)
          num = half*(temp+nu)
          temp0 = w(i, j, k+1, irho)
          temp = rlv(i, j, k+1)/temp0
          nupd = half*((rlvd(i, j, k+1)-temp*wd(i, j, k+1, irho))/temp0+&
&           nud)
          nup = half*(temp+nu)
          temp0 = num + (one+rsacb2)*nutm
          cdmd = cb3inv*(ttm*(numd+(one+rsacb2)*nutmd)+temp0*ttmd)
          cdm = cb3inv*(temp0*ttm)
          temp0 = nup + (one+rsacb2)*nutp
          cdpd = cb3inv*(ttp*(nupd+(one+rsacb2)*nutpd)+temp0*ttpd)
          cdp = cb3inv*(temp0*ttp)
          if (cdm + cam .lt. zero) then
            c1m = zero
            c1md = 0.0_8
          else
            c1md = cdmd + camd
            c1m = cdm + cam
          end if
          if (cdp + cap .lt. zero) then
            c1p = zero
            c1pd = 0.0_8
          else
            c1pd = cdpd + capd
            c1p = cdp + cap
          end if
          c10d = c1md + c1pd
          c10 = c1m + c1p
! update the residual for this cell and store the possible
! coefficients for the matrix in b1, c1 and d1.
          temp0 = w(i, j, k-1, itu1)
          temp = w(i, j, k+1, itu1)
          temp1 = w(i, j, k, itu1)
          scratchd(i, j, k, idvt) = scratchd(i, j, k, idvt) + temp0*c1md&
&           + c1m*wd(i, j, k-1, itu1) + temp*c1pd + c1p*wd(i, j, k+1, &
&           itu1) - temp1*c10d - c10*wd(i, j, k, itu1)
          scratch(i, j, k, idvt) = scratch(i, j, k, idvt) + c1m*temp0 + &
&           c1p*temp - c10*temp1
        end do
      end do
    end do
!
!       viscous terms in j-direction.
!
    do k=2,kl
      do j=2,jl
        do i=2,il
! compute the metrics in eta-direction, i.e. along the
! line j = constant.
          temp1 = one/vol(i, j, k)
          volid = -(temp1*vold(i, j, k)/vol(i, j, k))
          voli = temp1
          temp1 = two/(vol(i, j, k)+vol(i, j-1, k))
          volmid = -(temp1*(vold(i, j, k)+vold(i, j-1, k))/(vol(i, j, k)&
&           +vol(i, j-1, k)))
          volmi = temp1
          temp1 = two/(vol(i, j, k)+vol(i, j+1, k))
          volpid = -(temp1*(vold(i, j, k)+vold(i, j+1, k))/(vol(i, j, k)&
&           +vol(i, j+1, k)))
          volpi = temp1
          temp1 = sj(i, j-1, k, 1)
          xmd = volmi*sjd(i, j-1, k, 1) + temp1*volmid
          xm = temp1*volmi
          temp1 = sj(i, j-1, k, 2)
          ymd = volmi*sjd(i, j-1, k, 2) + temp1*volmid
          ym = temp1*volmi
          temp1 = sj(i, j-1, k, 3)
          zmd = volmi*sjd(i, j-1, k, 3) + temp1*volmid
          zm = temp1*volmi
          temp1 = sj(i, j, k, 1)
          xpd = volpi*sjd(i, j, k, 1) + temp1*volpid
          xp = temp1*volpi
          temp1 = sj(i, j, k, 2)
          ypd = volpi*sjd(i, j, k, 2) + temp1*volpid
          yp = temp1*volpi
          temp1 = sj(i, j, k, 3)
          zpd = volpi*sjd(i, j, k, 3) + temp1*volpid
          zp = temp1*volpi
          temp1 = sj(i, j, k, 1) + sj(i, j-1, k, 1)
          xad = half*(voli*(sjd(i, j, k, 1)+sjd(i, j-1, k, 1))+temp1*&
&           volid)
          xa = half*(temp1*voli)
          temp1 = sj(i, j, k, 2) + sj(i, j-1, k, 2)
          yad = half*(voli*(sjd(i, j, k, 2)+sjd(i, j-1, k, 2))+temp1*&
&           volid)
          ya = half*(temp1*voli)
          temp1 = sj(i, j, k, 3) + sj(i, j-1, k, 3)
          zad = half*(voli*(sjd(i, j, k, 3)+sjd(i, j-1, k, 3))+temp1*&
&           volid)
          za = half*(temp1*voli)
          ttmd = xa*xmd + xm*xad + ya*ymd + ym*yad + za*zmd + zm*zad
          ttm = xm*xa + ym*ya + zm*za
          ttpd = xa*xpd + xp*xad + ya*ypd + yp*yad + za*zpd + zp*zad
          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.
          cnudd = -(rsacb2*cb3inv*wd(i, j, k, itu1))
          cnud = -(rsacb2*w(i, j, k, itu1)*cb3inv)
          camd = cnud*ttmd + ttm*cnudd
          cam = ttm*cnud
          capd = cnud*ttpd + ttp*cnudd
          cap = ttp*cnud
          nutmd = half*(wd(i, j-1, k, itu1)+wd(i, j, k, itu1))
          nutm = half*(w(i, j-1, k, itu1)+w(i, j, k, itu1))
          nutpd = half*(wd(i, j+1, k, itu1)+wd(i, j, k, itu1))
          nutp = half*(w(i, j+1, k, itu1)+w(i, j, k, itu1))
          temp1 = w(i, j, k, irho)
          temp0 = rlv(i, j, k)/temp1
          nud = (rlvd(i, j, k)-temp0*wd(i, j, k, irho))/temp1
          nu = temp0
          temp1 = w(i, j-1, k, irho)
          temp0 = rlv(i, j-1, k)/temp1
          numd = half*((rlvd(i, j-1, k)-temp0*wd(i, j-1, k, irho))/temp1&
&           +nud)
          num = half*(temp0+nu)
          temp1 = w(i, j+1, k, irho)
          temp0 = rlv(i, j+1, k)/temp1
          nupd = half*((rlvd(i, j+1, k)-temp0*wd(i, j+1, k, irho))/temp1&
&           +nud)
          nup = half*(temp0+nu)
          temp1 = num + (one+rsacb2)*nutm
          cdmd = cb3inv*(ttm*(numd+(one+rsacb2)*nutmd)+temp1*ttmd)
          cdm = cb3inv*(temp1*ttm)
          temp1 = nup + (one+rsacb2)*nutp
          cdpd = cb3inv*(ttp*(nupd+(one+rsacb2)*nutpd)+temp1*ttpd)
          cdp = cb3inv*(temp1*ttp)
          if (cdm + cam .lt. zero) then
            c1m = zero
            c1md = 0.0_8
          else
            c1md = cdmd + camd
            c1m = cdm + cam
          end if
          if (cdp + cap .lt. zero) then
            c1p = zero
            c1pd = 0.0_8
          else
            c1pd = cdpd + capd
            c1p = cdp + cap
          end if
          c10d = c1md + c1pd
          c10 = c1m + c1p
! update the residual for this cell and store the possible
! coefficients for the matrix in b1, c1 and d1.
          temp1 = w(i, j-1, k, itu1)
          temp0 = w(i, j+1, k, itu1)
          temp = w(i, j, k, itu1)
          scratchd(i, j, k, idvt) = scratchd(i, j, k, idvt) + temp1*c1md&
&           + c1m*wd(i, j-1, k, itu1) + temp0*c1pd + c1p*wd(i, j+1, k, &
&           itu1) - temp*c10d - c10*wd(i, j, k, itu1)
          scratch(i, j, k, idvt) = scratch(i, j, k, idvt) + c1m*temp1 + &
&           c1p*temp0 - c10*temp
        end do
      end do
    end do
!
!       viscous terms in i-direction.
!
    do k=2,kl
      do j=2,jl
        do i=2,il
! compute the metrics in xi-direction, i.e. along the
! line i = constant.
          temp1 = one/vol(i, j, k)
          volid = -(temp1*vold(i, j, k)/vol(i, j, k))
          voli = temp1
          temp1 = two/(vol(i, j, k)+vol(i-1, j, k))
          volmid = -(temp1*(vold(i, j, k)+vold(i-1, j, k))/(vol(i, j, k)&
&           +vol(i-1, j, k)))
          volmi = temp1
          temp1 = two/(vol(i, j, k)+vol(i+1, j, k))
          volpid = -(temp1*(vold(i, j, k)+vold(i+1, j, k))/(vol(i, j, k)&
&           +vol(i+1, j, k)))
          volpi = temp1
          temp1 = si(i-1, j, k, 1)
          xmd = volmi*sid(i-1, j, k, 1) + temp1*volmid
          xm = temp1*volmi
          temp1 = si(i-1, j, k, 2)
          ymd = volmi*sid(i-1, j, k, 2) + temp1*volmid
          ym = temp1*volmi
          temp1 = si(i-1, j, k, 3)
          zmd = volmi*sid(i-1, j, k, 3) + temp1*volmid
          zm = temp1*volmi
          temp1 = si(i, j, k, 1)
          xpd = volpi*sid(i, j, k, 1) + temp1*volpid
          xp = temp1*volpi
          temp1 = si(i, j, k, 2)
          ypd = volpi*sid(i, j, k, 2) + temp1*volpid
          yp = temp1*volpi
          temp1 = si(i, j, k, 3)
          zpd = volpi*sid(i, j, k, 3) + temp1*volpid
          zp = temp1*volpi
          temp1 = si(i, j, k, 1) + si(i-1, j, k, 1)
          xad = half*(voli*(sid(i, j, k, 1)+sid(i-1, j, k, 1))+temp1*&
&           volid)
          xa = half*(temp1*voli)
          temp1 = si(i, j, k, 2) + si(i-1, j, k, 2)
          yad = half*(voli*(sid(i, j, k, 2)+sid(i-1, j, k, 2))+temp1*&
&           volid)
          ya = half*(temp1*voli)
          temp1 = si(i, j, k, 3) + si(i-1, j, k, 3)
          zad = half*(voli*(sid(i, j, k, 3)+sid(i-1, j, k, 3))+temp1*&
&           volid)
          za = half*(temp1*voli)
          ttmd = xa*xmd + xm*xad + ya*ymd + ym*yad + za*zmd + zm*zad
          ttm = xm*xa + ym*ya + zm*za
          ttpd = xa*xpd + xp*xad + ya*ypd + yp*yad + za*zpd + zp*zad
          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.
          cnudd = -(rsacb2*cb3inv*wd(i, j, k, itu1))
          cnud = -(rsacb2*w(i, j, k, itu1)*cb3inv)
          camd = cnud*ttmd + ttm*cnudd
          cam = ttm*cnud
          capd = cnud*ttpd + ttp*cnudd
          cap = ttp*cnud
          nutmd = half*(wd(i-1, j, k, itu1)+wd(i, j, k, itu1))
          nutm = half*(w(i-1, j, k, itu1)+w(i, j, k, itu1))
          nutpd = half*(wd(i+1, j, k, itu1)+wd(i, j, k, itu1))
          nutp = half*(w(i+1, j, k, itu1)+w(i, j, k, itu1))
          temp1 = w(i, j, k, irho)
          temp0 = rlv(i, j, k)/temp1
          nud = (rlvd(i, j, k)-temp0*wd(i, j, k, irho))/temp1
          nu = temp0
          temp1 = w(i-1, j, k, irho)
          temp0 = rlv(i-1, j, k)/temp1
          numd = half*((rlvd(i-1, j, k)-temp0*wd(i-1, j, k, irho))/temp1&
&           +nud)
          num = half*(temp0+nu)
          temp1 = w(i+1, j, k, irho)
          temp0 = rlv(i+1, j, k)/temp1
          nupd = half*((rlvd(i+1, j, k)-temp0*wd(i+1, j, k, irho))/temp1&
&           +nud)
          nup = half*(temp0+nu)
          temp1 = num + (one+rsacb2)*nutm
          cdmd = cb3inv*(ttm*(numd+(one+rsacb2)*nutmd)+temp1*ttmd)
          cdm = cb3inv*(temp1*ttm)
          temp1 = nup + (one+rsacb2)*nutp
          cdpd = cb3inv*(ttp*(nupd+(one+rsacb2)*nutpd)+temp1*ttpd)
          cdp = cb3inv*(temp1*ttp)
          if (cdm + cam .lt. zero) then
            c1m = zero
            c1md = 0.0_8
          else
            c1md = cdmd + camd
            c1m = cdm + cam
          end if
          if (cdp + cap .lt. zero) then
            c1p = zero
            c1pd = 0.0_8
          else
            c1pd = cdpd + capd
            c1p = cdp + cap
          end if
          c10d = c1md + c1pd
          c10 = c1m + c1p
! update the residual for this cell and store the possible
! coefficients for the matrix in b1, c1 and d1.
          temp1 = w(i-1, j, k, itu1)
          temp0 = w(i+1, j, k, itu1)
          temp = w(i, j, k, itu1)
          scratchd(i, j, k, idvt) = scratchd(i, j, k, idvt) + temp1*c1md&
&           + c1m*wd(i-1, j, k, itu1) + temp0*c1pd + c1p*wd(i+1, j, k, &
&           itu1) - temp*c10d - c10*wd(i, j, k, itu1)
          scratch(i, j, k, idvt) = scratch(i, j, k, idvt) + c1m*temp1 + &
&           c1p*temp0 - c10*temp
        end do
      end do
    end do
  end subroutine saviscous_d
 
  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 max
! set model constants
    cv13 = rsacv1**3
    kar2inv = one/rsak**2
    cw36 = rsacw3**6
    cb3inv = one/rsacb3
!
!       viscous terms in k-direction.
!
    do k=2,kl
      do j=2,jl
        do i=2,il
! 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
      end do
    end do
!
!       viscous terms in j-direction.
!
    do k=2,kl
      do j=2,jl
        do i=2,il
! 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
      end do
    end do
!
!       viscous terms in i-direction.
!
    do k=2,kl
      do j=2,jl
        do i=2,il
! 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 do
    end do
  end subroutine saviscous
 
!  differentiation of saresscale in forward (tangent) mode (with options i4 dr8 r8):
!   variations   of useful results: *dw
!   with respect to varying inputs: *dw *scratch
!   rw status of diff variables: *dw:in-out *scratch:in
!   plus diff mem management of: dw:in scratch:in
  subroutine saresscale_d()
!
!  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 real
    intrinsic max
    real(realtype) :: x1
    real(kind=realtype) :: temp
    do k=2,kl
      do j=2,jl
        do i=2,il
          x1 = real(iblank(i, j, k), realtype)
          if (x1 .lt. zero) then
            rblank = zero
          else
            rblank = x1
          end if
          temp = volref(i, j, k)*rblank
          dwd(i, j, k, itu1) = -(temp*scratchd(i, j, k, idvt))
          dw(i, j, k, itu1) = -(temp*scratch(i, j, k, idvt))
        end do
      end do
    end do
  end subroutine saresscale_d
 
  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 real
    intrinsic max
    real(realtype) :: x1
    do k=2,kl
      do j=2,jl
        do i=2,il
          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 do
    end do
  end subroutine saresscale
 
end module sa_d