pkg/mom_common/mom_calc_visc.F

f5bcbacdce Bayl*0001 #include "MOM_COMMON_OPTIONS.h"
14bb46b4fa Jean*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
f5bcbacdce Bayl*0005
f0501c53d1 Jean*0006 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0007 CBOP
                0008 C     !ROUTINE: MOM_CALC_VISC
b8d51be678 Bayl*0009
f0501c53d1 Jean*0010 C     !INTERFACE:
f5bcbacdce Bayl*0011       SUBROUTINE MOM_CALC_VISC(
                0012      I        bi,bj,k,
                0013      O        viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
f59d76b0dd Ed D*0014      I        hDiv,vort3,tension,strain,stretching,KE,hFacZ,
f5bcbacdce Bayl*0015      I        myThid)
                0016
f0501c53d1 Jean*0017 C     !DESCRIPTION:
b8d51be678 Bayl*0018 C     Calculate horizontal viscosities (L is typical grid width)
                0019 C     harmonic viscosity=
                0020 C       viscAh (or viscAhD on div pts and viscAhZ on zeta pts)
                0021 C       +0.25*L**2*viscAhGrid/deltaT
20e7e2b61c Bayl*0022 C       +sqrt((viscC2leith/pi)**6*grad(Vort3)**2
                0023 C             +(viscC2leithD/pi)**6*grad(hDiv)**2)*L**3
b8d51be678 Bayl*0024 C       +(viscC2smag/pi)**2*L**2*sqrt(Tension**2+Strain**2)
                0025 C
                0026 C     biharmonic viscosity=
                0027 C       viscA4 (or viscA4D on div pts and viscA4Z on zeta pts)
                0028 C       +0.25*0.125*L**4*viscA4Grid/deltaT (approx)
20e7e2b61c Bayl*0029 C       +0.125*L**5*sqrt((viscC4leith/pi)**6*grad(Vort3)**2
                0030 C                        +(viscC4leithD/pi)**6*grad(hDiv)**2)
b8d51be678 Bayl*0031 C       +0.125*L**4*(viscC4smag/pi)**2*sqrt(Tension**2+Strain**2)
                0032 C
                0033 C     Note that often 0.125*L**2 is the scale between harmonic and
                0034 C     biharmonic (see Griffies and Hallberg (2000))
                0035 C     This allows the same value of the coefficient to be used
                0036 C     for roughly similar results with biharmonic and harmonic
                0037 C
                0038 C     LIMITERS -- limit min and max values of viscosities
002795c988 Jean*0039 C     viscAhReMax is min value for grid point harmonic Reynolds num
                0040 C      harmonic viscosity>sqrt(2*KE)*L/viscAhReMax
b8d51be678 Bayl*0041 C
002795c988 Jean*0042 C     viscA4ReMax is min value for grid point biharmonic Reynolds num
                0043 C      biharmonic viscosity>sqrt(2*KE)*L**3/8/viscA4ReMax
b8d51be678 Bayl*0044 C
                0045 C     viscAhgridmax is CFL stability limiter for harmonic viscosity
                0046 C      harmonic viscosity<0.25*viscAhgridmax*L**2/deltaT
                0047 C
                0048 C     viscA4gridmax is CFL stability limiter for biharmonic viscosity
                0049 C      biharmonic viscosity<viscA4gridmax*L**4/32/deltaT (approx)
                0050 C
                0051 C     viscAhgridmin and viscA4gridmin are lower limits for viscosity:
420f2398b8 Chri*0052 C       harmonic viscosity>0.25*viscAhgridmin*L**2/deltaT
                0053 C       biharmonic viscosity>viscA4gridmin*L**4/32/deltaT (approx)
                0054
b8d51be678 Bayl*0055 C     RECOMMENDED VALUES
a8b99ae3ea Bayl*0056 C     viscC2Leith=1-3
                0057 C     viscC2LeithD=1-3
f59d76b0dd Ed D*0058 C     viscC2LeithQG=1
a8b99ae3ea Bayl*0059 C     viscC4Leith=1-3
                0060 C     viscC4LeithD=1.5-3
002795c988 Jean*0061 C     viscC2smag=2.2-4 (Griffies and Hallberg,2000)
b8d51be678 Bayl*0062 C               0.2-0.9 (Smagorinsky,1993)
002795c988 Jean*0063 C     viscC4smag=2.2-4 (Griffies and Hallberg,2000)
                0064 C     viscAhReMax>=1, (<2 suppresses a computational mode)
                0065 C     viscA4ReMax>=1, (<2 suppresses a computational mode)
b8d51be678 Bayl*0066 C     viscAhgridmax=1
                0067 C     viscA4gridmax=1
                0068 C     viscAhgrid<1
                0069 C     viscA4grid<1
                0070 C     viscAhgridmin<<1
                0071 C     viscA4gridmin<<1
f5bcbacdce Bayl*0072
f0501c53d1 Jean*0073 C     !USES:
                0074       IMPLICIT NONE
                0075
f5bcbacdce Bayl*0076 C     == Global variables ==
                0077 #include "SIZE.h"
                0078 #include "GRID.h"
                0079 #include "EEPARAMS.h"
                0080 #include "PARAMS.h"
f0501c53d1 Jean*0081 #include "MOM_VISC.h"
7c50f07931 Mart*0082 #ifdef ALLOW_AUTODIFF_TAMC
cc254a4876 Patr*0083 #include "tamc.h"
7c50f07931 Mart*0084 #endif /* ALLOW_AUTODIFF_TAMC */
f5bcbacdce Bayl*0085
f0501c53d1 Jean*0086 C     !INPUT/OUTPUT PARAMETERS:
                0087 C     myThid               :: my thread Id number
