pkg/gmredi/gmredi_calc_tensor.F

0c49347dc7 Alis*0001 #include "GMREDI_OPTIONS.h"
14e0496834 Jean*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
ee8a6f4ffb Jean*0005 #ifdef ALLOW_CTRL
                0006 # include "CTRL_OPTIONS.h"
                0007 #endif
e2259a1942 Jean*0008 #ifdef ALLOW_KPP
                0009 # include "KPP_OPTIONS.h"
                0010 #endif
472f570611 Jean*0011 #if ( defined GM_READ_K3D_REDI || defined GM_VISBECK_VARIABLE_K \
                0012       || defined ALLOW_GM_LEITH_QG || defined GM_BATES_K3D )
                0013 # define GM_USE_K3D_REDI
                0014 #endif
                0015 C no GM_GEOM_VARIABLE here since it is at a different location
                0016 #if ( defined GM_READ_K3D_GM   || defined GM_VISBECK_VARIABLE_K \
                0017       || defined ALLOW_GM_LEITH_QG || defined GM_BATES_K3D )
                0018 # define GM_USE_K3D_GM
                0019 #endif
0c49347dc7 Alis*0020
e2259a1942 Jean*0021 CBOP
                0022 C     !ROUTINE: GMREDI_CALC_TENSOR
                0023 C     !INTERFACE:
0c49347dc7 Alis*0024       SUBROUTINE GMREDI_CALC_TENSOR(
e2259a1942 Jean*0025      I             iMin, iMax, jMin, jMax,
0c49347dc7 Alis*0026      I             sigmaX, sigmaY, sigmaR,
e2259a1942 Jean*0027      I             bi, bj, myTime, myIter, myThid )
                0028
                0029 C     !DESCRIPTION: \bv
                0030 C     *==========================================================*
                0031 C     | SUBROUTINE GMREDI_CALC_TENSOR
                0032 C     | o Calculate tensor elements for GM/Redi tensor.
                0033 C     *==========================================================*
                0034 C     *==========================================================*
                0035 C     \ev
                0036
                0037 C     !USES:
0c49347dc7 Alis*0038       IMPLICIT NONE
                0039
                0040 C     == Global variables ==
                0041 #include "SIZE.h"
                0042 #include "GRID.h"
                0043 #include "DYNVARS.h"
                0044 #include "EEPARAMS.h"
                0045 #include "PARAMS.h"
                0046 #include "GMREDI.h"
ee8a6f4ffb Jean*0047 #ifdef ALLOW_CTRL
                0048 # include "CTRL_FIELDS.h"
                0049 #endif
e2259a1942 Jean*0050 #ifdef ALLOW_KPP
                0051 # include "KPP.h"
                0052 #endif
0c49347dc7 Alis*0053
b6b11b9b2f Patr*0054 #ifdef ALLOW_AUTODIFF_TAMC
a4576c7cde Juli*0055 # include "tamc.h"
b6b11b9b2f Patr*0056 #endif /* ALLOW_AUTODIFF_TAMC */
                0057
e2259a1942 Jean*0058 C     !INPUT/OUTPUT PARAMETERS:
                0059 C     bi, bj    :: tile indices
                0060 C     myTime    :: Current time in simulation
                0061 C     myIter    :: Current iteration number in simulation
                0062 C     myThid    :: My Thread Id. number
0c49347dc7 Alis*0063 C
e2259a1942 Jean*0064       INTEGER iMin,iMax,jMin,jMax
ee8a6f4ffb Jean*0065       _RL sigmaX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0066       _RL sigmaY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0067       _RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
e2259a1942 Jean*0068       INTEGER bi, bj
                0069       _RL     myTime
                0070       INTEGER myIter
0c49347dc7 Alis*0071       INTEGER myThid
e2259a1942 Jean*0072 CEOP
0c49347dc7 Alis*0073
                0074 #ifdef ALLOW_GMREDI
