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