f5bcbacdce Bayl*0088       INTEGER bi,bj,k
                0089       _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0090       _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0091       _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0092       _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0093       _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0094       _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0095       _RL tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0096       _RL strain(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
f59d76b0dd Ed D*0097       _RL stretching(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
f5bcbacdce Bayl*0098       _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0099       _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0100       INTEGER myThid
f0501c53d1 Jean*0101 CEOP
f5bcbacdce Bayl*0102
f0501c53d1 Jean*0103 C     !LOCAL VARIABLES:
55479d97c4 Jean*0104       INTEGER i,j
edc5c0b6b1 Jean*0105 #ifdef ALLOW_NONHYDROSTATIC
f0501c53d1 Jean*0106       _RL shiftAh, shiftA4
edc5c0b6b1 Jean*0107 #endif
0f82b218cd Jean*0108 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0109 C     kkey   :: tape key (depends on levels and tiles)
                0110 C     ijkkey :: tape key (depends on i,j-indices, levels, and tiles)
                0111       INTEGER kkey, ijkkey
0f82b218cd Jean*0112 #endif
b8d51be678 Bayl*0113       _RL smag2fac, smag4fac
20e7e2b61c Bayl*0114       _RL leith2fac, leith4fac
                0115       _RL leithD2fac, leithD4fac
f59d76b0dd Ed D*0116       _RL leithQG2fac
1e9e8678dd Bayl*0117       _RL viscAhRe_max, viscA4Re_max
f59d76b0dd Ed D*0118       _RL Alin, grdVrt, grdDiv, keZpt
220a9e88b5 Jean*0119       _RL deepFac3, deepFac4
f0501c53d1 Jean*0120       _RL L2, L3, L5, L2rdt, L4rdt, recip_dt
b8d51be678 Bayl*0121       _RL Uscl,U4scl
e47a30e113 Jean*0122       _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0123       _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8def1fc7b9 Jean*0124       _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0125       _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
b8d51be678 Bayl*0126       _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0127       _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0128       _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0129       _RL viscA4_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0130       _RL viscAh_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0131       _RL viscAh_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0132       _RL viscA4_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0133       _RL viscA4_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0134       _RL viscAh_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0135       _RL viscAh_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0136       _RL viscA4_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0137       _RL viscA4_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0138       _RL viscAh_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0139       _RL viscAh_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0140       _RL viscA4_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0141       _RL viscA4_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
06244a5e4f Jean*0142 #ifdef ALLOW_LEITH_QG
f59d76b0dd Ed D*0143       _RL viscAh_ZLthQG(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0144       _RL viscAh_DLthQG(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
aa6b2555c8 Jean*0145       _RL sqargQG
06244a5e4f Jean*0146 #endif
b8d51be678 Bayl*0147       _RL viscAh_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0148       _RL viscAh_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0149       _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0150       _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
aa6b2555c8 Jean*0151       _RL sqargAh, sqargA4, sqargAhD, sqargA4D, sqargSmag
f59d76b0dd Ed D*0152       LOGICAL calcLeith, calcSmag, calcLeithQG
f5bcbacdce Bayl*0153
cc254a4876 Patr*0154 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0155       kkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0156       kkey = k  + (kkey-1)*Nr
cc254a4876 Patr*0157 #endif /* ALLOW_AUTODIFF_TAMC */
                0158
002795c988 Jean*0159 C--   Set flags which are used in this S/R and elsewhere :
f0501c53d1 Jean*0160 C     useVariableVisc, useHarmonicVisc and useBiharmonicVisc
                0161 C     are now set early on (in S/R SET_PARAMS)
f5bcbacdce Bayl*0162
ba2530fec3 Jean*0163 c     IF ( useVariableVisc ) THEN
002795c988 Jean*0164 C---- variable viscosity :
                0165
f0501c53d1 Jean*0166        recip_dt = 1. _d 0
14bb46b4fa Jean*0167        IF ( deltaTMom.NE.0. ) recip_dt = 1. _d 0/deltaTMom
220a9e88b5 Jean*0168        deepFac3 = deepFac2C(k)*deepFacC(k)
                0169        deepFac4 = deepFac2C(k)*deepFac2C(k)
f0501c53d1 Jean*0170
                0171        IF ( useHarmonicVisc .AND. viscAhReMax.NE.0. ) THEN
002795c988 Jean*0172         viscAhRe_max=SQRT(2. _d 0)/viscAhReMax
                0173        ELSE
                0174         viscAhRe_max=0. _d 0
                0175        ENDIF
                0176
f0501c53d1 Jean*0177        IF ( useBiharmonicVisc .AND. viscA4ReMax.NE.0. ) THEN
002795c988 Jean*0178         viscA4Re_max=0.125 _d 0*SQRT(2. _d 0)/viscA4ReMax
                0179        ELSE
                0180         viscA4Re_max=0. _d 0
                0181        ENDIF
b8d51be678 Bayl*0182
f59d76b0dd Ed D*0183        calcLeithQG = (viscC2LeithQG.NE.zeroRL)
002795c988 Jean*0184        calcLeith=
b8d51be678 Bayl*0185      &      (viscC2leith.NE.0.)
                0186      &  .OR.(viscC2leithD.NE.0.)
                0187      &  .OR.(viscC4leith.NE.0.)
                0188      &  .OR.(viscC4leithD.NE.0.)
f59d76b0dd Ed D*0189      &  .OR. calcLeithQG
b8d51be678 Bayl*0190
002795c988 Jean*0191        calcSmag=
b8d51be678 Bayl*0192      &      (viscC2smag.NE.0.)
                0193      &  .OR.(viscC4smag.NE.0.)
                0194
002795c988 Jean*0195        IF (calcSmag) THEN
b8d51be678 Bayl*0196         smag2fac=(viscC2smag/pi)**2
bb1bd5219d Jean*0197         smag4fac=0.125 _d 0*(viscC4smag/pi)**2
002795c988 Jean*0198        ELSE
bb1bd5219d Jean*0199         smag2fac=0. _d 0
                0200         smag4fac=0. _d 0
002795c988 Jean*0201        ENDIF
f5bcbacdce Bayl*0202
002795c988 Jean*0203        IF (calcLeith) THEN
20e7e2b61c Bayl*0204         IF (useFullLeith) THEN
f59d76b0dd Ed D*0205 C       Uses correct calculation for gradients, but might not work on cube sphere
e2b2bf4d90 Bayl*0206          leith2fac =(viscC2leith /pi)**6
20e7e2b61c Bayl*0207          leithD2fac=(viscC2leithD/pi)**6
f59d76b0dd Ed D*0208          leithQG2fac = (viscC2LeithQG/pi)**6
e2b2bf4d90 Bayl*0209          leith4fac =0.015625 _d 0*(viscC4leith /pi)**6
20e7e2b61c Bayl*0210          leithD4fac=0.015625 _d 0*(viscC4leithD/pi)**6
                0211         ELSE
f59d76b0dd Ed D*0212 C       Uses approximate gradients, but works on cube sphere. No reason to use this
                0213 C        unless `useFullLeith` fails for your setup.
e2b2bf4d90 Bayl*0214          leith2fac =(viscC2leith /pi)**3
20e7e2b61c Bayl*0215          leithD2fac=(viscC2leithD/pi)**3
f59d76b0dd Ed D*0216          leithQG2fac = (viscC2LeithQG/pi)**3
e2b2bf4d90 Bayl*0217          leith4fac =0.125 _d 0*(viscC4leith /pi)**3
                0218          leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3
20e7e2b61c Bayl*0219         ENDIF
002795c988 Jean*0220        ELSE
20e7e2b61c Bayl*0221         leith2fac=0. _d 0
                0222         leith4fac=0. _d 0
f59d76b0dd Ed D*0223         leithQG2fac=0. _d 0
20e7e2b61c Bayl*0224         leithD2fac=0. _d 0
                0225         leithD4fac=0. _d 0
002795c988 Jean*0226        ENDIF
20e7e2b61c Bayl*0227
f0501c53d1 Jean*0228        DO j=1-OLy,sNy+OLy
                0229         DO i=1-OLx,sNx+OLx
ba2530fec3 Jean*0230 C-    viscosity arrays have been initialised everywhere before calling this S/R
                0231 c         viscAh_D(i,j) = viscAhD
                0232 c         viscAh_Z(i,j) = viscAhZ
                0233 c         viscA4_D(i,j) = viscA4D
                0234 c         viscA4_Z(i,j) = viscA4Z
f0501c53d1 Jean*0235
06244a5e4f Jean*0236           viscAh_DLth(i,j) = 0. _d 0
812e613baf Patr*0237           viscAh_ZLth(i,j) = 0. _d 0
06244a5e4f Jean*0238           viscA4_DLth(i,j) = 0. _d 0
812e613baf Patr*0239           viscA4_ZLth(i,j) = 0. _d 0
06244a5e4f Jean*0240           viscAh_DLthD(i,j)= 0. _d 0
812e613baf Patr*0241           viscAh_ZLthD(i,j)= 0. _d 0
06244a5e4f Jean*0242           viscA4_DLthD(i,j)= 0. _d 0
812e613baf Patr*0243           viscA4_ZLthD(i,j)= 0. _d 0
06244a5e4f Jean*0244 #ifdef ALLOW_LEITH_QG
                0245           viscAh_DLthQG(i,j) = 0. _d 0
                0246           viscAh_ZLthQG(i,j) = 0. _d 0
                0247 #endif
                0248
                0249           viscAh_DSmg(i,j) = 0. _d 0
                0250           viscAh_ZSmg(i,j) = 0. _d 0
                0251           viscA4_DSmg(i,j) = 0. _d 0
                0252           viscA4_ZSmg(i,j) = 0. _d 0
812e613baf Patr*0253         ENDDO
002795c988 Jean*0254        ENDDO
e47a30e113 Jean*0255
f59d76b0dd Ed D*0256 C-    Initialise to zero gradient of vorticity and divergence:
f0501c53d1 Jean*0257        DO j=1-OLy,sNy+OLy
                0258         DO i=1-OLx,sNx+OLx
e47a30e113 Jean*0259           divDx(i,j) = 0.
                0260           divDy(i,j) = 0.
8def1fc7b9 Jean*0261           vrtDx(i,j) = 0.
                0262           vrtDy(i,j) = 0.
e47a30e113 Jean*0263         ENDDO
                0264        ENDDO
8def1fc7b9 Jean*0265
ba2530fec3 Jean*0266        IF ( calcLeith ) THEN
                0267 C--   horizontal gradient of horizontal divergence:
e47a30e113 Jean*0268 C-       gradient in x direction:
                0269          IF (useCubedSphereExchange) THEN
                0270 C        to compute d/dx(hDiv), fill corners with appropriate values:
93e3461d85 Jean*0271            CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
1891130b05 Jean*0272      &                                hDiv, bi,bj, myThid )
e47a30e113 Jean*0273          ENDIF
f0501c53d1 Jean*0274          DO j=2-OLy,sNy+OLy-1
                0275           DO i=2-OLx,sNx+OLx-1
220a9e88b5 Jean*0276             divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))
                0277      &                  *recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
e47a30e113 Jean*0278           ENDDO
                0279          ENDDO
                0280
                0281 C-       gradient in y direction:
                0282          IF (useCubedSphereExchange) THEN
                0283 C        to compute d/dy(hDiv), fill corners with appropriate values:
93e3461d85 Jean*0284            CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
1891130b05 Jean*0285      &                                hDiv, bi,bj, myThid )
e47a30e113 Jean*0286          ENDIF
f0501c53d1 Jean*0287          DO j=2-OLy,sNy+OLy-1
                0288           DO i=2-OLx,sNx+OLx-1