e2259a1942 Jean*0075 C     !LOCAL VARIABLES:
8233d0ceb9 Jean*0076 C     maskFk    :: 2-D mask for vertical interface k (between level k-1 & k)
2ae58a73ff Jean*0077       INTEGER i,j,k
8233d0ceb9 Jean*0078       _RS maskFk(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
ee8a6f4ffb Jean*0079       _RL SlopeX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0080       _RL SlopeY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0081       _RL dSigmaDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0082       _RL dSigmaDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0083       _RL dSigmaDr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0084       _RL SlopeSqr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085       _RL taperFct(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0086       _RL ldd97_LrhoC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0087       _RL ldd97_LrhoW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0088       _RL ldd97_LrhoS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
5b172de0d2 Jean*0089       _RL Cspd, LrhoInf, LrhoSup, fCoriLoc, rDepth
0c12568530 mjlo*0090       _RL Kgm_tmp, Kredi_tmp, isopycK, bolus_K
0c49347dc7 Alis*0091
ee8a6f4ffb Jean*0092       INTEGER kLow_W (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0093       INTEGER kLow_S (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0094       _RL locMixLayer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0095       _RL baseSlope  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0096       _RL hTransLay  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0097       _RL recipLambda(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a4576c7cde Juli*0098       INTEGER km1
                0099 #if ( defined GM_NON_UNITY_DIAGONAL || defined GM_EXTRA_DIAGONAL )
2ae58a73ff Jean*0100       INTEGER kp1
                0101       _RL maskp1
                0102 #endif
472f570611 Jean*0103 #ifdef GM_USE_K3D_REDI
                0104       _RL locK3dRedi(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0105 #endif
                0106 #ifdef GM_USE_K3D_GM
                0107       _RL locK3dGM(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0108 #else
                0109       _RL locK3dGM(1)
                0110 #endif
e2259a1942 Jean*0111
0c49347dc7 Alis*0112 #ifdef GM_VISBECK_VARIABLE_K
8233d0ceb9 Jean*0113       INTEGER ks
                0114       _RL dSigmaH, dSigmaR, Sloc, rTop
                0115 c     _RL M2loc
9bee368eff Jean*0116       _RL recipMaxSlope
ea71059281 Jean*0117       _RL deltaH, integrDepth
                0118       _RL N2loc, SNloc
9bee368eff Jean*0119 #endif /* GM_VISBECK_VARIABLE_K */
0c49347dc7 Alis*0120
066e0d5e64 Jean*0121 #ifdef ALLOW_DIAGNOSTICS
                0122       LOGICAL  doDiagRediFlx
                0123       LOGICAL  DIAGNOSTICS_IS_ON
                0124       EXTERNAL DIAGNOSTICS_IS_ON
a4576c7cde Juli*0125 # if ( defined GM_NON_UNITY_DIAGONAL || defined GM_EXTRA_DIAGONAL )
066e0d5e64 Jean*0126       _RL dTdz
                0127       _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a4576c7cde Juli*0128 # endif
2ae58a73ff Jean*0129 #endif /* ALLOW_DIAGNOSTICS */
7c50f07931 Mart*0130 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0131 C     tkey :: tape key (depends on tiles)
                0132 C     kkey :: tape key (depends on levels and tiles)
                0133       INTEGER tkey, kkey
7c50f07931 Mart*0134 #endif
066e0d5e64 Jean*0135
549d1a8d8c Jean*0136 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0137
b6b11b9b2f Patr*0138 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0139       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
b6b11b9b2f Patr*0140 #endif /* ALLOW_AUTODIFF_TAMC */
                0141
066e0d5e64 Jean*0142 #ifdef ALLOW_DIAGNOSTICS
                0143       doDiagRediFlx = .FALSE.
                0144       IF ( useDiagnostics ) THEN
                0145         doDiagRediFlx = DIAGNOSTICS_IS_ON('GM_KuzTz', myThid )
b0e49a1609 Jean*0146         doDiagRediFlx = doDiagRediFlx .OR.
066e0d5e64 Jean*0147      &                  DIAGNOSTICS_IS_ON('GM_KvzTz', myThid )
                0148       ENDIF
                0149 #endif
b0e49a1609 Jean*0150
2092dbb101 Patr*0151 #ifdef GM_VISBECK_VARIABLE_K
9bee368eff Jean*0152       recipMaxSlope = 0. _d 0
                0153       IF ( GM_Visbeck_maxSlope.GT.0. _d 0 ) THEN
                0154         recipMaxSlope = 1. _d 0 / GM_Visbeck_maxSlope
                0155       ENDIF
ee8a6f4ffb Jean*0156       DO j=1-OLy,sNy+OLy
                0157        DO i=1-OLx,sNx+OLx
2092dbb101 Patr*0158         VisbeckK(i,j,bi,bj) = 0. _d 0
                0159        ENDDO
                0160       ENDDO
                0161 #endif
                0162
549d1a8d8c Jean*0163 C--   set ldd97_Lrho (for tapering scheme ldd97):
e2259a1942 Jean*0164       IF ( GM_taper_scheme.EQ.'ldd97' .OR.
                0165      &     GM_taper_scheme.EQ.'fm07' ) THEN
549d1a8d8c Jean*0166        Cspd = 2. _d 0
                0167        LrhoInf = 15. _d 3
                0168        LrhoSup = 100. _d 3
                0169 C-     Tracer point location (center):
ee8a6f4ffb Jean*0170        DO j=1-OLy,sNy+OLy
                0171         DO i=1-OLx,sNx+OLx
549d1a8d8c Jean*0172          IF (fCori(i,j,bi,bj).NE.0.) THEN
                0173            ldd97_LrhoC(i,j) = Cspd/ABS(fCori(i,j,bi,bj))
                0174          ELSE
                0175            ldd97_LrhoC(i,j) = LrhoSup
                0176          ENDIF
                0177          ldd97_LrhoC(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoC(i,j),LrhoSup))
                0178         ENDDO
                0179        ENDDO
                0180 C-     U point location (West):
ee8a6f4ffb Jean*0181        DO j=1-OLy,sNy+OLy
                0182         kLow_W(1-OLx,j) = 0
                0183         ldd97_LrhoW(1-OLx,j) = LrhoSup
                0184         DO i=1-OLx+1,sNx+OLx
e2259a1942 Jean*0185          kLow_W(i,j) = MIN(kLowC(i-1,j,bi,bj),kLowC(i,j,bi,bj))
549d1a8d8c Jean*0186          fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
5b172de0d2 Jean*0187          IF ( fCoriLoc.NE.zeroRL ) THEN
549d1a8d8c Jean*0188            ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
                0189          ELSE
                0190            ldd97_LrhoW(i,j) = LrhoSup
                0191          ENDIF
                0192          ldd97_LrhoW(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoW(i,j),LrhoSup))
                0193         ENDDO
                0194        ENDDO
                0195 C-     V point location (South):
ee8a6f4ffb Jean*0196        DO i=1-OLx+1,sNx+OLx
                0197          kLow_S(i,1-OLy) = 0
                0198          ldd97_LrhoS(i,1-OLy) = LrhoSup
549d1a8d8c Jean*0199        ENDDO
ee8a6f4ffb Jean*0200        DO j=1-OLy+1,sNy+OLy
                0201         DO i=1-OLx,sNx+OLx
e2259a1942 Jean*0202          kLow_S(i,j) = MIN(kLowC(i,j-1,bi,bj),kLowC(i,j,bi,bj))
549d1a8d8c Jean*0203          fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
5b172de0d2 Jean*0204          IF ( fCoriLoc.NE.zeroRL ) THEN
549d1a8d8c Jean*0205            ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
                0206          ELSE
                0207            ldd97_LrhoS(i,j) = LrhoSup
                0208          ENDIF
                0209          ldd97_LrhoS(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoS(i,j),LrhoSup))
                0210         ENDDO
                0211        ENDDO
                0212       ELSE
                0213 C-    Just initialize to zero (not use anyway)
ee8a6f4ffb Jean*0214        DO j=1-OLy,sNy+OLy
                0215         DO i=1-OLx,sNx+OLx
549d1a8d8c Jean*0216           ldd97_LrhoC(i,j) = 0. _d 0
                0217           ldd97_LrhoW(i,j) = 0. _d 0
                0218           ldd97_LrhoS(i,j) = 0. _d 0
                0219         ENDDO
                0220        ENDDO
                0221       ENDIF
                0222
050b4366e6 Jean*0223 #ifdef GM_BOLUS_ADVEC
                0224       DO k=1,Nr
ee8a6f4ffb Jean*0225        DO j=1-OLy,sNy+OLy
                0226         DO i=1-OLx,sNx+OLx
050b4366e6 Jean*0227          GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
                0228          GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
                0229         ENDDO
                0230        ENDDO
                0231       ENDDO
                0232 #endif /* GM_BOLUS_ADVEC */
0c12568530 mjlo*0233 #ifdef GM_USE_K3D_REDI
                0234       DO k=1,Nr
                0235        DO j=1-OLy,sNy+OLy
                0236         DO i=1-OLx,sNx+OLx
9c1d5d8ec0 Jean*0237          locK3dRedi(i,j,k) =
0c12568530 mjlo*0238 # ifdef GM_READ_K3D_REDI
9c1d5d8ec0 Jean*0239      &       GM_inpK3dRedi(i,j,k,bi,bj)
0c12568530 mjlo*0240 # else
9c1d5d8ec0 Jean*0241      &       GM_isopycK*GM_isoFac1d(k)*GM_isoFac2d(i,j,bi,bj)
0c12568530 mjlo*0242 # endif
                0243         ENDDO
                0244        ENDDO
                0245       ENDDO
                0246 #endif /* GM_USE_K3D_REDI */
                0247 #ifdef GM_USE_K3D_GM
                0248       DO k=1,Nr
                0249        DO j=1-OLy,sNy+OLy
                0250         DO i=1-OLx,sNx+OLx
9c1d5d8ec0 Jean*0251          locK3dGM(i,j,k) =
0c12568530 mjlo*0252 # ifdef GM_READ_K3D_GM
9c1d5d8ec0 Jean*0253      &       GM_inpK3dGM(i,j,k,bi,bj)
0c12568530 mjlo*0254 # else
9c1d5d8ec0 Jean*0255      &       GM_background_K*GM_bolFac1d(k)*GM_bolFac2d(i,j,bi,bj)
0c12568530 mjlo*0256 # endif
                0257         ENDDO
                0258        ENDDO
                0259       ENDDO
                0260 #else
                0261       locK3dGM(1) = 0. _d 0
                0262 #endif /* GM_USE_K3D_GM */
14e0496834 Jean*0263 #ifdef ALLOW_AUTODIFF
050b4366e6 Jean*0264       DO k=1,Nr
ee8a6f4ffb Jean*0265        DO j=1-OLy,sNy+OLy
                0266         DO i=1-OLx,sNx+OLx
050b4366e6 Jean*0267          Kwx(i,j,k,bi,bj)  = 0. _d 0
                0268          Kwy(i,j,k,bi,bj)  = 0. _d 0
                0269          Kwz(i,j,k,bi,bj)  = 0. _d 0
                0270          Kux(i,j,k,bi,bj)  = 0. _d 0
                0271          Kvy(i,j,k,bi,bj)  = 0. _d 0
                0272 # ifdef GM_EXTRA_DIAGONAL
                0273          Kuz(i,j,k,bi,bj)  = 0. _d 0
                0274          Kvz(i,j,k,bi,bj)  = 0. _d 0
                0275 # endif
                0276         ENDDO
                0277        ENDDO
                0278       ENDDO
14e0496834 Jean*0279 #endif /* ALLOW_AUTODIFF */
0c49347dc7 Alis*0280
050b4366e6 Jean*0281 C--   Initialise Mixed Layer related array:
ee8a6f4ffb Jean*0282       DO j=1-OLy,sNy+OLy
                0283        DO i=1-OLx,sNx+OLx
5755dbe390 Patr*0284          hTransLay(i,j) = R_low(i,j,bi,bj)
                0285          baseSlope(i,j) =  0. _d 0
                0286          recipLambda(i,j) = 0. _d 0
                0287          locMixLayer(i,j) = 0. _d 0
                0288        ENDDO
                0289       ENDDO
e2259a1942 Jean*0290 #ifdef ALLOW_KPP
                0291       IF ( useKPP ) THEN
ee8a6f4ffb Jean*0292        DO j=1-OLy,sNy+OLy
                0293         DO i=1-OLx,sNx+OLx
e2259a1942 Jean*0294          locMixLayer(i,j) = KPPhbl(i,j,bi,bj)
                0295         ENDDO
                0296        ENDDO
                0297       ELSE
                0298 #else
                0299       IF ( .TRUE. ) THEN
                0300 #endif
ee8a6f4ffb Jean*0301        DO j=1-OLy,sNy+OLy
                0302         DO i=1-OLx,sNx+OLx
e2259a1942 Jean*0303          locMixLayer(i,j) = hMixLayer(i,j,bi,bj)
                0304         ENDDO
                0305        ENDDO
                0306       ENDIF
                0307
05118ac017 Jean*0308 #ifdef GM_BATES_K3D
                0309       IF (GM_useBatesK3d) THEN
7ea279c259 Jean*0310 C     Calculate the 3D diffusivity as per Bates et al. (JPO, 2014)
                0311         CALL TIMER_START('GMREDI_CALC_BATES_K [GMREDI_CALC_TENSOR]',
0d1e4b948d Mich*0312      &                    myThid)
7ea279c259 Jean*0313         CALL GMREDI_CALC_BATES_K(
0d1e4b948d Mich*0314      I       iMin, iMax, jMin, jMax,
                0315      I       sigmaX, sigmaY, sigmaR,
7ea279c259 Jean*0316      I       bi, bj, myTime, myIter, myThid )
                0317         CALL TIMER_STOP('GMREDI_CALC_BATES_K [GMREDI_CALC_TENSOR]',
0d1e4b948d Mich*0318      &                  myThid)
0c12568530 mjlo*0319 # ifndef GM_BATES_PASSIVE
                0320         DO k=1,Nr
                0321          DO j=1-OLy,sNy+OLy
                0322           DO i=1-OLx,sNx+OLx
                0323            locK3dRedi(i,j,k) = locK3dRedi(i,j,k)
                0324      &          + GM_BatesK3d(i,j,k,bi,bj)*GM_isoFac_calcK
                0325            locK3dGM(i,j,k)   = locK3dGM(i,j,k)
                0326      &          + GM_BatesK3d(i,j,k,bi,bj)
                0327           ENDDO
                0328          ENDDO
                0329         ENDDO
                0330 # endif
0d1e4b948d Mich*0331       ENDIF
                0332 #endif
                0333
f59d76b0dd Ed D*0334 #ifdef ALLOW_GM_LEITH_QG
e25acdb1f2 Jean*0335 # ifdef ALLOW_AUTODIFF
                0336       DO k=1,Nr
8233d0ceb9 Jean*0337        DO j=1-OLy,sNy+OLy
                0338         DO i=1-OLx,sNx+OLx
                0339           GM_LeithQG_K(i,j,k,bi,bj) = 0. _d 0
f59d76b0dd Ed D*0340         ENDDO
8233d0ceb9 Jean*0341        ENDDO
e25acdb1f2 Jean*0342       ENDDO
                0343 # endif
                0344       IF ( GM_useLeithQG ) THEN
                0345 C   Use Leith QG viscosity in GMRedi scheme
8233d0ceb9 Jean*0346         CALL GMREDI_CALC_QGLEITH(
f59d76b0dd Ed D*0347      O                GM_LeithQG_K(1-OLx,1-OLy,1,bi,bj),
                0348      I                bi, bj, myTime, myIter, myThid )
0c12568530 mjlo*0349         DO k=1,Nr
                0350          DO j=1-OLy,sNy+OLy
                0351           DO i=1-OLx,sNx+OLx
                0352            locK3dRedi(i,j,k) = locK3dRedi(i,j,k)
                0353      &          + GM_LeithQG_K(i,j,k,bi,bj)*GM_isoFac_calcK
                0354            locK3dGM(i,j,k)   = locK3dGM(i,j,k)
                0355      &          + GM_LeithQG_K(i,j,k,bi,bj)
                0356           ENDDO
                0357          ENDDO
                0358         ENDDO
f59d76b0dd Ed D*0359       ENDIF
                0360 #endif /* ALLOW_GM_LEITH_QG */
                0361
a4576c7cde Juli*0362 #ifdef GM_GEOM_VARIABLE_K
                0363 C     Calculate the kgm as per the GEOMETRIC formulation
                0364 C        (e.g. Mak et al, 2018; see also Marshall et al, 2012)
                0365       IF ( GM_useGEOM ) THEN
                0366         DO k=1,Nr
                0367           DO j=1-OLy,sNy+OLy
                0368             DO i=1-OLx,sNx+OLx
                0369               GEOM_K3d(i,j,k,bi,bj) = 0. _d 0
                0370             ENDDO
                0371           ENDDO
                0372         ENDDO
                0373         CALL GMREDI_CALC_GEOM(
                0374      I              sigmaX, sigmaY, sigmaR,
                0375      I              bi, bj, myTime, myIter, myThid )
                0376       ENDIF
                0377 #endif /* GM_GEOM_VARIABLE_K */
                0378
050b4366e6 Jean*0379 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0380 C-- 1rst loop on k : compute Tensor Coeff. at W points.
                0381
e2259a1942 Jean*0382       DO k=Nr,2,-1
                0383
ee8a6f4ffb Jean*0384        DO j=1-OLy,sNy+OLy
                0385         DO i=1-OLx,sNx+OLx
8233d0ceb9 Jean*0386 #ifdef ALLOW_AUTODIFF
b6b11b9b2f Patr*0387          SlopeX(i,j)       = 0. _d 0
                0388          SlopeY(i,j)       = 0. _d 0
a4576c7cde Juli*0389          dSigmaDx(i,j)    = 0. _d 0
2092dbb101 Patr*0390          dSigmaDy(i,j)     = 0. _d 0
549d1a8d8c Jean*0391          dSigmaDr(i,j)     = 0. _d 0
b6b11b9b2f Patr*0392          SlopeSqr(i,j)     = 0. _d 0
                0393          taperFct(i,j)     = 0. _d 0
8233d0ceb9 Jean*0394 #endif /* ALLOW_AUTODIFF */
                0395          maskFk(i,j)       = maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj)
b6b11b9b2f Patr*0396         ENDDO
                0397        ENDDO
                0398
ee8a6f4ffb Jean*0399        DO j=1-OLy+1,sNy+OLy-1
                0400         DO i=1-OLx+1,sNx+OLx-1
0c49347dc7 Alis*0401 C      Gradient of Sigma at rVel points
b0e49a1609 Jean*0402          dSigmaDx(i,j)=op25*( sigmaX(i+1,j,k-1)+sigmaX(i,j,k-1)
                0403      &                       +sigmaX(i+1,j, k )+sigmaX(i,j, k )
8233d0ceb9 Jean*0404      &                      )*maskFk(i,j)
b0e49a1609 Jean*0405          dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)
                0406      &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )
8233d0ceb9 Jean*0407      &                      )*maskFk(i,j)
9bee368eff Jean*0408 c        dSigmaDr(i,j)=sigmaR(i,j,k)
b0e49a1609 Jean*0409         ENDDO
0c49347dc7 Alis*0410        ENDDO
                0411
