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