220a9e88b5 Jean*0289             divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))
                0290      &                  *recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
e47a30e113 Jean*0291           ENDDO
                0292          ENDDO
8def1fc7b9 Jean*0293
ba2530fec3 Jean*0294 C--   horizontal gradient of vertical vorticity:
8def1fc7b9 Jean*0295 C-       gradient in x direction:
f0501c53d1 Jean*0296          DO j=2-OLy,sNy+OLy
                0297           DO i=2-OLx,sNx+OLx-1
8def1fc7b9 Jean*0298             vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
220a9e88b5 Jean*0299      &                  *recip_dxG(i,j,bi,bj)*recip_deepFacC(k)
8def1fc7b9 Jean*0300      &                  *maskS(i,j,k,bi,bj)
55479d97c4 Jean*0301 #ifdef ALLOW_OBCS
                0302      &                  *maskInS(i,j,bi,bj)
                0303 #endif
8def1fc7b9 Jean*0304           ENDDO
                0305          ENDDO
                0306 C-       gradient in y direction:
f0501c53d1 Jean*0307          DO j=2-OLy,sNy+OLy-1
                0308           DO i=2-OLx,sNx+OLx
8def1fc7b9 Jean*0309             vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
220a9e88b5 Jean*0310      &                  *recip_dyG(i,j,bi,bj)*recip_deepFacC(k)
8def1fc7b9 Jean*0311      &                  *maskW(i,j,k,bi,bj)
55479d97c4 Jean*0312 #ifdef ALLOW_OBCS
                0313      &                  *maskInW(i,j,bi,bj)
                0314 #endif