b0e49a1609 Jean*0412 #ifdef GM_VISBECK_VARIABLE_K
d1c10e3d93 mjlo*0413 # ifdef ALLOW_AUTODIFF_TAMC
                0414 CADJ INCOMPLETE dSigmaDr
                0415 # endif
5b172de0d2 Jean*0416        IF ( GM_Visbeck_alpha.GT.zeroRL .AND.
b0e49a1609 Jean*0417      &      -rC(k-1).LT.GM_Visbeck_depth ) THEN
                0418
ee8a6f4ffb Jean*0419         DO j=1-OLy,sNy+OLy
                0420          DO i=1-OLx,sNx+OLx
5b172de0d2 Jean*0421           dSigmaDr(i,j) = MAX( gravitySign*sigmaR(i,j,k), zeroRL )
9bee368eff Jean*0422          ENDDO
                0423         ENDDO
                0424
b0e49a1609 Jean*0425 C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N
                0426 C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.
                0427
                0428 C       Calculate terms for mean Richardson number which is used
                0429 C       in the "variable K" parameterisaton:
                0430 C       compute depth average from surface down to the bottom or
                0431 C       GM_Visbeck_depth, whatever is the shallower.
                0432
8233d0ceb9 Jean*0433         rTop = rF(1)
ee8a6f4ffb Jean*0434         DO j=1-OLy+1,sNy+OLy-1
                0435          DO i=1-OLx+1,sNx+OLx-1
8233d0ceb9 Jean*0436           IF ( maskFk(i,j).NE.zeroRS ) THEN
                0437            ks = kSurfC(i,j,bi,bj)
                0438 C-      If dry-cell @ top (e.g., under ice-shelf), integrate over Visbeck_depth
                0439 C       as relative to top solid BC (as below) or stop integral if deeper (in
                0440 C       absolute sense) than Visbeck_depth (just comment "rTop =" line below):
                0441            rTop = Ro_surf(i,j,bi,bj)
                0442            integrDepth = rTop - rC( kLowC(i,j,bi,bj) )
b0e49a1609 Jean*0443 C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
                0444            integrDepth = MIN( integrDepth, GM_Visbeck_depth )
9bee368eff Jean*0445 C-      to recover "old-visbeck" form with Visbeck_minDepth = Visbeck_depth
                0446            integrDepth = MAX( integrDepth, GM_Visbeck_minDepth )
b0e49a1609 Jean*0447 C       Distance between level center above and the integration depth
8233d0ceb9 Jean*0448            deltaH = integrDepth - rTop + rC(k-1)
b0e49a1609 Jean*0449 C       If negative then we are below the integration level
                0450 C       (cannot be the case with 2 conditions on maskC & -rC(k-1))
                0451 C       If positive we limit this to the distance from center above
                0452            deltaH = MIN( deltaH, drC(k) )
                0453 C       Now we convert deltaH to a non-dimensional fraction
8233d0ceb9 Jean*0454            deltaH = deltaH/( integrDepth - rTop + rC(ks) )
b0e49a1609 Jean*0455
ea71059281 Jean*0456 C--      compute: ( M^2 * S )^1/2   (= S*N since S=M^2/N^2 )
9bee368eff Jean*0457 C        a 5 points average gives a more "homogeneous" formulation
                0458 C        (same stencil and same weights as for dSigmaH calculation)
                0459            dSigmaR = ( dSigmaDr(i,j)*4. _d 0
                0460      &               + dSigmaDr(i-1,j)
                0461      &               + dSigmaDr(i+1,j)
                0462      &               + dSigmaDr(i,j-1)
                0463      &               + dSigmaDr(i,j+1)
                0464      &               )/( 4. _d 0
8233d0ceb9 Jean*0465      &                 + maskFk(i-1,j)
                0466      &                 + maskFk(i+1,j)
                0467      &                 + maskFk(i,j-1)
                0468      &                 + maskFk(i,j+1)
9bee368eff Jean*0469      &                 )
b0e49a1609 Jean*0470            dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
                0471      &             + dSigmaDy(i,j)*dSigmaDy(i,j)
                0472            IF ( dSigmaH .GT. 0. _d 0 ) THEN
                0473              dSigmaH = SQRT( dSigmaH )