8def1fc7b9 Jean*0315           ENDDO
                0316          ENDDO
                0317
06244a5e4f Jean*0318 #ifdef ALLOW_LEITH_QG
                0319         IF ( calcLeithQG ) THEN
f59d76b0dd Ed D*0320 C      horizontal gradient of vorticity and vortex stretching:
                0321 C        In the case of using QG Leith, we want to add a term
                0322 C        before calculating vector magnitude, so add to the
                0323 C        values just calculated.
                0324 C        gradient in x direction:
                0325          DO j=2-OLy,sNy+OLy
                0326           DO i=2-OLx,sNx+OLx-1
                0327 C        Average d/dx of stretching onto V-points to match vrtDX
                0328             vrtDx(i,j) = vrtDx(i,j)
                0329      &                 + halfRL*halfRL*
                0330      &                   ( (stretching(i+1,j)-stretching(i,j))
                0331      &                     *recip_dxC(i+1,j,bi,bj)*recip_deepFacC(k)
                0332      &                   + (stretching(i,j)-stretching(i-1,j))
                0333      &                     *recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
                0334      &                   + (stretching(i+1,j-1)-stretching(i,j-1))
                0335      &                     *recip_dxC(i,j-1,bi,bj)*recip_deepFacC(k)
                0336      &                   + (stretching(i,j-1)-stretching(i-1,j-1))
                0337      &                     *recip_dxC(i-1,j-1,bi,bj)*recip_deepFacC(k)
                0338      &                   )*maskS(i,j,k,bi,bj)
                0339 #ifdef ALLOW_OBCS
                0340      &                    *maskInS(i,j,bi,bj)
                0341 #endif
                0342           ENDDO
                0343          ENDDO
                0344 C-       gradient in y direction:
                0345          DO j=2-OLy,sNy+OLy-1
                0346           DO i=2-OLx,sNx+OLx
                0347 C        Average d/dy of stretching onto U-points to match vrtDy
                0348             vrtDy(i,j) = vrtDy(i,j)
                0349      &                 + halfRL*halfRL*
                0350      &                   ( (stretching(i,j+1)-stretching(i,j))
                0351      &                     *recip_dyC(i,j+1,bi,bj)*recip_deepFacC(k)
                0352      &                   + (stretching(i,j)-stretching(i,j-1))
                0353      &                     *recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
                0354      &                   + (stretching(i-1,j+1)-stretching(i-1,j))
                0355      &                     *recip_dyC(i-1,j+1,bi,bj)*recip_deepFacC(k)
                0356      &                   + (stretching(i-1,j)-stretching(i-1,j-1))
                0357      &                     *recip_dyC(i-1,j,bi,bj)*recip_deepFacC(k)
                0358      &                   )*maskW(i,j,k,bi,bj)
                0359 #ifdef ALLOW_OBCS
                0360      &                    *maskInW(i,j,bi,bj)
                0361 #endif
                0362           ENDDO
                0363          ENDDO
                0364 C      end if calcLeithQG
                0365         ENDIF
06244a5e4f Jean*0366 #endif /* ALLOW_LEITH_QG */
f59d76b0dd Ed D*0367
ba2530fec3 Jean*0368 C--   end if calcLeith
e47a30e113 Jean*0369        ENDIF
                0370
f0501c53d1 Jean*0371        DO j=2-OLy,sNy+OLy-1
                0372         DO i=2-OLx,sNx+OLx-1
f5bcbacdce Bayl*0373 CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
b8d51be678 Bayl*0374
14bb46b4fa Jean*0375 #ifdef ALLOW_AUTODIFF_TAMC
ec3bfd0c34 Patr*0376 # ifndef AUTODIFF_DISABLE_LEITH
edb6656069 Mart*0377          ijkkey = i+OLx + (j+OLy-1)*(sNx+2*OLx)
                0378      &                   + (kkey-1)*(sNx+2*OLx)*(sNy+2*OLy)
                0379 CADJ STORE viscA4_ZSmg(i,j)=comlev1_mom_ijk_loop,key=ijkkey,byte=isbyte
                0380 CADJ STORE viscAh_ZSmg(i,j)=comlev1_mom_ijk_loop,key=ijkkey,byte=isbyte
ec3bfd0c34 Patr*0381 # endif
14bb46b4fa Jean*0382 #endif /* ALLOW_AUTODIFF_TAMC */
ec3bfd0c34 Patr*0383
002795c988 Jean*0384 C These are (powers of) length scales
220a9e88b5 Jean*0385          L2 = L2_D(i,j,bi,bj)*deepFac2C(k)
0c0d21fb5c Davi*0386          L2rdt = 0.25 _d 0*recip_dt*L2
220a9e88b5 Jean*0387          L3 = L3_D(i,j,bi,bj)*deepFac3
                0388          L4rdt = L4rdt_D(i,j,bi,bj)*deepFac4
0c0d21fb5c Davi*0389          L5 = (L2*L3)
b8d51be678 Bayl*0390
5c4dff8b20 Mart*0391 #ifndef AUTODIFF_DISABLE_REYNOLDS_SCALE
b8d51be678 Bayl*0392 C Velocity Reynolds Scale
8d4c48ddd0 Jean*0393          IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
002795c988 Jean*0394            Uscl=SQRT(KE(i,j)*L2)*viscAhRe_max
8d4c48ddd0 Jean*0395          ELSE
                0396            Uscl=0.
                0397          ENDIF
                0398          IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
002795c988 Jean*0399            U4scl=SQRT(KE(i,j))*L3*viscA4Re_max
8d4c48ddd0 Jean*0400          ELSE
                0401            U4scl=0.
                0402          ENDIF
5c4dff8b20 Mart*0403 #endif /* ndef AUTODIFF_DISABLE_REYNOLDS_SCALE */
b8d51be678 Bayl*0404
cc254a4876 Patr*0405 #ifndef AUTODIFF_DISABLE_LEITH
002795c988 Jean*0406          IF (useFullLeith.AND.calcLeith) THEN
f5bcbacdce Bayl*0407 C This is the vector magnitude of the vorticity gradient squared
8def1fc7b9 Jean*0408           grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
                0409      &                        + vrtDx(i,j)*vrtDx(i,j) )
                0410      &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
                0411      &                        + vrtDy(i,j)*vrtDy(i,j) )  )
f5bcbacdce Bayl*0412
                0413 C This is the vector magnitude of grad (div.v) squared
                0414 C Using it in Leith serves to damp instabilities in w.
e47a30e113 Jean*0415           grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j)
                0416      &                        + divDx(i,j)*divDx(i,j) )
                0417      &                     + (divDy(i,j+1)*divDy(i,j+1)
                0418      &                        + divDy(i,j)*divDy(i,j) )  )
b8d51be678 Bayl*0419
70ceddf25f Timo*0420           sqargAh  = leith2fac*grdVrt+leithD2fac*grdDiv
                0421           sqargA4  = leith4fac*grdVrt+leithD4fac*grdDiv
                0422           sqargAhD = leithD2fac*grdDiv
                0423           sqargA4D = leithD4fac*grdDiv
06244a5e4f Jean*0424 #ifdef ALLOW_LEITH_QG
70ceddf25f Timo*0425           sqargQG  = leithQG2fac*(grdVrt+grdDiv)
                0426 #endif
                0427
                0428 #ifdef ALLOW_AUTODIFF
                0429 C Avoid derivative of SQRT(0)
                0430           IF (sqargAh .GT.0. _d 0) sqargAh = SQRT(sqargAh)
                0431           IF (sqargA4 .GT.0. _d 0) sqargA4 = SQRT(sqargA4)
                0432           IF (sqargAhD .GT.0. _d 0) sqargAhD = SQRT(sqargAhD)
                0433           IF (sqargA4D .GT.0. _d 0) sqargA4D = SQRT(sqargA4D)
                0434 # ifdef ALLOW_LEITH_QG
                0435           IF (sqargQG .GT.0. _d 0) sqargQG = SQRT(sqargQG)
                0436 # endif
                0437 #else /* ALLOW_AUTODIFF */
                0438           sqargAh = SQRT(sqargAh)
                0439           sqargA4 = SQRT(sqargA4)
                0440           sqargAhD = SQRT(sqargAhD)
                0441           sqargA4D = SQRT(sqargA4D)
                0442 # ifdef ALLOW_LEITH_QG
                0443           sqargQG = SQRT(sqargQG)
                0444 # endif
                0445 #endif /* ALLOW_AUTODIFF */
                0446           viscAh_DLth(i,j) = sqargAh  * L3
                0447           viscA4_DLth(i,j) = sqargA4  * L5
                0448           viscAh_DLthd(i,j)= sqargAhD * L3
                0449           viscA4_DLthd(i,j)= sqargA4D * L5
                0450 #ifdef ALLOW_LEITH_QG
                0451           viscAh_DLthQG(i,j)=sqargQG  * L3
06244a5e4f Jean*0452 #endif
f59d76b0dd Ed D*0453
002795c988 Jean*0454          ELSEIF (calcLeith) THEN
f0501c53d1 Jean*0455 C but this approximation will work on cube (and differs by as much as 4X)
002795c988 Jean*0456           grdVrt=MAX( ABS(vrtDx(i,j+1)), ABS(vrtDx(i,j)) )
                0457           grdVrt=MAX( grdVrt, ABS(vrtDy(i+1,j)) )
                0458           grdVrt=MAX( grdVrt, ABS(vrtDy(i,j))   )
b8d51be678 Bayl*0459
f0501c53d1 Jean*0460 C This approximation is good to the same order as above...
002795c988 Jean*0461           grdDiv=MAX( ABS(divDx(i+1,j)), ABS(divDx(i,j)) )
                0462           grdDiv=MAX( grdDiv, ABS(divDy(i,j+1)) )
                0463           grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
f5bcbacdce Bayl*0464
06244a5e4f Jean*0465           viscAh_DLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
                0466           viscA4_DLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
                0467           viscAh_DLthD(i,j)=((leithD2fac*grdDiv))*L3
                0468           viscA4_DLthD(i,j)=((leithD4fac*grdDiv))*L5
                0469 #ifdef ALLOW_LEITH_QG
f59d76b0dd Ed D*0470           viscAh_DLthQG(i,j)=leithQG2fac*(grdVrt + grdDiv)*L3
06244a5e4f Jean*0471 #endif
f59d76b0dd Ed D*0472
f5bcbacdce Bayl*0473          ELSE
06244a5e4f Jean*0474           viscAh_DLth(i,j)=0. _d 0
                0475           viscA4_DLth(i,j)=0. _d 0
                0476           viscAh_DLthD(i,j)=0. _d 0
                0477           viscA4_DLthD(i,j)=0. _d 0
                0478 #ifdef ALLOW_LEITH_QG
f59d76b0dd Ed D*0479           viscAh_DLthQG(i,j)=0. _d 0
06244a5e4f Jean*0480 #endif
f5bcbacdce Bayl*0481          ENDIF
                0482
002795c988 Jean*0483          IF (calcSmag) THEN
70ceddf25f Timo*0484           sqargSmag = tension(i,j)**2
bb1bd5219d Jean*0485      &       +0.25 _d 0*(strain(i+1, j )**2+strain( i ,j+1)**2
70ceddf25f Timo*0486      &                  +strain(i  , j )**2+strain(i+1,j+1)**2)
                0487 #ifdef ALLOW_AUTODIFF
                0488 C Avoid derivative of SQRT(0)
                0489           IF (sqargSmag.GT.0. _d 0) sqargSmag = SQRT(sqargSmag)
                0490 #else
                0491           sqargSmag = SQRT(sqargSmag)
                0492 #endif
                0493           viscAh_DSmg(i,j)=L2*sqargSmag
b8d51be678 Bayl*0494           viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
                0495           viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)
f5bcbacdce Bayl*0496          ELSE
bb1bd5219d Jean*0497           viscAh_DSmg(i,j)=0. _d 0
                0498           viscA4_DSmg(i,j)=0. _d 0