d1c10e3d93 mjlo*0474 # ifdef ALLOW_AUTODIFF_TAMC
                0475 C     This gets rid of a TAF-recomputation warning for dSigmaH that is
                0476 C     not associated with any actual recomputations. TAF is probably
                0477 C     just confused by the cascade of if-statements.
                0478            ENDIF
                0479            IF ( dSigmaH .GT. 0. _d 0 ) THEN
                0480 # endif
9bee368eff Jean*0481 C-       compute slope, limited by GM_Visbeck_maxSlope:
5b172de0d2 Jean*0482              IF ( dSigmaR.GT.dSigmaH*recipMaxSlope ) THEN
                0483               Sloc = dSigmaH / dSigmaR
b0e49a1609 Jean*0484              ELSE
9bee368eff Jean*0485               Sloc = GM_Visbeck_maxSlope
b0e49a1609 Jean*0486              ENDIF
8233d0ceb9 Jean*0487 c            M2loc = gravity*recip_rhoConst*dSigmaH
5b172de0d2 Jean*0488              N2loc = gravity*recip_rhoConst*dSigmaR
ea71059281 Jean*0489              IF ( N2loc.GT.0. _d 0 ) THEN
                0490                SNloc = Sloc*SQRT(N2loc)
                0491              ELSE
                0492                SNloc = 0. _d 0
                0493              ENDIF
b0e49a1609 Jean*0494            ELSE
                0495              SNloc = 0. _d 0
                0496            ENDIF
                0497            VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
                0498      &       +deltaH*GM_Visbeck_alpha
                0499      &              *GM_Visbeck_length*GM_Visbeck_length*SNloc
                0500           ENDIF
                0501          ENDDO
                0502         ENDDO
                0503        ENDIF
                0504 #endif /* GM_VISBECK_VARIABLE_K */
5b172de0d2 Jean*0505 C     Z-Coords: change sign of vertical Sigma gradient (always downward)
                0506 C     to match stratification sign (i.e., positive if stratified)
ee8a6f4ffb Jean*0507        DO j=1-OLy,sNy+OLy
                0508         DO i=1-OLx,sNx+OLx
5b172de0d2 Jean*0509          dSigmaDr(i,j) = gravitySign*sigmaR(i,j,k)
9bee368eff Jean*0510         ENDDO
                0511        ENDDO
                0512
                0513 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0514        kkey = k + (tkey-1)*Nr
9bee368eff Jean*0515 CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
                0516 CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
                0517 CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
361543e926 Jean*0518 # ifndef GM_EXCLUDE_FM07_TAP
9bee368eff Jean*0519 CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
                0520 CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
                0521 CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
361543e926 Jean*0522 # endif
9bee368eff Jean*0523 #endif /* ALLOW_AUTODIFF_TAMC */
b0e49a1609 Jean*0524
5b172de0d2 Jean*0525 C     set "rDepth" (= depth from the surface, in rUnit) for 'ldd97' tapering
                0526        IF ( usingZcoords ) THEN
                0527         rDepth = rF(1) - rF(k)
                0528        ELSE
                0529         rDepth = rF(k) - rF(Nr+1)
                0530        ENDIF
0c49347dc7 Alis*0531 C     Calculate slopes for use in tensor, taper and/or clip
b0e49a1609 Jean*0532        CALL GMREDI_SLOPE_LIMIT(
549d1a8d8c Jean*0533      O             SlopeX, SlopeY,
e9fd580bc4 Jean*0534      O             SlopeSqr, taperFct,
e2259a1942 Jean*0535      U             hTransLay, baseSlope, recipLambda,
549d1a8d8c Jean*0536      U             dSigmaDr,
                0537      I             dSigmaDx, dSigmaDy,
5b172de0d2 Jean*0538      I             ldd97_LrhoC, locMixLayer, rDepth, rF,
ee8a6f4ffb Jean*0539      I             kLowC(1-OLx,1-OLy,bi,bj),
5b172de0d2 Jean*0540      I             3, k, bi, bj, myTime, myIter, myThid )
0c49347dc7 Alis*0541
8233d0ceb9 Jean*0542 #ifdef GMREDI_MASK_SLOPES
ee8a6f4ffb Jean*0543        DO j=1-OLy+1,sNy+OLy-1
                0544         DO i=1-OLx+1,sNx+OLx-1
b0e49a1609 Jean*0545 C      Mask Iso-neutral slopes
8233d0ceb9 Jean*0546          SlopeX(i,j)   = SlopeX(i,j)*maskFk(i,j)
                0547          SlopeY(i,j)   = SlopeY(i,j)*maskFk(i,j)
                0548          SlopeSqr(i,j) = SlopeSqr(i,j)*maskFk(i,j)
b0e49a1609 Jean*0549         ENDDO
b6b11b9b2f Patr*0550        ENDDO
8233d0ceb9 Jean*0551 #endif
b6b11b9b2f Patr*0552
                0553 #ifdef ALLOW_AUTODIFF_TAMC
9cb619cfcd Patr*0554 CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
                0555 CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
20b8842d78 Patr*0556 CADJ STORE SlopeSqr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
549d1a8d8c Jean*0557 CADJ STORE dSigmaDr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
                0558 CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
b6b11b9b2f Patr*0559 #endif /* ALLOW_AUTODIFF_TAMC */
                0560
5b172de0d2 Jean*0561 C     Components of Redi/GM tensor, first step with just the slope:
                0562 C     P-Coords: since slope comp. have same sign as in Z-coord, need to
                0563 C     flip sign of extra-diagonal component (see manual sec. 8.4.1.2.4)
ee8a6f4ffb Jean*0564        DO j=1-OLy+1,sNy+OLy-1
                0565         DO i=1-OLx+1,sNx+OLx-1
5b172de0d2 Jean*0566           Kwx(i,j,k,bi,bj) = -gravitySign*SlopeX(i,j)*taperFct(i,j)
                0567           Kwy(i,j,k,bi,bj) = -gravitySign*SlopeY(i,j)*taperFct(i,j)
                0568           Kwz(i,j,k,bi,bj) = SlopeSqr(i,j)*taperFct(i,j)
e2259a1942 Jean*0569         ENDDO
                0570        ENDDO
0c49347dc7 Alis*0571
f42e64b3e7 Jean*0572 C-- end 1rst loop on vertical level index k
                0573       ENDDO
                0574
0c49347dc7 Alis*0575 #ifdef GM_VISBECK_VARIABLE_K
a4576c7cde Juli*0576 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0577 CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=tkey, byte=isbyte
a4576c7cde Juli*0578 # endif
5b172de0d2 Jean*0579       IF ( GM_Visbeck_alpha.GT.zeroRL ) THEN
94a8024bbe Jean*0580 C-     Limit range that Kappa-GM can take
ee8a6f4ffb Jean*0581        DO j=1-OLy+1,sNy+OLy-1
                0582         DO i=1-OLx+1,sNx+OLx-1
f42e64b3e7 Jean*0583          VisbeckK(i,j,bi,bj)=
9bee368eff Jean*0584      &       MIN( MAX( VisbeckK(i,j,bi,bj), GM_Visbeck_minVal_K ),
                0585      &            GM_Visbeck_maxVal_K )
f42e64b3e7 Jean*0586         ENDDO
0c49347dc7 Alis*0587        ENDDO
0c12568530 mjlo*0588        DO k=1,Nr
                0589         DO j=1-OLy,sNy+OLy
                0590          DO i=1-OLx,sNx+OLx
                0591           locK3dRedi(i,j,k) = locK3dRedi(i,j,k)
                0592      &         + VisbeckK(i,j,bi,bj)*GM_isoFac_calcK
                0593           locK3dGM(i,j,k)   = locK3dGM(i,j,k)
                0594      &         + VisbeckK(i,j,bi,bj)
                0595          ENDDO
                0596         ENDDO
                0597        ENDDO
f42e64b3e7 Jean*0598       ENDIF
d1c10e3d93 mjlo*0599 C     If you want to avoid recomputing the above loop, you can uncomment
0c12568530 mjlo*0600 C     the following lines.
d1c10e3d93 mjlo*0601 c# ifdef ALLOW_AUTODIFF_TAMC
                0602 cCADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=tkey, byte=isbyte
                0603 c# endif
f42e64b3e7 Jean*0604 #endif /* GM_VISBECK_VARIABLE_K */
0c49347dc7 Alis*0605
472f570611 Jean*0606 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
050b4366e6 Jean*0607 C--   Calculate Stream-Functions used in Advective Form:
                0608
                0609 #ifdef GM_BOLUS_ADVEC
                0610       IF (GM_AdvForm) THEN