f5bcbacdce Bayl*0499          ENDIF
cc254a4876 Patr*0500 #endif /* AUTODIFF_DISABLE_LEITH */
f5bcbacdce Bayl*0501
                0502 C  Harmonic on Div.u points
b8d51be678 Bayl*0503          Alin=viscAhD+viscAhGrid*L2rdt
06244a5e4f Jean*0504      &         + viscAh_DLth(i,j)+viscAh_DSmg(i,j)
                0505 #ifdef ALLOW_LEITH_QG
                0506      &         + viscAh_DLthQG(i,j)
                0507 #endif
d062b1342f Gael*0508 #ifdef ALLOW_3D_VISCAH
06244a5e4f Jean*0509      &         + viscAhDfld(i,j,k,bi,bj)
4240547d2d Mart*0510 # ifdef AUTODIFF_ALLOW_VISCFACADJ
4c01acf062 Gael*0511      &          *viscFacAdj
4240547d2d Mart*0512 # endif /* AUTODIFF_ALLOW_VISCFACADJ */
                0513 #endif /* ALLOW_3D_VISCAH */
002795c988 Jean*0514          viscAh_DMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
                0515          viscAh_D(i,j)=MAX(viscAh_DMin(i,j),Alin)
                0516          viscAh_DMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
                0517          viscAh_D(i,j)=MIN(viscAh_DMax(i,j),viscAh_D(i,j))
f5bcbacdce Bayl*0518
                0519 C  BiHarmonic on Div.u points
b8d51be678 Bayl*0520          Alin=viscA4D+viscA4Grid*L4rdt
06244a5e4f Jean*0521      &         + viscA4_DLth(i,j)+viscA4_DSmg(i,j)
d062b1342f Gael*0522 #ifdef ALLOW_3D_VISCA4
06244a5e4f Jean*0523      &         + viscA4Dfld(i,j,k,bi,bj)
4240547d2d Mart*0524 # ifdef AUTODIFF_ALLOW_VISCFACADJ
4c01acf062 Gael*0525      &          *viscFacAdj
4240547d2d Mart*0526 # endif /* AUTODIFF_ALLOW_VISCFACADJ */
                0527 #endif /* ALLOW_3D_VISCA4 */
002795c988 Jean*0528          viscA4_DMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
                0529          viscA4_D(i,j)=MAX(viscA4_DMin(i,j),Alin)
                0530          viscA4_DMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
                0531          viscA4_D(i,j)=MIN(viscA4_DMax(i,j),viscA4_D(i,j))
f5bcbacdce Bayl*0532
                0533 CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
002795c988 Jean*0534 C These are (powers of) length scales
220a9e88b5 Jean*0535          L2 = L2_Z(i,j,bi,bj)*deepFac2C(k)
0c0d21fb5c Davi*0536          L2rdt = 0.25 _d 0*recip_dt*L2
220a9e88b5 Jean*0537          L3 = L3_Z(i,j,bi,bj)*deepFac3
                0538          L4rdt = L4rdt_Z(i,j,bi,bj)*deepFac4
0c0d21fb5c Davi*0539          L5 = (L2*L3)
b8d51be678 Bayl*0540
5c4dff8b20 Mart*0541 #ifndef AUTODIFF_DISABLE_REYNOLDS_SCALE
8d4c48ddd0 Jean*0542 C Velocity Reynolds Scale (Pb here at CS-grid corners !)
                0543          IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
                0544            keZpt=0.25 _d 0*( (KE(i,j)+KE(i-1,j-1))
                0545      &                      +(KE(i-1,j)+KE(i,j-1)) )
                0546            IF ( keZpt.GT.0. ) THEN
002795c988 Jean*0547              Uscl = SQRT(keZpt*L2)*viscAhRe_max
                0548              U4scl= SQRT(keZpt)*L3*viscA4Re_max
8d4c48ddd0 Jean*0549            ELSE
                0550              Uscl =0.
                0551              U4scl=0.
                0552            ENDIF
                0553          ELSE
                0554            Uscl =0.
                0555            U4scl=0.
                0556          ENDIF
5c4dff8b20 Mart*0557 #endif /* ndef AUTODIFF_DISABLE_REYNOLDS_SCALE */
f5bcbacdce Bayl*0558
cc254a4876 Patr*0559 #ifndef AUTODIFF_DISABLE_LEITH
f5bcbacdce Bayl*0560 C This is the vector magnitude of the vorticity gradient squared
002795c988 Jean*0561          IF (useFullLeith.AND.calcLeith) THEN
8def1fc7b9 Jean*0562           grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
                0563      &                        + vrtDx(i,j)*vrtDx(i,j) )
                0564      &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
                0565      &                        + vrtDy(i,j)*vrtDy(i,j) )  )
f5bcbacdce Bayl*0566
                0567 C This is the vector magnitude of grad(div.v) squared
e47a30e113 Jean*0568           grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
                0569      &                        + divDx(i,j)*divDx(i,j) )
                0570      &                     + (divDy(i-1,j)*divDy(i-1,j)
                0571      &                        + divDy(i,j)*divDy(i,j) )  )
b8d51be678 Bayl*0572
70ceddf25f Timo*0573           sqargAh  = leith2fac*grdVrt+leithD2fac*grdDiv
                0574           sqargA4  = leith4fac*grdVrt+leithD4fac*grdDiv
                0575           sqargAhD = leithD2fac*grdDiv
                0576           sqargA4D = leithD4fac*grdDiv
06244a5e4f Jean*0577 #ifdef ALLOW_LEITH_QG
70ceddf25f Timo*0578           sqargQG  = leithQG2fac*(grdVrt+grdDiv)
                0579 #endif
                0580 #ifdef ALLOW_AUTODIFF
                0581 C Avoid derivative of SQRT(0)
                0582           IF (sqargAh .GT.0. _d 0) sqargAh = SQRT(sqargAh)
                0583           IF (sqargA4 .GT.0. _d 0) sqargA4 = SQRT(sqargA4)
                0584           IF (sqargAhD .GT.0. _d 0) sqargAhD = SQRT(sqargAhD)
                0585           IF (sqargA4D .GT.0. _d 0) sqargA4D = SQRT(sqargA4D)
                0586 # ifdef ALLOW_LEITH_QG
                0587           IF (sqargQG .GT.0. _d 0) sqargQG = SQRT(sqargQG)
                0588 # endif
                0589 #else /* ALLOW_AUTODIFF */
                0590           sqargAh = SQRT(sqargAh)
                0591           sqargA4 = SQRT(sqargA4)
                0592           sqargAhD = SQRT(sqargAhD)
                0593           sqargA4D = SQRT(sqargA4D)
                0594 # ifdef ALLOW_LEITH_QG
                0595           sqargQG = SQRT(sqargQG)
                0596 # endif
                0597 #endif /* ALLOW_AUTODIFF */
                0598           viscAh_ZLth(i,j) = sqargAh  * L3
                0599           viscA4_ZLth(i,j) = sqargA4  * L5
                0600           viscAh_ZLthd(i,j)= sqargAhD * L3
                0601           viscA4_ZLthd(i,j)= sqargA4D * L5
                0602 #ifdef ALLOW_LEITH_QG
                0603           viscAh_ZLthQG(i,j)=sqargQG  * L3
06244a5e4f Jean*0604 #endif
b8d51be678 Bayl*0605
002795c988 Jean*0606          ELSEIF (calcLeith) THEN
f5bcbacdce Bayl*0607 C but this approximation will work on cube (and differs by 4X)
002795c988 Jean*0608           grdVrt=MAX( ABS(vrtDx(i-1,j)), ABS(vrtDx(i,j)) )
                0609           grdVrt=MAX( grdVrt, ABS(vrtDy(i,j-1)) )
                0610           grdVrt=MAX( grdVrt, ABS(vrtDy(i,j))   )
b8d51be678 Bayl*0611
002795c988 Jean*0612           grdDiv=MAX( ABS(divDx(i,j)), ABS(divDx(i,j-1)) )
                0613           grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
                0614           grdDiv=MAX( grdDiv, ABS(divDy(i-1,j)) )
b8d51be678 Bayl*0615
20e7e2b61c Bayl*0616           viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
                0617           viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
                0618           viscAh_ZLthD(i,j)=(leithD2fac*grdDiv)*L3
                0619           viscA4_ZLthD(i,j)=(leithD4fac*grdDiv)*L5
06244a5e4f Jean*0620 #ifdef ALLOW_LEITH_QG
f59d76b0dd Ed D*0621           viscAh_ZLthQG(i,j)=leithQG2fac*(grdVrt + grdDiv)*L3
06244a5e4f Jean*0622 #endif
f5bcbacdce Bayl*0623          ELSE
bb1bd5219d Jean*0624           viscAh_ZLth(i,j)=0. _d 0
                0625           viscA4_ZLth(i,j)=0. _d 0
                0626           viscAh_ZLthD(i,j)=0. _d 0
                0627           viscA4_ZLthD(i,j)=0. _d 0
06244a5e4f Jean*0628 #ifdef ALLOW_LEITH_QG
f59d76b0dd Ed D*0629           viscAh_ZLthQG(i,j)=0. _d 0
06244a5e4f Jean*0630 #endif
f5bcbacdce Bayl*0631          ENDIF
                0632
002795c988 Jean*0633          IF (calcSmag) THEN
70ceddf25f Timo*0634           sqargSmag = strain(i,j)**2
bb1bd5219d Jean*0635      &        +0.25 _d 0*(tension( i , j )**2+tension( i ,j-1)**2
70ceddf25f Timo*0636      &                   +tension(i-1, j )**2+tension(i-1,j-1)**2)
                0637 #ifdef ALLOW_AUTODIFF
                0638 C Avoid derivative of SQRT(0)
                0639           IF (sqargSmag.GT.0. _d 0) sqargSmag = SQRT(sqargSmag)
                0640 #else
                0641           sqargSmag = SQRT(sqargSmag)
                0642 #endif
                0643           viscAh_ZSmg(i,j)=L2*sqargSmag
b8d51be678 Bayl*0644           viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
                0645           viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
f5bcbacdce Bayl*0646          ENDIF
cc254a4876 Patr*0647 #endif /* AUTODIFF_DISABLE_LEITH */
f5bcbacdce Bayl*0648
                0649 C  Harmonic on Zeta points
b8d51be678 Bayl*0650          Alin=viscAhZ+viscAhGrid*L2rdt
06244a5e4f Jean*0651      &         + viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
                0652 #ifdef ALLOW_LEITH_QG
                0653      &         + viscAh_ZLthQG(i,j)
                0654 #endif
55479d97c4 Jean*0655 #ifdef ALLOW_3D_VISCAH
06244a5e4f Jean*0656      &         + viscAhZfld(i,j,k,bi,bj)
4240547d2d Mart*0657 # ifdef AUTODIFF_ALLOW_VISCFACADJ
                0658      &          *viscFacAdj
                0659 # endif /* AUTODIFF_ALLOW_VISCFACADJ */
55479d97c4 Jean*0660 #endif
002795c988 Jean*0661          viscAh_ZMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
                0662          viscAh_Z(i,j)=MAX(viscAh_ZMin(i,j),Alin)
                0663          viscAh_ZMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
                0664          viscAh_Z(i,j)=MIN(viscAh_ZMax(i,j),viscAh_Z(i,j))
b8d51be678 Bayl*0665
                0666 C  BiHarmonic on Zeta points
                0667          Alin=viscA4Z+viscA4Grid*L4rdt
06244a5e4f Jean*0668      &         + viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
d062b1342f Gael*0669 #ifdef ALLOW_3D_VISCA4
06244a5e4f Jean*0670      &         + viscA4Zfld(i,j,k,bi,bj)
4240547d2d Mart*0671 # ifdef AUTODIFF_ALLOW_VISCFACADJ
                0672      &          *viscFacAdj
                0673 # endif /* AUTODIFF_ALLOW_VISCFACADJ */
55479d97c4 Jean*0674 #endif
002795c988 Jean*0675          viscA4_ZMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
                0676          viscA4_Z(i,j)=MAX(viscA4_ZMin(i,j),Alin)
                0677          viscA4_ZMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
                0678          viscA4_Z(i,j)=MIN(viscA4_ZMax(i,j),viscA4_Z(i,j))