a4576c7cde Juli*0611 # ifdef GM_BOLUS_BVP
050b4366e6 Jean*0612        IF (GM_UseBVP) THEN
                0613         CALL GMREDI_CALC_PSI_BVP(
                0614      I             bi, bj, iMin, iMax, jMin, jMax,
                0615      I             sigmaX, sigmaY, sigmaR,
                0616      I             myThid )
                0617        ELSE
a4576c7cde Juli*0618 # endif
d1eb4480d0 Jean*0619 # ifndef GM_BATES_PASSIVE
05118ac017 Jean*0620          IF ( .NOT.GM_useBatesK3d ) THEN
a4576c7cde Juli*0621 # endif
05118ac017 Jean*0622 C          If using BatesK3d PsiX and PsiY are calculated in GMREDI_CALC_BATES_K
d1eb4480d0 Jean*0623            CALL GMREDI_CALC_PSI_BOLUS(
050b4366e6 Jean*0624      I              bi, bj, iMin, iMax, jMin, jMax,
472f570611 Jean*0625      I              sigmaX, sigmaY, sigmaR, locK3dGM,
050b4366e6 Jean*0626      I              ldd97_LrhoW, ldd97_LrhoS,
                0627      I              myThid )
a4576c7cde Juli*0628 # ifndef GM_BATES_PASSIVE
0d1e4b948d Mich*0629          ENDIF
a4576c7cde Juli*0630 # endif
                0631 # ifdef GM_BOLUS_BVP
050b4366e6 Jean*0632        ENDIF
a4576c7cde Juli*0633 # endif
050b4366e6 Jean*0634       ENDIF
361543e926 Jean*0635 #endif /* GM_BOLUS_ADVEC */
050b4366e6 Jean*0636
                0637 #ifndef GM_EXCLUDE_SUBMESO
                0638       IF ( GM_useSubMeso .AND. GM_AdvForm ) THEN
                0639         CALL SUBMESO_CALC_PSI(
                0640      I              bi, bj, iMin, iMax, jMin, jMax,
                0641      I              sigmaX, sigmaY, sigmaR,
                0642      I              locMixLayer,
                0643      I              myIter, myThid )
472f570611 Jean*0644 c    U              locK3dRedi
050b4366e6 Jean*0645       ENDIF
                0646 #endif /* ndef GM_EXCLUDE_SUBMESO */
0c49347dc7 Alis*0647
472f570611 Jean*0648 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0649 C     Update components of Redi/GM tensor, second step: multiply by Diffusivity
                0650       DO k=1,Nr
                0651 #ifdef ALLOW_AUTODIFF_TAMC
                0652 C to avoid few minor recomputations (with KAPREDI_CONTROL, KAPGM_CONTROL ?)
                0653        kkey = k + (tkey-1)*Nr
                0654 CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
                0655 CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
                0656 CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
                0657 #endif
                0658        km1 = MAX(k-1,1)
                0659        isopycK = GM_isopycK
                0660      &         *(GM_isoFac1d(km1)+GM_isoFac1d(k))*op5
                0661        bolus_K = GM_background_K
                0662      &         *(GM_bolFac1d(km1)+GM_bolFac1d(k))*op5
                0663        DO j=1-OLy+1,sNy+OLy-1
                0664         DO i=1-OLx+1,sNx+OLx-1
0c12568530 mjlo*0665          Kredi_tmp =
472f570611 Jean*0666 #ifdef GM_USE_K3D_REDI
                0667      &      op5*( locK3dRedi(i,j,km1) + locK3dRedi(i,j,k) )
                0668 #else
                0669      &      isopycK*GM_isoFac2d(i,j,bi,bj)
                0670 #endif
0c12568530 mjlo*0671          Kgm_tmp = Kredi_tmp
472f570611 Jean*0672 #ifdef GM_USE_K3D_GM
                0673      &    + GM_skewflx*op5*( locK3dGM(i,j,km1)+locK3dGM(i,j,k) )
                0674 #else
                0675      &    + GM_skewflx*bolus_K*GM_bolFac2d(i,j,bi,bj)
                0676 #endif
                0677 #ifdef GM_GEOM_VARIABLE_K
                0678      &    + GM_skewflx*GEOM_K3d(i,j,k,bi,bj)
                0679 #endif
                0680          Kwx(i,j,k,bi,bj) = Kgm_tmp*Kwx(i,j,k,bi,bj)
                0681          Kwy(i,j,k,bi,bj) = Kgm_tmp*Kwy(i,j,k,bi,bj)
0c12568530 mjlo*0682          Kwz(i,j,k,bi,bj) = Kredi_tmp*Kwz(i,j,k,bi,bj)
472f570611 Jean*0683         ENDDO
                0684        ENDDO
                0685       ENDDO
                0686
                0687 #ifdef ALLOW_DIAGNOSTICS
                0688       IF ( useDiagnostics .AND. GM_taper_scheme.EQ.'fm07' ) THEN
                0689        CALL DIAGNOSTICS_FILL( hTransLay, 'GM_hTrsL', 0,1,2,bi,bj,myThid)
                0690        CALL DIAGNOSTICS_FILL( baseSlope, 'GM_baseS', 0,1,2,bi,bj,myThid)
                0691        CALL DIAGNOSTICS_FILL(recipLambda,'GM_rLamb', 0,1,2,bi,bj,myThid)
                0692       ENDIF
                0693 #endif /* ALLOW_DIAGNOSTICS */
                0694
                0695 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
a4576c7cde Juli*0696 #if ( defined GM_NON_UNITY_DIAGONAL || defined GM_EXTRA_DIAGONAL )
472f570611 Jean*0697
e2259a1942 Jean*0698 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0699 C-- 2nd  k loop : compute Tensor Coeff. at U point
                0700
a4576c7cde Juli*0701 # ifdef ALLOW_KPP
e2259a1942 Jean*0702       IF ( useKPP ) THEN
ee8a6f4ffb Jean*0703        DO j=1-OLy,sNy+OLy
                0704         DO i=2-OLx,sNx+OLx
e2259a1942 Jean*0705          locMixLayer(i,j) = ( KPPhbl(i-1,j,bi,bj)
                0706      &                      + KPPhbl( i ,j,bi,bj) )*op5
                0707         ENDDO
                0708        ENDDO
                0709       ELSE
a4576c7cde Juli*0710 # else
e2259a1942 Jean*0711       IF ( .TRUE. ) THEN
a4576c7cde Juli*0712 # endif
ee8a6f4ffb Jean*0713        DO j=1-OLy,sNy+OLy
                0714         DO i=2-OLx,sNx+OLx
e2259a1942 Jean*0715          locMixLayer(i,j) = ( hMixLayer(i-1,j,bi,bj)
                0716      &                      + hMixLayer( i ,j,bi,bj) )*op5
                0717         ENDDO
                0718        ENDDO
                0719       ENDIF
ee8a6f4ffb Jean*0720       DO j=1-OLy,sNy+OLy
                0721        DO i=1-OLx,sNx+OLx
e2259a1942 Jean*0722          hTransLay(i,j) =  0.
                0723          baseSlope(i,j) =  0.
                0724          recipLambda(i,j)= 0.
                0725        ENDDO
ee8a6f4ffb Jean*0726        DO i=2-OLx,sNx+OLx
e2259a1942 Jean*0727          hTransLay(i,j) = MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
                0728        ENDDO
                0729       ENDDO
                0730
                0731       DO k=Nr,1,-1
                0732        kp1 = MIN(Nr,k+1)
                0733        maskp1 = 1. _d 0
                0734        IF (k.GE.Nr) maskp1 = 0. _d 0
a4576c7cde Juli*0735 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0736        kkey = k + (tkey-1)*Nr
a4576c7cde Juli*0737 # endif
0c49347dc7 Alis*0738
                0739 C     Gradient of Sigma at U points
ee8a6f4ffb Jean*0740        DO j=1-OLy+1,sNy+OLy-1
                0741         DO i=1-OLx+1,sNx+OLx-1
b0e49a1609 Jean*0742          dSigmaDx(i,j)=sigmaX(i,j,k)
                0743      &                       *_maskW(i,j,k,bi,bj)
                0744          dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k)+sigmaY(i,j+1,k)
                0745      &                       +sigmaY(i-1, j ,k)+sigmaY(i, j ,k)
                0746      &                      )*_maskW(i,j,k,bi,bj)
                0747          dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k )+sigmaR(i,j, k )
                0748      &                      +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
5b172de0d2 Jean*0749      &                      )*_maskW(i,j,k,bi,bj)*gravitySign
b0e49a1609 Jean*0750         ENDDO
0c49347dc7 Alis*0751        ENDDO
                0752
a4576c7cde Juli*0753 # ifdef ALLOW_AUTODIFF_TAMC
2092dbb101 Patr*0754 CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
                0755 CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
5755dbe390 Patr*0756 CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
a4576c7cde Juli*0757 #  ifndef GM_EXCLUDE_FM07_TAP
5755dbe390 Patr*0758 CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
                0759 CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
                0760 CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