f5bcbacdce Bayl*0679         ENDDO
                0680        ENDDO
002795c988 Jean*0681
f0501c53d1 Jean*0682 #ifdef ALLOW_NONHYDROSTATIC
                0683        IF ( nonHydrostatic ) THEN
ba2530fec3 Jean*0684 C--   Pass Viscosities to calc_gw (if constant, not necessary)
f0501c53d1 Jean*0685
                0686         IF ( k.LT.Nr ) THEN
                0687 C     Prepare for next level (next call)
                0688          DO j=1-OLy,sNy+OLy
                0689           DO i=1-OLx,sNx+OLx
                0690             viscAh_W(i,j,k+1,bi,bj) = halfRL*viscAh_D(i,j)
                0691             viscA4_W(i,j,k+1,bi,bj) = halfRL*viscA4_D(i,j)
                0692           ENDDO
                0693          ENDDO
                0694         ENDIF
                0695
                0696         shiftAh = viscAhW - viscAhD
                0697         shiftA4 = viscA4W - viscA4D
                0698         IF ( k.EQ.1 ) THEN
                0699 C     These values dont get used
                0700          DO j=1-OLy,sNy+OLy
                0701           DO i=1-OLx,sNx+OLx
                0702             viscAh_W(i,j,k,bi,bj) = shiftAh + viscAh_D(i,j)
                0703             viscA4_W(i,j,k,bi,bj) = shiftA4 + viscA4_D(i,j)
                0704           ENDDO
                0705          ENDDO
                0706         ELSE
                0707 C     Note that previous call of this function has already added half.
                0708          DO j=1-OLy,sNy+OLy
                0709           DO i=1-OLx,sNx+OLx
                0710             viscAh_W(i,j,k,bi,bj) = shiftAh + viscAh_W(i,j,k,bi,bj)
                0711      &                                      + halfRL*viscAh_D(i,j)
                0712             viscA4_W(i,j,k,bi,bj) = shiftA4 + viscA4_W(i,j,k,bi,bj)
                0713      &                                      + halfRL*viscA4_D(i,j)
                0714           ENDDO
                0715          ENDDO
                0716         ENDIF
                0717
                0718        ENDIF
                0719 #endif /* ALLOW_NONHYDROSTATIC */
                0720
ba2530fec3 Jean*0721 c     ELSE
f0501c53d1 Jean*0722 C---- use constant viscosity (useVariableVisc=F):
ba2530fec3 Jean*0723 c      DO j=1-OLy,sNy+OLy
                0724 c       DO i=1-OLx,sNx+OLx
                0725 c        viscAh_D(i,j) = viscAhD
                0726 c        viscAh_Z(i,j) = viscAhZ
                0727 c        viscA4_D(i,j) = viscA4D
                0728 c        viscA4_Z(i,j) = viscA4Z
                0729 c       ENDDO
                0730 c      ENDDO
002795c988 Jean*0731 C---- variable/constant viscosity : end if/else block
ba2530fec3 Jean*0732 c     ENDIF
f5bcbacdce Bayl*0733
                0734 #ifdef ALLOW_DIAGNOSTICS
                0735       IF (useDiagnostics) THEN
                0736        CALL DIAGNOSTICS_FILL(viscAh_D,'VISCAHD ',k,1,2,bi,bj,myThid)
                0737        CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
                0738        CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
                0739        CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
b8d51be678 Bayl*0740
                0741        CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
                0742        CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)
                0743        CALL DIAGNOSTICS_FILL(viscAh_ZMax,'VAHZMAX ',k,1,2,bi,bj,myThid)
                0744        CALL DIAGNOSTICS_FILL(viscA4_ZMax,'VA4ZMAX ',k,1,2,bi,bj,myThid)
                0745
                0746        CALL DIAGNOSTICS_FILL(viscAh_DMin,'VAHDMIN ',k,1,2,bi,bj,myThid)
                0747        CALL DIAGNOSTICS_FILL(viscA4_DMin,'VA4DMIN ',k,1,2,bi,bj,myThid)
                0748        CALL DIAGNOSTICS_FILL(viscAh_ZMin,'VAHZMIN ',k,1,2,bi,bj,myThid)
                0749        CALL DIAGNOSTICS_FILL(viscA4_ZMin,'VA4ZMIN ',k,1,2,bi,bj,myThid)
                0750
                0751        CALL DIAGNOSTICS_FILL(viscAh_DLth,'VAHDLTH ',k,1,2,bi,bj,myThid)
                0752        CALL DIAGNOSTICS_FILL(viscA4_DLth,'VA4DLTH ',k,1,2,bi,bj,myThid)
                0753        CALL DIAGNOSTICS_FILL(viscAh_ZLth,'VAHZLTH ',k,1,2,bi,bj,myThid)
                0754        CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)
                0755
ba2530fec3 Jean*0756        CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD',
                0757      &                       k,1,2,bi,bj,myThid)
                0758        CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD',
                0759      &                       k,1,2,bi,bj,myThid)
                0760        CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD',
                0761      &                       k,1,2,bi,bj,myThid)
                0762        CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD',
                0763      &                       k,1,2,bi,bj,myThid)
06244a5e4f Jean*0764 #ifdef ALLOW_LEITH_QG
f59d76b0dd Ed D*0765        CALL DIAGNOSTICS_FILL(viscAh_DLthQG,'VAHDLTHQ',
                0766      &                       k,1,2,bi,bj,myThid)
                0767        CALL DIAGNOSTICS_FILL(viscAh_ZLthQG,'VAHZLTHQ',
                0768      &                       k,1,2,bi,bj,myThid)
06244a5e4f Jean*0769 #endif
b8d51be678 Bayl*0770        CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)
                0771        CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)
                0772        CALL DIAGNOSTICS_FILL(viscAh_ZSmg,'VAHZSMAG',k,1,2,bi,bj,myThid)
                0773        CALL DIAGNOSTICS_FILL(viscA4_ZSmg,'VA4ZSMAG',k,1,2,bi,bj,myThid)
f5bcbacdce Bayl*0774       ENDIF
06244a5e4f Jean*0775 #endif /* ALLOW_DIAGNOSTICS */
f5bcbacdce Bayl*0776
                0777       RETURN
                0778       END