a4576c7cde Juli*0761 #  endif
                0762 # endif /* ALLOW_AUTODIFF_TAMC */
2092dbb101 Patr*0763
5b172de0d2 Jean*0764 C     set "rDepth" (= depth from the surface, in rUnit) for 'ldd97' tapering
                0765        IF ( usingZcoords ) THEN
                0766         rDepth = rF(1) - rC(k)
                0767        ELSE
                0768         rDepth = rC(k) - rF(Nr+1)
                0769        ENDIF
0c49347dc7 Alis*0770 C     Calculate slopes for use in tensor, taper and/or clip
b0e49a1609 Jean*0771        CALL GMREDI_SLOPE_LIMIT(
549d1a8d8c Jean*0772      O             SlopeX, SlopeY,
e9fd580bc4 Jean*0773      O             SlopeSqr, taperFct,
e2259a1942 Jean*0774      U             hTransLay, baseSlope, recipLambda,
549d1a8d8c Jean*0775      U             dSigmaDr,
                0776      I             dSigmaDx, dSigmaDy,
5b172de0d2 Jean*0777      I             ldd97_LrhoW, locMixLayer, rDepth, rC,
e2259a1942 Jean*0778      I             kLow_W,
5b172de0d2 Jean*0779      I             1, k, bi, bj, myTime, myIter, myThid )
0c49347dc7 Alis*0780
a4576c7cde Juli*0781 # ifdef ALLOW_AUTODIFF_TAMC
549d1a8d8c Jean*0782 CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
a4576c7cde Juli*0783 # endif /* ALLOW_AUTODIFF_TAMC */
9cb619cfcd Patr*0784
a4576c7cde Juli*0785 # ifdef GM_NON_UNITY_DIAGONAL
b0e49a1609 Jean*0786 c      IF ( GM_nonUnitDiag ) THEN
ee8a6f4ffb Jean*0787         DO j=1-OLy+1,sNy+OLy-1
                0788          DO i=1-OLx+1,sNx+OLx-1
f42e64b3e7 Jean*0789           Kux(i,j,k,bi,bj) =
472f570611 Jean*0790 #  ifdef GM_USE_K3D_REDI
                0791      &      op5*( locK3dRedi(i-1,j,k) + locK3dRedi(i,j,k) )
a4576c7cde Juli*0792 #  else
472f570611 Jean*0793      &      GM_isopycK*GM_isoFac1d(k)
f6de204bec Jean*0794      &        *op5*(GM_isoFac2d(i-1,j,bi,bj)+GM_isoFac2d(i,j,bi,bj))
a4576c7cde Juli*0795 #  endif
472f570611 Jean*0796      &     *taperFct(i,j)
b6b11b9b2f Patr*0797          ENDDO
                0798         ENDDO
a4576c7cde Juli*0799 #  if ( defined ALLOW_AUTODIFF_TAMC && defined GM_EXCLUDE_CLIPPING )
2092dbb101 Patr*0800 CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
a4576c7cde Juli*0801 #  endif
ee8a6f4ffb Jean*0802         DO j=1-OLy+1,sNy+OLy-1
                0803          DO i=1-OLx+1,sNx+OLx-1
f42e64b3e7 Jean*0804           Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
                0805          ENDDO
                0806         ENDDO
b0e49a1609 Jean*0807 c      ENDIF
a4576c7cde Juli*0808 # endif /* GM_NON_UNITY_DIAGONAL */
f42e64b3e7 Jean*0809
a4576c7cde Juli*0810 # ifdef GM_EXTRA_DIAGONAL
2092dbb101 Patr*0811
a4576c7cde Juli*0812 #  ifdef ALLOW_AUTODIFF_TAMC
2092dbb101 Patr*0813 CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
a4576c7cde Juli*0814 #  endif /* ALLOW_AUTODIFF_TAMC */
e2259a1942 Jean*0815        IF ( GM_ExtraDiag ) THEN
ee8a6f4ffb Jean*0816         DO j=1-OLy+1,sNy+OLy-1
                0817          DO i=1-OLx+1,sNx+OLx-1
5b172de0d2 Jean*0818           Kuz(i,j,k,bi,bj) = -gravitySign*
472f570611 Jean*0819 #  ifdef GM_USE_K3D_REDI
                0820      &     ( op5*( locK3dRedi(i-1,j,k) + locK3dRedi(i,j,k) )
a4576c7cde Juli*0821 #  else
f6de204bec Jean*0822      &     ( GM_isopycK*GM_isoFac1d(k)
                0823      &        *op5*(GM_isoFac2d(i-1,j,bi,bj)+GM_isoFac2d(i,j,bi,bj))
a4576c7cde Juli*0824 #  endif
472f570611 Jean*0825 #  ifdef GM_USE_K3D_GM
                0826      &     - GM_skewflx*op5*( locK3dGM(i-1,j,k) + locK3dGM(i,j,k) )
a4576c7cde Juli*0827 #  else
f6de204bec Jean*0828      &     - GM_skewflx*GM_background_K*GM_bolFac1d(k)
                0829      &        *op5*(GM_bolFac2d(i-1,j,bi,bj)+GM_bolFac2d(i,j,bi,bj))
a4576c7cde Juli*0830 #  endif
                0831 #  ifdef GM_GEOM_VARIABLE_K
                0832      &     - GM_skewflx*op25*( ( GEOM_K3d(i-1,j, k, bi,bj)
                0833      &                         + GEOM_K3d( i, j, k, bi,bj) )
                0834      &                       + ( GEOM_K3d(i-1,j,kp1,bi,bj)
                0835      &                         + GEOM_K3d( i, j,kp1,bi,bj) ) )
                0836 #  endif
f42e64b3e7 Jean*0837      &     )*SlopeX(i,j)*taperFct(i,j)
                0838          ENDDO
                0839         ENDDO
796b5e35f7 Jean*0840 c      ELSE
                0841 c       DO j=1-OLy+1,sNy+OLy-1
                0842 c        DO i=1-OLx+1,sNx+OLx-1
                0843 c         Kuz(i,j,k,bi,bj) = 0. _d 0
                0844 c        ENDDO
                0845 c       ENDDO
b0e49a1609 Jean*0846        ENDIF
a4576c7cde Juli*0847 # endif /* GM_EXTRA_DIAGONAL */
0c49347dc7 Alis*0848
a4576c7cde Juli*0849 # ifdef ALLOW_DIAGNOSTICS
b0e49a1609 Jean*0850        IF (doDiagRediFlx) THEN
066e0d5e64 Jean*0851         km1 = MAX(k-1,1)
                0852         DO j=1,sNy
                0853          DO i=1,sNx+1
                0854 C         store in tmp1k Kuz_Redi
5b172de0d2 Jean*0855           tmp1k(i,j) = -gravitySign*
472f570611 Jean*0856 #  ifdef GM_USE_K3D_REDI
                0857      &      op5*( locK3dRedi(i-1,j,k)+ locK3dRedi(i,j,k) )
a4576c7cde Juli*0858 #  else
472f570611 Jean*0859      &      GM_isopycK*GM_isoFac1d(k)
                0860      &       *op5*(GM_isoFac2d(i-1,j,bi,bj)+GM_isoFac2d(i,j,bi,bj))
a4576c7cde Juli*0861 #  endif
472f570611 Jean*0862      &     *SlopeX(i,j)*taperFct(i,j)
066e0d5e64 Jean*0863          ENDDO
                0864         ENDDO
                0865         DO j=1,sNy
                0866          DO i=1,sNx+1
                0867 C-        Vertical gradients interpolated to U points
                0868           dTdz = (
                0869      &     +recip_drC(k)*
8233d0ceb9 Jean*0870      &       ( maskC(i-1,j,km1,bi,bj)*maskC(i-1,j,k,bi,bj)*
066e0d5e64 Jean*0871      &           (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
8233d0ceb9 Jean*0872      &        +maskC( i ,j,km1,bi,bj)*maskC( i ,j,k,bi,bj)*
066e0d5e64 Jean*0873      &           (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
                0874      &       )
                0875      &     +recip_drC(kp1)*
8233d0ceb9 Jean*0876      &       ( maskC(i-1,j,k,bi,bj)*maskC(i-1,j,kp1,bi,bj)*
066e0d5e64 Jean*0877      &           (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
8233d0ceb9 Jean*0878      &        +maskC( i ,j,k,bi,bj)*maskC( i ,j,kp1,bi,bj)*
066e0d5e64 Jean*0879      &           (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
                0880      &       )      ) * 0.25 _d 0
a67797e4f0 Jean*0881            tmp1k(i,j) = dyG(i,j,bi,bj) * deepFacC(k)
                0882      &                * drF(k) * _hFacW(i,j,k,bi,bj)
066e0d5e64 Jean*0883      &                * tmp1k(i,j) * dTdz
                0884          ENDDO
                0885         ENDDO
                0886         CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
b0e49a1609 Jean*0887        ENDIF
a4576c7cde Juli*0888 # endif /* ALLOW_DIAGNOSTICS */
066e0d5e64 Jean*0889
e2259a1942 Jean*0890 C-- end 2nd  loop on vertical level index k
                0891       ENDDO
                0892
                0893 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0894 C-- 3rd  k loop : compute Tensor Coeff. at V point
                0895
a4576c7cde Juli*0896 # ifdef ALLOW_KPP
e2259a1942 Jean*0897       IF ( useKPP ) THEN
ee8a6f4ffb Jean*0898        DO j=2-OLy,sNy+OLy
                0899         DO i=1-OLx,sNx+OLx
e2259a1942 Jean*0900          locMixLayer(i,j) = ( KPPhbl(i,j-1,bi,bj)
                0901      &                      + KPPhbl(i, j ,bi,bj) )*op5
                0902         ENDDO
                0903        ENDDO
                0904       ELSE
a4576c7cde Juli*0905 # else
e2259a1942 Jean*0906       IF ( .TRUE. ) THEN
a4576c7cde Juli*0907 # endif
ee8a6f4ffb Jean*0908        DO j=2-OLy,sNy+OLy
                0909         DO i=1-OLx,sNx+OLx
e2259a1942 Jean*0910          locMixLayer(i,j) = ( hMixLayer(i,j-1,bi,bj)
                0911      &                      + hMixLayer(i, j ,bi,bj) )*op5
                0912         ENDDO
                0913        ENDDO
                0914       ENDIF
ee8a6f4ffb Jean*0915       DO j=1-OLy,sNy+OLy
                0916        DO i=1-OLx,sNx+OLx
e2259a1942 Jean*0917          hTransLay(i,j) =  0.
                0918          baseSlope(i,j) =  0.
                0919          recipLambda(i,j)= 0.
                0920        ENDDO
                0921       ENDDO
ee8a6f4ffb Jean*0922       DO j=2-OLy,sNy+OLy
                0923        DO i=1-OLx,sNx+OLx
e2259a1942 Jean*0924          hTransLay(i,j) = MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
                0925        ENDDO
                0926       ENDDO
                0927
0c49347dc7 Alis*0928 C     Gradient of Sigma at V points
e2259a1942 Jean*0929       DO k=Nr,1,-1
                0930        kp1 = MIN(Nr,k+1)
                0931        maskp1 = 1. _d 0
                0932        IF (k.GE.Nr) maskp1 = 0. _d 0
a4576c7cde Juli*0933 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0934        kkey = k + (tkey-1)*Nr
a4576c7cde Juli*0935 # endif
e2259a1942 Jean*0936
ee8a6f4ffb Jean*0937        DO j=1-OLy+1,sNy+OLy-1
                0938         DO i=1-OLx+1,sNx+OLx-1
b0e49a1609 Jean*0939          dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
                0940      &                       +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k)
                0941      &                      )*_maskS(i,j,k,bi,bj)
                0942          dSigmaDy(i,j)=sigmaY(i,j,k)
                0943      &                       *_maskS(i,j,k,bi,bj)
                0944          dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
                0945      &                      +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
5b172de0d2 Jean*0946      &                      )*_maskS(i,j,k,bi,bj)*gravitySign
b0e49a1609 Jean*0947         ENDDO
0c49347dc7 Alis*0948        ENDDO
                0949
a4576c7cde Juli*0950 # ifdef ALLOW_AUTODIFF_TAMC
2092dbb101 Patr*0951 CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
                0952 CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
5755dbe390 Patr*0953 CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
a4576c7cde Juli*0954 #  ifndef GM_EXCLUDE_FM07_TAP
5755dbe390 Patr*0955 CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
                0956 CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
                0957 CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
a4576c7cde Juli*0958 #  endif
                0959 # endif /* ALLOW_AUTODIFF_TAMC */
2092dbb101 Patr*0960
5b172de0d2 Jean*0961 C     set "rDepth" (= depth from the surface, in rUnit) for 'ldd97' tapering
                0962        IF ( usingZcoords ) THEN
                0963         rDepth = rF(1) - rC(k)
                0964        ELSE
                0965         rDepth = rC(k) - rF(Nr+1)
                0966        ENDIF
0c49347dc7 Alis*0967 C     Calculate slopes for use in tensor, taper and/or clip
b0e49a1609 Jean*0968        CALL GMREDI_SLOPE_LIMIT(
549d1a8d8c Jean*0969      O             SlopeX, SlopeY,
e9fd580bc4 Jean*0970      O             SlopeSqr, taperFct,
e2259a1942 Jean*0971      U             hTransLay, baseSlope, recipLambda,
549d1a8d8c Jean*0972      U             dSigmaDr,
                0973      I             dSigmaDx, dSigmaDy,
5b172de0d2 Jean*0974      I             ldd97_LrhoS, locMixLayer, rDepth, rC,
e2259a1942 Jean*0975      I             kLow_S,
5b172de0d2 Jean*0976      I             2, k, bi, bj, myTime, myIter, myThid )
0c49347dc7 Alis*0977
a4576c7cde Juli*0978 # ifdef ALLOW_AUTODIFF_TAMC
361543e926 Jean*0979 CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
a4576c7cde Juli*0980 # endif /* ALLOW_AUTODIFF_TAMC */
9cb619cfcd Patr*0981
a4576c7cde Juli*0982 # ifdef GM_NON_UNITY_DIAGONAL
b0e49a1609 Jean*0983 c      IF ( GM_nonUnitDiag ) THEN
ee8a6f4ffb Jean*0984         DO j=1-OLy+1,sNy+OLy-1
                0985          DO i=1-OLx+1,sNx+OLx-1
f42e64b3e7 Jean*0986           Kvy(i,j,k,bi,bj) =
472f570611 Jean*0987 #  ifdef GM_USE_K3D_REDI
                0988      &      op5*( locK3dRedi(i,j-1,k) + locK3dRedi(i,j,k) )
a4576c7cde Juli*0989 #  else
472f570611 Jean*0990      &      GM_isopycK*GM_isoFac1d(k)
f6de204bec Jean*0991      &        *op5*(GM_isoFac2d(i,j-1,bi,bj)+GM_isoFac2d(i,j,bi,bj))
a4576c7cde Juli*0992 #  endif
472f570611 Jean*0993      &     *taperFct(i,j)
b6b11b9b2f Patr*0994          ENDDO
                0995         ENDDO
a4576c7cde Juli*0996 #  if ( defined ALLOW_AUTODIFF_TAMC && defined GM_EXCLUDE_CLIPPING )
2092dbb101 Patr*0997 CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
a4576c7cde Juli*0998 #  endif
ee8a6f4ffb Jean*0999         DO j=1-OLy+1,sNy+OLy-1
                1000          DO i=1-OLx+1,sNx+OLx-1
f42e64b3e7 Jean*1001           Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
                1002          ENDDO
                1003         ENDDO
b0e49a1609 Jean*1004 c      ENDIF
a4576c7cde Juli*1005 # endif /* GM_NON_UNITY_DIAGONAL */
f42e64b3e7 Jean*1006
a4576c7cde Juli*1007 # ifdef GM_EXTRA_DIAGONAL
2092dbb101 Patr*1008
a4576c7cde Juli*1009 #  ifdef ALLOW_AUTODIFF_TAMC
2092dbb101 Patr*1010 CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
a4576c7cde Juli*1011 #  endif /* ALLOW_AUTODIFF_TAMC */
e2259a1942 Jean*1012        IF ( GM_ExtraDiag ) THEN
ee8a6f4ffb Jean*1013         DO j=1-OLy+1,sNy+OLy-1
                1014          DO i=1-OLx+1,sNx+OLx-1
5b172de0d2 Jean*1015           Kvz(i,j,k,bi,bj) = -gravitySign*
472f570611 Jean*1016 #  ifdef GM_USE_K3D_REDI
                1017      &     ( op5*( locK3dRedi(i,j-1,k) + locK3dRedi(i,j,k) )
a4576c7cde Juli*1018 #  else
f6de204bec Jean*1019      &     ( GM_isopycK*GM_isoFac1d(k)
                1020      &        *op5*(GM_isoFac2d(i,j-1,bi,bj)+GM_isoFac2d(i,j,bi,bj))
a4576c7cde Juli*1021 #  endif
472f570611 Jean*1022 #  ifdef GM_USE_K3D_GM
                1023      &     - GM_skewflx*op5*( locK3dGM(i,j-1,k) + locK3dGM(i,j,k) )
a4576c7cde Juli*1024 #  else
f6de204bec Jean*1025      &     - GM_skewflx*GM_background_K*GM_bolFac1d(k)
                1026      &        *op5*(GM_bolFac2d(i,j-1,bi,bj)+GM_bolFac2d(i,j,bi,bj))
a4576c7cde Juli*1027 #  endif
                1028 #  ifdef GM_GEOM_VARIABLE_K
                1029      &     - GM_skewflx*op25*( ( GEOM_K3d(i,j-1, k, bi,bj)
                1030      &                         + GEOM_K3d(i, j,  k, bi,bj) )
                1031      &                       + ( GEOM_K3d(i,j-1,kp1,bi,bj)
                1032      &                         + GEOM_K3d(i, j, kp1,bi,bj) ) )
                1033 #  endif
f42e64b3e7 Jean*1034      &     )*SlopeY(i,j)*taperFct(i,j)
                1035          ENDDO
                1036         ENDDO
796b5e35f7 Jean*1037 c      ELSE
                1038 c       DO j=1-OLy+1,sNy+OLy-1
                1039 c        DO i=1-OLx+1,sNx+OLx-1
                1040 c         Kvz(i,j,k,bi,bj) = 0. _d 0
                1041 c        ENDDO
                1042 c       ENDDO
b0e49a1609 Jean*1043        ENDIF
a4576c7cde Juli*1044 # endif /* GM_EXTRA_DIAGONAL */
f42e64b3e7 Jean*1045
a4576c7cde Juli*1046 # ifdef ALLOW_DIAGNOSTICS
b0e49a1609 Jean*1047        IF (doDiagRediFlx) THEN
81880fdab4 Davi*1048         km1 = MAX(k-1,1)
066e0d5e64 Jean*1049         DO j=1,sNy+1
                1050          DO i=1,sNx
                1051 C         store in tmp1k Kvz_Redi
5b172de0d2 Jean*1052           tmp1k(i,j) = -gravitySign*
472f570611 Jean*1053 #  ifdef GM_USE_K3D_REDI
                1054      &      op5*( locK3dRedi(i,j-1,k) + locK3dRedi(i,j,k) )
a4576c7cde Juli*1055 #  else
472f570611 Jean*1056      &      GM_isopycK*GM_isoFac1d(k)
                1057      &       *op5*(GM_isoFac2d(i,j-1,bi,bj)+GM_isoFac2d(i,j,bi,bj))
a4576c7cde Juli*1058 #  endif
472f570611 Jean*1059      &     *SlopeY(i,j)*taperFct(i,j)
066e0d5e64 Jean*1060          ENDDO
                1061         ENDDO
                1062         DO j=1,sNy+1
                1063          DO i=1,sNx
                1064 C-        Vertical gradients interpolated to U points
                1065           dTdz = (
                1066      &     +recip_drC(k)*
8233d0ceb9 Jean*1067      &       ( maskC(i,j-1,km1,bi,bj)*maskC(i,j-1,k,bi,bj)*
066e0d5e64 Jean*1068      &           (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
8233d0ceb9 Jean*1069      &        +maskC(i, j ,km1,bi,bj)*maskC(i, j ,k,bi,bj)*
066e0d5e64 Jean*1070      &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
                1071      &       )
                1072      &     +recip_drC(kp1)*
8233d0ceb9 Jean*1073      &       ( maskC(i,j-1,kp1,bi,bj)*maskC(i,j-1,k,bi,bj)*
066e0d5e64 Jean*1074      &           (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
8233d0ceb9 Jean*1075      &        +maskC(i, j ,kp1,bi,bj)*maskC(i, j ,k,bi,bj)*
066e0d5e64 Jean*1076      &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
                1077      &       )      ) * 0.25 _d 0
a67797e4f0 Jean*1078            tmp1k(i,j) = dxG(i,j,bi,bj) * deepFacC(k)
                1079      &                * drF(k) * _hFacS(i,j,k,bi,bj)
066e0d5e64 Jean*1080      &                * tmp1k(i,j) * dTdz
                1081          ENDDO
                1082         ENDDO
                1083         CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
b0e49a1609 Jean*1084        ENDIF
a4576c7cde Juli*1085 # endif /* ALLOW_DIAGNOSTICS */
066e0d5e64 Jean*1086
e2259a1942 Jean*1087 C-- end 3rd  loop on vertical level index k
f42e64b3e7 Jean*1088       ENDDO
0c49347dc7 Alis*1089
e2259a1942 Jean*1090 #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
472f570611 Jean*1091 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
e2259a1942 Jean*1092
796b5e35f7 Jean*1093 #ifndef GM_NON_UNITY_DIAGONAL
f5509be190 Mart*1094 C--   keep a simplified setting here (used to be inside gmredi_x/ytransport)
                1095 C     for backeard compatibility + for trying to generate a simpler adjoint code
796b5e35f7 Jean*1096       DO k=1,Nr
                1097        DO j=1-OLy+1,sNy+OLy-1
                1098         DO i=1-OLx+1,sNx+OLx-1
472f570611 Jean*1099           Kux(i,j,k,bi,bj) =
                1100 # ifdef GM_USE_K3D_REDI
                1101      &      op5*( locK3dRedi(i-1,j,k) + locK3dRedi(i,j,k) )
a4576c7cde Juli*1102 # else
472f570611 Jean*1103      &      GM_isopycK
a4576c7cde Juli*1104 # endif
796b5e35f7 Jean*1105         ENDDO
                1106        ENDDO
                1107        DO j=1-OLy+1,sNy+OLy-1
                1108         DO i=1-OLx+1,sNx+OLx-1
472f570611 Jean*1109           Kvy(i,j,k,bi,bj) =
                1110 # ifdef GM_USE_K3D_REDI
                1111      &      op5*( locK3dRedi(i,j-1,k) + locK3dRedi(i,j,k) )
a4576c7cde Juli*1112 # else
472f570611 Jean*1113      &      GM_isopycK
a4576c7cde Juli*1114 # endif
796b5e35f7 Jean*1115         ENDDO
                1116        ENDDO
                1117       ENDDO
                1118 #endif /* ndef GM_NON_UNITY_DIAGONAL */
                1119
d29d98918f Jean*1120 #ifdef ALLOW_DIAGNOSTICS
                1121       IF ( useDiagnostics ) THEN
472f570611 Jean*1122         CALL GMREDI_DIAGNOSTICS_FILL( bi, bj, myThid )
d29d98918f Jean*1123       ENDIF
                1124 #endif /* ALLOW_DIAGNOSTICS */
                1125
0c49347dc7 Alis*1126 #endif /* ALLOW_GMREDI */
                1127
                1128       RETURN
                1129       END
b58589f5c2 Patr*1130
8d8488bf10 Jean*1131 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
b58589f5c2 Patr*1132
3a15bf3a95 Jean*1133 CBOP
                1134 C     !ROUTINE: GMREDI_CALC_TENSOR_DUMMY
                1135 C     !INTERFACE:
b58589f5c2 Patr*1136       SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
e2259a1942 Jean*1137      I             iMin, iMax, jMin, jMax,
b58589f5c2 Patr*1138      I             sigmaX, sigmaY, sigmaR,
e2259a1942 Jean*1139      I             bi, bj, myTime, myIter, myThid )
3a15bf3a95 Jean*1140
                1141 C     !DESCRIPTION: \bv
                1142 C     *==========================================================*
                1143 C     | SUBROUTINE GMREDI_CALC_TENSOR_DUMMY
                1144 C     | o Calculate tensor elements for GM/Redi tensor.
                1145 C     *==========================================================*
                1146 C     \ev
                1147
                1148 C     !USES:
b58589f5c2 Patr*1149       IMPLICIT NONE
                1150
                1151 C     == Global variables ==
                1152 #include "SIZE.h"
                1153 #include "EEPARAMS.h"
                1154 #include "GMREDI.h"
                1155
3a15bf3a95 Jean*1156 C     !INPUT/OUTPUT PARAMETERS:
ee8a6f4ffb Jean*1157       _RL sigmaX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                1158       _RL sigmaY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                1159       _RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
e2259a1942 Jean*1160       INTEGER iMin,iMax,jMin,jMax
                1161       INTEGER bi, bj
                1162       _RL     myTime
                1163       INTEGER myIter
b58589f5c2 Patr*1164       INTEGER myThid
3a15bf3a95 Jean*1165 CEOP
b58589f5c2 Patr*1166
                1167 #ifdef ALLOW_GMREDI
3a15bf3a95 Jean*1168 C     !LOCAL VARIABLES:
e2259a1942 Jean*1169       INTEGER i, j, k
                1170
f42e64b3e7 Jean*1171       DO k=1,Nr
ee8a6f4ffb Jean*1172        DO j=1-OLy+1,sNy+OLy-1
                1173         DO i=1-OLx+1,sNx+OLx-1
796b5e35f7 Jean*1174          Kwx(i,j,k,bi,bj) = 0. _d 0
                1175          Kwy(i,j,k,bi,bj) = 0. _d 0
                1176          Kwz(i,j,k,bi,bj) = 0. _d 0
                1177          Kux(i,j,k,bi,bj) = 0. _d 0
                1178          Kvy(i,j,k,bi,bj) = 0. _d 0
                1179 # ifdef GM_EXTRA_DIAGONAL
                1180          Kuz(i,j,k,bi,bj) = 0. _d 0
                1181          Kvz(i,j,k,bi,bj) = 0. _d 0
                1182 # endif
f42e64b3e7 Jean*1183         ENDDO
b58589f5c2 Patr*1184        ENDDO
                1185       ENDDO
                1186 #endif /* ALLOW_GMREDI */
                1187
f42e64b3e7 Jean*1188       RETURN
                1189       END