pkg/gmredi/gmredi_slope_limit.F

0c49347dc7 Alis*0001 #include "GMREDI_OPTIONS.h"
14e0496834 Jean*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
0c49347dc7 Alis*0005
e2259a1942 Jean*0006 CBOP
                0007 C     !ROUTINE: GMREDI_SLOPE_LIMIT
                0008 C     !INTERFACE:
0c49347dc7 Alis*0009       SUBROUTINE GMREDI_SLOPE_LIMIT(
549d1a8d8c Jean*0010      O             SlopeX, SlopeY,
e9fd580bc4 Jean*0011      O             SlopeSqr, taperFct,
e2259a1942 Jean*0012      U             hTransLay, baseSlope, recipLambda,
549d1a8d8c Jean*0013      U             dSigmaDr,
                0014      I             dSigmaDx, dSigmaDy,
5b172de0d2 Jean*0015      I             Lrho, hMixLay, rDepth, depthZ, kLow,
                0016      I             kPos, k, bi, bj, myTime, myIter, myThid )
e2259a1942 Jean*0017
                0018 C     !DESCRIPTION: \bv
                0019 C     *==========================================================*
                0020 C     | SUBROUTINE GMREDI_SLOPE_LIMIT
                0021 C     | o Calculate slopes for use in GM/Redi tensor
                0022 C     *==========================================================*
                0023 C     | On entry:
5b172de0d2 Jean*0024 C     |            dSigmaDr     contains the downward d/dz Sigma
e2259a1942 Jean*0025 C     |            dSigmaDx/Dy  contains X/Y gradients of sigma
                0026 C     |            Lrho
                0027 C     |            hMixLay      mixed layer depth (> 0)
5b172de0d2 Jean*0028 C     |            rDepth       depth (> 0) in r-Unit from the surface
                0029 C     |            depthZ       contains the depth (< 0 !) [m]
e2259a1942 Jean*0030 C    U             hTransLay    transition layer depth (> 0)
                0031 C    U             baseSlope, recipLambda,
                0032 C     | On exit:
5b172de0d2 Jean*0033 C     |            dSigmaDr     contains the effective downward dSig/dz
e2259a1942 Jean*0034 C     |            SlopeX/Y     contains X/Y slopes
                0035 C     |            SlopeSqr     contains Sx^2+Sy^2
                0036 C     |            taperFct     contains tapering funct. value ;
                0037 C     |                         = 1 when using no tapering
                0038 C    U             hTransLay    transition layer depth (> 0)
                0039 C    U             baseSlope, recipLambda
                0040 C     *==========================================================*
                0041 C     \ev
                0042
                0043 C     !USES:
0c49347dc7 Alis*0044       IMPLICIT NONE
                0045
                0046 C     == Global variables ==
                0047 #include "SIZE.h"
b6b11b9b2f Patr*0048 #include "GRID.h"
0c49347dc7 Alis*0049 #include "EEPARAMS.h"
                0050 #include "GMREDI.h"
                0051 #include "PARAMS.h"
e673dad091 Patr*0052 #ifdef ALLOW_AUTODIFF_TAMC
                0053 #include "tamc.h"
                0054 #endif /* ALLOW_AUTODIFF_TAMC */
0c49347dc7 Alis*0055
549d1a8d8c Jean*0056 C     !INPUT PARAMETERS:
e2259a1942 Jean*0057 C     dSigmaDx :: zonal      gradient of density
                0058 C     dSigmaDy :: meridional gradient of density
                0059 C     Lrho     ::
                0060 C     hMixLay  :: mixed layer thickness (> 0) [m]
5b172de0d2 Jean*0061 C     rDepth   :: depth (> 0) in r-Unit from the surface
e2259a1942 Jean*0062 C     depthZ   :: model discretized depth (< 0) [m]
                0063 C     kLow     :: bottom level index 2-D array
5b172de0d2 Jean*0064 C     kPos     :: grid-cell location: 1,2,3 : at U,V,W location
e2259a1942 Jean*0065 C     k        :: level index
                0066 C     bi, bj   :: tile indices
                0067 C     myTime   :: time in simulation
                0068 C     myIter   :: iteration number in simulation
                0069 C     myThid   :: My Thread Id. number
549d1a8d8c Jean*0070 C     !OUTPUT PARAMETERS:
e2259a1942 Jean*0071 C     SlopeX      :: isopycnal slope in X direction
                0072 C     SlopeY      :: isopycnal slope in Y direction
                0073 C     SlopeSqr    :: square of isopycnal slope
                0074 C     taperFct    :: tapering function
                0075 C     hTransLay   :: depth of the base of the transition layer (> 0) [m]
                0076 C     baseSlope   :: slope at the the base of the transition layer
                0077 C     recipLambda :: Slope vertical gradient at Trans. Layer Base (=recip.Lambda)
5b172de0d2 Jean*0078 C     dSigmaDr    :: downward gradient of neutral density
14e0496834 Jean*0079       _RL SlopeX     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0080       _RL SlopeY     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0081       _RL SlopeSqr   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0082       _RL taperFct   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0083       _RL hTransLay  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0084       _RL baseSlope  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085       _RL recipLambda(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0086       _RL dSigmaDr   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0087       _RL dSigmaDx   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0088       _RL dSigmaDy   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0089       _RL Lrho       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0090       _RL hMixLay    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
5b172de0d2 Jean*0091       _RL rDepth
e2259a1942 Jean*0092       _RS depthZ(*)
14e0496834 Jean*0093       INTEGER kLow   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
5b172de0d2 Jean*0094       INTEGER kPos, k, bi, bj
e2259a1942 Jean*0095       _RL     myTime
                0096       INTEGER myIter
                0097       INTEGER myThid
0c49347dc7 Alis*0098
                0099 #ifdef ALLOW_GMREDI
e2259a1942 Jean*0100 C     !LOCAL VARIABLES:
0c49347dc7 Alis*0101       _RL fpi
413f546150 Jean*0102       PARAMETER( fpi = PI )
9f5240b52a Jean*0103       INTEGER i, j
                0104       _RL f1, Smod, f2, Rnondim
                0105       _RL maxSlopeSqr
5b172de0d2 Jean*0106       _RL convSlopeUnit, loc_rMaxSlope
9f5240b52a Jean*0107       _RL dSigmMod   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0108       _RL dRdSigmaLtd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0109       _RL tmpFld     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0110 #ifndef GM_EXCLUDE_FM07_TAP
25ef785e26 Jean*0111       _RL dTransLay, rLambMin, DoverLamb
ac0d9d59a8 Jean*0112       _RL taperFctLoc, taperFctHat
25ef785e26 Jean*0113       _RL minTransLay
9f5240b52a Jean*0114       _RL SlopeMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0115       _RL locVar  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0116 #endif
                0117 #ifdef GMREDI_WITH_STABLE_ADJOINT
                0118       _RL slopeSqTmp, slopeTmp, slopeMax
5a6997640b Mart*0119 c# ifdef ALLOW_AUTODIFF_TAMC
                0120 c      INTEGER ikey
                0121 c# endif
7c50f07931 Mart*0122 #endif
df4facdcff Jean*0123 C-   need to put this one in GM namelist :
                0124       _RL GM_bigSlope
                0125       GM_bigSlope = 1. _d +02
5b172de0d2 Jean*0126 CEOP
                0127
                0128       IF (kPos.EQ.3 ) THEN
                0129        GM_bigSlope   = GM_bigSlope*wUnit2rVel(k)
                0130        maxSlopeSqr   = GM_maxSlope*GM_maxSlope
                0131      &                *wUnit2rVel(k)*wUnit2rVel(k)
                0132        convSlopeUnit = rVel2wUnit(k)
                0133       ELSE
                0134        GM_bigSlope   = GM_bigSlope*z2rUnit(k)
                0135        maxSlopeSqr   = GM_maxSlope*GM_maxSlope
                0136      &                *z2rUnit(k)*z2rUnit(k)
                0137        convSlopeUnit = rUnit2z(k)
                0138       ENDIF
                0139       loc_rMaxSlope  = GM_rMaxSlope*convSlopeUnit
df4facdcff Jean*0140
e673dad091 Patr*0141 #ifdef ALLOW_AUTODIFF_TAMC
5b172de0d2 Jean*0142 C     TAF thinks for some reason that rDepth is an active variable.
9d44c9cacc Mart*0143 C     While this does not make the adjoint code wrong, the resulting
                0144 C     code inhibits vectorization in some cases so we tell TAF here
5b172de0d2 Jean*0145 C     that rDepth is actually a passive grid variable that needs no adjoint
                0146 CADJ PASSIVE rDepth
113b21cd45 Mart*0147 C     Without this store directive, TAF generates an extra field anyway
                0148 C     so we do it here explicitly with a local tape and live without the
                0149 C     corresponding warning.
                0150 CADJ INIT loctape_gm = COMMON, 1
7448700841 Mart*0151 C     Note that this directive is not necessary in any verification
                0152 C     experiment, so we leave it commented out for now.
                0153 cCADJ STORE dSigmaDr = loctape_gm
5a6997640b Mart*0154 C
                0155 C     Alternatively one can use a global tape, to be defined in
                0156 C     the_main_loop.F. This option is still here but commented out with
                0157 C     "c" (also in the_main_loop.F), in case we decide to go for larger
                0158 C     memory overhead and fewer recomputations
                0159 c      ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0160 c      ikey = kPos + (k - 1  + (ikey - 1)*Nr)*3
                0161 cCADJ STORE dSigmaDr = comlev1_gmredi_slope, key=ikey, byte=isbyte
e673dad091 Patr*0162 #endif /* ALLOW_AUTODIFF_TAMC */
                0163
14e0496834 Jean*0164       DO j=1-OLy+1,sNy+OLy-1
                0165        DO i=1-OLx+1,sNx+OLx-1
e2259a1942 Jean*0166         dSigmMod(i,j)    = 0. _d 0
ac0d9d59a8 Jean*0167         tmpFld(i,j)      = 0. _d 0
b6b11b9b2f Patr*0168        ENDDO
                0169       ENDDO
                0170
e9fd580bc4 Jean*0171       IF (GM_taper_scheme.EQ.'orig' .OR.
                0172      &    GM_taper_scheme.EQ.'clipping') THEN
0c49347dc7 Alis*0173
bd3f833a36 Jean*0174 #ifdef GM_EXCLUDE_CLIPPING
                0175
                0176         STOP 'Need to compile without "#define GM_EXCLUDE_CLIPPING"'
                0177
                0178 #else  /* GM_EXCLUDE_CLIPPING */
2092dbb101 Patr*0179
0c49347dc7 Alis*0180 C-      Original implementation in mitgcmuv
                0181 C       (this turns out to be the same as Cox slope clipping)
e9fd580bc4 Jean*0182
                0183 C-      Cox 1987 "Slope clipping"
14e0496834 Jean*0184         DO j=1-OLy+1,sNy+OLy-1
                0185          DO i=1-OLx+1,sNx+OLx-1
ac0d9d59a8 Jean*0186           tmpFld(i,j) = dSigmaDx(i,j)*dSigmaDx(i,j)
e2259a1942 Jean*0187      &                + dSigmaDy(i,j)*dSigmaDy(i,j)
5b172de0d2 Jean*0188           IF ( tmpFld(i,j) .EQ. zeroRL ) THEN
e2259a1942 Jean*0189            dSigmMod(i,j) = 0. _d 0
b6b11b9b2f Patr*0190           ELSE
ac0d9d59a8 Jean*0191            dSigmMod(i,j) = SQRT( tmpFld(i,j) )
b6b11b9b2f Patr*0192           ENDIF
e673dad091 Patr*0193          ENDDO
                0194         ENDDO
0c49347dc7 Alis*0195
14e0496834 Jean*0196         DO j=1-OLy+1,sNy+OLy-1
                0197          DO i=1-OLx+1,sNx+OLx-1
5b172de0d2 Jean*0198           IF ( dSigmMod(i,j) .NE. zeroRL ) THEN
                0199            tmpFld(i,j) = dSigmMod(i,j)*loc_rMaxSlope
                0200            IF ( dSigmaDr(i,j) .LE. tmpFld(i,j) )
ac0d9d59a8 Jean*0201      &          dSigmaDr(i,j) = tmpFld(i,j)
b6b11b9b2f Patr*0202           ENDIF
e673dad091 Patr*0203          ENDDO
                0204         ENDDO
                0205
14e0496834 Jean*0206         DO j=1-OLy+1,sNy+OLy-1
                0207          DO i=1-OLx+1,sNx+OLx-1
5b172de0d2 Jean*0208           IF ( dSigmMod(i,j) .EQ. zeroRL ) THEN
b6b11b9b2f Patr*0209            SlopeX(i,j) = 0. _d 0
                0210            SlopeY(i,j) = 0. _d 0
                0211           ELSE
e2259a1942 Jean*0212            dRdSigmaLtd(i,j) = 1. _d 0/( dSigmaDr(i,j) )
5b172de0d2 Jean*0213            SlopeX(i,j) = dSigmaDx(i,j)*dRdSigmaLtd(i,j)
                0214            SlopeY(i,j) = dSigmaDy(i,j)*dRdSigmaLtd(i,j)
b6b11b9b2f Patr*0215           ENDIF
e673dad091 Patr*0216          ENDDO
                0217         ENDDO
                0218
14e0496834 Jean*0219         DO j=1-OLy+1,sNy+OLy-1
                0220          DO i=1-OLx+1,sNx+OLx-1
8233d0ceb9 Jean*0221           SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j)
                0222      &                  + SlopeY(i,j)*SlopeY(i,j)
                0223           taperFct(i,j) = 1. _d 0
0c49347dc7 Alis*0224          ENDDO
                0225         ENDDO
                0226
bd3f833a36 Jean*0227 #endif /* GM_EXCLUDE_CLIPPING */
2092dbb101 Patr*0228
e2259a1942 Jean*0229       ELSEIF (GM_taper_scheme.EQ.'fm07' ) THEN
                0230 C--   Ferrari & Mc.Williams, 2007:
                0231
5755dbe390 Patr*0232 #ifdef GM_EXCLUDE_FM07_TAP
                0233
                0234         STOP 'Need to compile without "#define GM_EXCLUDE_FM07_TAP"'
                0235
                0236 #else  /* GM_EXCLUDE_FM07_TAP */
                0237
e2259a1942 Jean*0238 C-    a) Calculate separately slope magnitude (SlopeMod)
                0239 C        and slope horiz. direction (Slope_X,Y : normalized)
14e0496834 Jean*0240         DO j=1-OLy+1,sNy+OLy-1
                0241          DO i=1-OLx+1,sNx+OLx-1
e2259a1942 Jean*0242           IF ( k.GT.kLow(i,j) ) THEN
                0243 C-        Bottom or below:
                0244             SlopeX  (i,j) = 0. _d 0
                0245             SlopeY  (i,j) = 0. _d 0
                0246             SlopeMod(i,j) = 0. _d 0
                0247             taperFct(i,j) = 0. _d 0
                0248           ELSE
                0249 C-        Above bottom:
5b172de0d2 Jean*0250            IF ( dSigmaDr(i,j).LE. GM_Small_Number )
                0251      &          dSigmaDr(i,j) = GM_Small_Number
ac0d9d59a8 Jean*0252            tmpFld(i,j) = dSigmaDx(i,j)*dSigmaDx(i,j)
                0253      &                 + dSigmaDy(i,j)*dSigmaDy(i,j)
5b172de0d2 Jean*0254            IF ( tmpFld(i,j).GT.zeroRL ) THEN
ac0d9d59a8 Jean*0255             locVar(i,j) = SQRT( tmpFld(i,j) )
5755dbe390 Patr*0256             SlopeX  (i,j) = dSigmaDx(i,j)/locVar(i,j)
                0257             SlopeY  (i,j) = dSigmaDy(i,j)/locVar(i,j)
5b172de0d2 Jean*0258             SlopeMod(i,j) = locVar(i,j)/dSigmaDr(i,j)
e2259a1942 Jean*0259             taperFct(i,j) = 1. _d 0
                0260            ELSE
                0261             SlopeX  (i,j) = 0. _d 0
                0262             SlopeY  (i,j) = 0. _d 0
                0263             SlopeMod(i,j) = 0. _d 0
                0264             taperFct(i,j) = 0. _d 0
                0265            ENDIF
                0266           ENDIF
                0267          ENDDO
                0268         ENDDO
                0269
                0270 C-    b) Set Transition Layer Depth:
ac0d9d59a8 Jean*0271         IF ( k.EQ.1 ) THEN
86e2f7bc1f Jean*0272           minTransLay = GM_facTrL2dz*( depthZ(k) - depthZ(k+1) )
ac0d9d59a8 Jean*0273         ELSE
86e2f7bc1f Jean*0274           minTransLay = GM_facTrL2dz*( depthZ(k-1) - depthZ(k) )
ac0d9d59a8 Jean*0275         ENDIF
14e0496834 Jean*0276         DO j=1-OLy+1,sNy+OLy-1
                0277          DO i=1-OLx+1,sNx+OLx-1
e2259a1942 Jean*0278           IF ( hTransLay(i,j).LE.0. _d 0 ) THEN
                0279 C-     previously below the transition layer
ac0d9d59a8 Jean*0280             tmpFld(i,j) = Lrho(i,j)*SlopeMod(i,j)
                0281 C-     ensure that transition layer is at least as thick than current level and
                0282 C      not thicker than the larger of the 2 : maxTransLay and facTrL2ML*hMixLay
e2259a1942 Jean*0283             dTransLay =
ac0d9d59a8 Jean*0284      &        MIN( MAX( tmpFld(i,j), minTransLay ),
86e2f7bc1f Jean*0285      &             MAX( GM_facTrL2ML*hMixLay(i,j), GM_maxTransLay ) )
e2259a1942 Jean*0286             IF ( k.GE.kLow(i,j) ) THEN
                0287 C-     bottom & below & 1rst level above the bottom :
                0288               recipLambda(i,j) = 0. _d 0
ac0d9d59a8 Jean*0289               baseSlope(i,j)   = SlopeMod(i,j)
e2259a1942 Jean*0290 C-- note: do not catch the 1rst level/interface (k=kLow) above the bottom
                0291 C         since no baseSlope has yet been stored (= 0); but might fit
                0292 C         well into transition layer criteria (if locally not stratified)
5b172de0d2 Jean*0293             ELSEIF ( dTransLay+hMixLay(i,j)+depthZ(k) .GE. zeroRL ) THEN
e2259a1942 Jean*0294 C-     Found the transition Layer : set depth of base of Transition layer
                0295               hTransLay(i,j) = -depthZ(k+1)
                0296 C      and compute inverse length scale "1/Lambda" from slope vert. grad
5b172de0d2 Jean*0297               IF ( baseSlope(i,j).GT.zeroRL ) THEN
ac0d9d59a8 Jean*0298                 recipLambda(i,j) = recipLambda(i,j)
                0299      &                           / MIN( baseSlope(i,j), GM_maxSlope )
e2259a1942 Jean*0300               ELSE
                0301                 recipLambda(i,j) = 0. _d 0
                0302               ENDIF
25ef785e26 Jean*0303 C      slope^2 & Kwz should remain > 0 : prevent too large negative D/lambda
5b172de0d2 Jean*0304               IF ( hMixLay(i,j)+depthZ(k+1).LT.zeroRL ) THEN
25ef785e26 Jean*0305                 rLambMin = 1. _d 0 /( hMixLay(i,j)+depthZ(k+1) )
                0306                 recipLambda(i,j) = MAX( recipLambda(i,j), rLambMin )
                0307               ENDIF
e2259a1942 Jean*0308             ELSE
                0309 C-     Still below Trans. layer: store slope & slope vert. grad.
ac0d9d59a8 Jean*0310               recipLambda(i,j) = ( MIN( SlopeMod(i,j), GM_maxSlope )
                0311      &                           - MIN( baseSlope(i,j), GM_maxSlope )
                0312      &                           ) / ( depthZ(k) - depthZ(k+1) )
e2259a1942 Jean*0313               baseSlope(i,j)   = SlopeMod(i,j)
                0314             ENDIF
                0315           ENDIF
                0316          ENDDO
                0317         ENDDO
                0318
                0319 C-    c) Set Slope component according to vertical position
                0320 C      (in Mixed-Layer / in Transition Layer / below Transition Layer)
14e0496834 Jean*0321         DO j=1-OLy+1,sNy+OLy-1
                0322          DO i=1-OLx+1,sNx+OLx-1
e2259a1942 Jean*0323           IF ( hTransLay(i,j).GT.0. _d 0 ) THEN
                0324 C-        already above base of transition layer:
                0325
ac0d9d59a8 Jean*0326             DoverLamb = (hTransLay(i,j)-hMixLay(i,j))*recipLambda(i,j)
e2259a1942 Jean*0327             IF ( -depthZ(k).LE.hMixLay(i,j) ) THEN
                0328 C-        compute tapering function -G(z) in the mixed Layer:
                0329               taperFctLoc =
                0330      &          ( -depthZ(k)/(hTransLay(i,j)+hMixLay(i,j))
ac0d9d59a8 Jean*0331      &            *( 2. _d 0 + DoverLamb )
e2259a1942 Jean*0332      &          )
                0333 C-        compute tapering function -G^(z) in the mixed Layer
                0334               taperFctHat =
                0335      &          ( -depthZ(k)/(hTransLay(i,j)+hMixLay(i,j))
                0336      &            *  2. _d 0
ac0d9d59a8 Jean*0337      &            *( 1. _d 0 + DoverLamb )
e2259a1942 Jean*0338      &          )
                0339             ELSE
                0340 C-        compute tapering function -G(z) in the transition Layer:
                0341               taperFctLoc =
                0342      &          ( -depthZ(k)/(hTransLay(i,j)+hMixLay(i,j))
ac0d9d59a8 Jean*0343      &            *( 2. _d 0 + DoverLamb )
e2259a1942 Jean*0344      &          )
                0345      &        - ( (depthZ(k)+hMixLay(i,j))*(depthZ(k)+hMixLay(i,j))
                0346      &            /( hTransLay(i,j)*hTransLay(i,j)
                0347      &               - hMixLay(i,j)*hMixLay(i,j)  )
                0348      &            *( 1. _d 0 + hTransLay(i,j)*recipLambda(i,j) )
                0349      &          )
                0350 C-        compute tapering function -G^(z) in the transition Layer:
                0351               taperFctHat =
                0352      &          ( -depthZ(k)/(hTransLay(i,j)+hMixLay(i,j))
                0353      &            *  2. _d 0
ac0d9d59a8 Jean*0354      &            *( 1. _d 0 + DoverLamb )
e2259a1942 Jean*0355      &          )
                0356      &        - ( (depthZ(k)+hMixLay(i,j))*(depthZ(k)+hMixLay(i,j))
                0357      &            /( hTransLay(i,j)*hTransLay(i,j)
                0358      &               - hMixLay(i,j)*hMixLay(i,j)  )
                0359      &            *( 1. _d 0 + hTransLay(i,j)*recipLambda(i,j)*2. _d 0 )
                0360      &          )
                0361             ENDIF
                0362 C-        modify the slope (above bottom of transition layer):
                0363 c           Smod = baseSlope(i,j)
                0364 C-        safer to limit the slope (even if it might never exceed GM_maxSlope)
                0365             Smod = MIN( baseSlope(i,j), GM_maxSlope )
                0366             SlopeX(i,j) = SlopeX(i,j)*Smod*taperFctLoc
                0367             SlopeY(i,j) = SlopeY(i,j)*Smod*taperFctLoc
                0368 c           SlopeSqr(i,j) = Smod*Smod*taperFctHat
df4facdcff Jean*0369 c           SlopeSqr(i,j) = baseSlope(i,j)*Smod*taperFctHat
                0370             SlopeSqr(i,j) = MIN( baseSlope(i,j), GM_bigSlope )
                0371      &                     *Smod*taperFctHat
e2259a1942 Jean*0372
                0373           ELSE
                0374 C--       Still below base of transition layer:
                0375 c           Smod = SlopeMod(i,j)
                0376 C-        safer to limit the slope:
                0377             Smod = MIN( SlopeMod(i,j), GM_maxSlope )
                0378             SlopeX(i,j) = SlopeX(i,j)*Smod
                0379             SlopeY(i,j) = SlopeY(i,j)*Smod
                0380 c           SlopeSqr(i,j) = Smod*Smod
df4facdcff Jean*0381 c           SlopeSqr(i,j) = SlopeMod(i,j)*Smod
                0382             SlopeSqr(i,j) = MIN( SlopeMod(i,j), GM_bigSlope )
                0383      &                     *Smod
e2259a1942 Jean*0384
                0385 C--       end if baseSlope > 0 / else => above/below base of Trans. Layer
                0386           ENDIF
                0387
                0388          ENDDO
                0389         ENDDO
                0390
5755dbe390 Patr*0391 #endif  /* GM_EXCLUDE_FM07_TAP */
                0392
8233d0ceb9 Jean*0393       ELSEIF (GM_taper_scheme.EQ.'ac02') THEN
2092dbb101 Patr*0394
bd3f833a36 Jean*0395 #ifdef GM_EXCLUDE_AC02_TAP
2092dbb101 Patr*0396
bd3f833a36 Jean*0397         STOP 'Need to compile without "#define GM_EXCLUDE_AC02_TAP"'
                0398
                0399 #else  /* GM_EXCLUDE_AC02_TAP */
2092dbb101 Patr*0400
bd3f833a36 Jean*0401 C-      New Scheme (A. & C. 2002): relax part of the small slope approximation
e2259a1942 Jean*0402 C         compute the true slope (no approximation)
bd3f833a36 Jean*0403 C         but still neglect Kxy & Kyx (assumed to be zero)
2092dbb101 Patr*0404
14e0496834 Jean*0405         DO j=1-OLy+1,sNy+OLy-1
                0406          DO i=1-OLx+1,sNx+OLx-1
5b172de0d2 Jean*0407           dRdSigmaLtd(i,j) = ( dSigmaDx(i,j)*dSigmaDx(i,j)
                0408      &                       + dSigmaDy(i,j)*dSigmaDy(i,j)
                0409      &                       )*convSlopeUnit*convSlopeUnit
                0410      &                       + dSigmaDr(i,j)*dSigmaDr(i,j)
20b8842d78 Patr*0411           taperFct(i,j) = 1. _d 0
e2259a1942 Jean*0412
5b172de0d2 Jean*0413           IF ( dRdSigmaLtd(i,j).NE.zeroRL ) THEN
8233d0ceb9 Jean*0414              dRdSigmaLtd(i,j) = 1. _d 0 / ( dRdSigmaLtd(i,j) )
                0415              SlopeSqr(i,j) = ( dSigmaDx(i,j)*dSigmaDx(i,j)
                0416      &                       + dSigmaDy(i,j)*dSigmaDy(i,j)
                0417      &                       )*dRdSigmaLtd(i,j)
5b172de0d2 Jean*0418              SlopeX(i,j) = dSigmaDx(i,j)
                0419      &                    *dRdSigmaLtd(i,j)*dSigmaDr(i,j)
                0420              SlopeY(i,j) = dSigmaDy(i,j)
                0421      &                    *dRdSigmaLtd(i,j)*dSigmaDr(i,j)
8233d0ceb9 Jean*0422           ELSE
                0423              SlopeSqr(i,j) = 0. _d 0
                0424              SlopeX(i,j) = 0. _d 0
                0425              SlopeY(i,j) = 0. _d 0
2092dbb101 Patr*0426           ENDIF
0f5add5564 Jean*0427 #ifndef ALLOW_AUTODIFF_TAMC
04522666de Ed H*0428 cph-- this part does not adjoint well
9cb619cfcd Patr*0429           IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND.
                0430      &         SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
                0431            taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
8233d0ceb9 Jean*0432           ELSEIF ( SlopeSqr(i,j) .GE. GM_slopeSqCutoff ) THEN
9cb619cfcd Patr*0433            taperFct(i,j) = 0. _d 0
                0434           ENDIF
                0435 #endif
2092dbb101 Patr*0436          ENDDO
                0437         ENDDO
                0438
bd3f833a36 Jean*0439 #endif /* GM_EXCLUDE_AC02_TAP */
2092dbb101 Patr*0440
e9fd580bc4 Jean*0441       ELSE
0c49347dc7 Alis*0442
bd3f833a36 Jean*0443 #ifdef GM_EXCLUDE_TAPERING
                0444
                0445         STOP 'Need to compile without "#define GM_EXCLUDE_TAPERING"'
                0446
                0447 #else  /* GM_EXCLUDE_TAPERING */
2092dbb101 Patr*0448
b6b11b9b2f Patr*0449 C----------------------------------------------------------------------
                0450
e673dad091 Patr*0451 C- Compute the slope, no clipping, but avoid reverse slope in negatively
5b172de0d2 Jean*0452 C                                  stratified (dSigmaDr < 0) region :
e673dad091 Patr*0453
                0454 #ifdef ALLOW_AUTODIFF_TAMC
113b21cd45 Mart*0455 C     Without this store directive, TAF generates an extra field anyway
                0456 C     so we do it here explicitly with a local tape and live without the
                0457 C     corresponding warning.
                0458 CADJ STORE dSigmaDr = loctape_gm
e673dad091 Patr*0459 #endif
                0460
14e0496834 Jean*0461         DO j=1-OLy+1,sNy+OLy-1
                0462          DO i=1-OLx+1,sNx+OLx-1
5b172de0d2 Jean*0463           IF ( dSigmaDr(i,j) .NE. zeroRL ) THEN
                0464            IF ( dSigmaDr(i,j) .LE. GM_Small_Number )
                0465      &         dSigmaDr(i,j) = GM_Small_Number
b6b11b9b2f Patr*0466           ENDIF
e673dad091 Patr*0467          ENDDO
                0468         ENDDO
                0469
14e0496834 Jean*0470         DO j=1-OLy+1,sNy+OLy-1
                0471          DO i=1-OLx+1,sNx+OLx-1
5b172de0d2 Jean*0472           IF ( dSigmaDr(i,j) .EQ. zeroRL ) THEN
                0473            IF ( dSigmaDx(i,j) .NE. zeroRL ) THEN
df4facdcff Jean*0474             SlopeX(i,j) = SIGN( GM_bigSlope, dSigmaDx(i,j) )
15dd2a0860 Patr*0475            ELSE
                0476             SlopeX(i,j) = 0. _d 0
                0477            ENDIF
5b172de0d2 Jean*0478            IF ( dSigmaDy(i,j) .NE. zeroRL ) THEN
df4facdcff Jean*0479             SlopeY(i,j) = SIGN( GM_bigSlope, dSigmaDy(i,j) )
15dd2a0860 Patr*0480            ELSE
                0481             SlopeY(i,j) = 0. _d 0
                0482            ENDIF
b6b11b9b2f Patr*0483           ELSE
549d1a8d8c Jean*0484            dRdSigmaLtd(i,j) = 1. _d 0 / dSigmaDr(i,j)
5b172de0d2 Jean*0485            SlopeX(i,j) = dSigmaDx(i,j)*dRdSigmaLtd(i,j)
                0486            SlopeY(i,j) = dSigmaDy(i,j)*dRdSigmaLtd(i,j)
                0487 c          SlopeX(i,j) = dSigmaDx(i,j)/dSigmaDr(i,j)
                0488 c          SlopeY(i,j) = dSigmaDy(i,j)/dSigmaDr(i,j)
b6b11b9b2f Patr*0489           ENDIF
e673dad091 Patr*0490          ENDDO
                0491         ENDDO
0c49347dc7 Alis*0492
e673dad091 Patr*0493 #ifdef ALLOW_AUTODIFF_TAMC
5a6997640b Mart*0494 cCADJ STORE slopeX(:,:) = comlev1_gmredi_slope, key=ikey, byte=isbyte
                0495 cCADJ STORE slopeY(:,:) = comlev1_gmredi_slope, key=ikey, byte=isbyte
                0496 CADJ STORE SlopeX = loctape_gm
                0497 CADJ STORE SlopeY = loctape_gm
e673dad091 Patr*0498 #endif
0c49347dc7 Alis*0499
14e0496834 Jean*0500         DO j=1-OLy+1,sNy+OLy-1
                0501          DO i=1-OLx+1,sNx+OLx-1
b6b11b9b2f Patr*0502           SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j)
                0503      &                   +SlopeY(i,j)*SlopeY(i,j)
                0504           taperFct(i,j) = 1. _d 0
5a6997640b Mart*0505 #ifdef ALLOW_AUTODIFF_TAMC
                0506          ENDDO
                0507         ENDDO
                0508 cCADJ STORE SlopeSqr = comlev1_gmredi_slope, key=ikey, byte=isbyte
                0509 CADJ STORE SlopeSqr = loctape_gm
                0510         DO j=1-OLy+1,sNy+OLy-1
                0511          DO i=1-OLx+1,sNx+OLx-1
                0512 #endif
8233d0ceb9 Jean*0513           IF ( SlopeSqr(i,j) .GE. GM_slopeSqCutoff ) THEN
                0514             slopeSqr(i,j) = GM_slopeSqCutoff
                0515             taperFct(i,j) = 0. _d 0
00e514b3a3 Jean*0516           ENDIF
0c49347dc7 Alis*0517          ENDDO
                0518         ENDDO
                0519
e9fd580bc4 Jean*0520 C- Compute the tapering function for the GM+Redi tensor :
0c49347dc7 Alis*0521
e9fd580bc4 Jean*0522        IF (GM_taper_scheme.EQ.'linear') THEN
                0523
                0524 C-      Simplest adiabatic tapering = Smax/Slope (linear)
14e0496834 Jean*0525         DO j=1-OLy+1,sNy+OLy-1
                0526          DO i=1-OLx+1,sNx+OLx-1
0c49347dc7 Alis*0527
5b172de0d2 Jean*0528           IF ( SlopeSqr(i,j) .EQ. zeroRL ) THEN
b6b11b9b2f Patr*0529            taperFct(i,j) = 1. _d 0
8233d0ceb9 Jean*0530           ELSEIF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND.
bd3f833a36 Jean*0531      &             SlopeSqr(i,j) .LT. GM_slopeSqCutoff )  THEN
e2259a1942 Jean*0532            taperFct(i,j) = SQRT(maxSlopeSqr / SlopeSqr(i,j))
df4facdcff Jean*0533            slopeSqr(i,j) = MIN( slopeSqr(i,j),GM_bigSlope*GM_bigSlope )
b6b11b9b2f Patr*0534           ENDIF
0c49347dc7 Alis*0535
e9fd580bc4 Jean*0536          ENDDO
                0537         ENDDO
                0538
                0539        ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
                0540
                0541 C-      Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
14e0496834 Jean*0542         DO j=1-OLy+1,sNy+OLy-1
                0543          DO i=1-OLx+1,sNx+OLx-1
e9fd580bc4 Jean*0544
5b172de0d2 Jean*0545           IF ( SlopeSqr(i,j) .EQ. zeroRL ) THEN
b6b11b9b2f Patr*0546            taperFct(i,j) = 1. _d 0
8233d0ceb9 Jean*0547           ELSEIF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND.
bd3f833a36 Jean*0548      &             SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
b6b11b9b2f Patr*0549            taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
                0550           ENDIF
0c49347dc7 Alis*0551
                0552          ENDDO
                0553         ENDDO
                0554
e9fd580bc4 Jean*0555        ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
0c49347dc7 Alis*0556
                0557 C-      Danabasoglu and McWilliams, J. Clim. 1995
14e0496834 Jean*0558         DO j=1-OLy+1,sNy+OLy-1
                0559          DO i=1-OLx+1,sNx+OLx-1
0c49347dc7 Alis*0560
5b172de0d2 Jean*0561           IF ( SlopeSqr(i,j) .EQ. zeroRL ) THEN
b6b11b9b2f Patr*0562            taperFct(i,j) = 1. _d 0
8233d0ceb9 Jean*0563           ELSEIF ( SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
5b172de0d2 Jean*0564            Smod = SQRT(SlopeSqr(i,j))*convSlopeUnit
8233d0ceb9 Jean*0565            taperFct(i,j) = op5*( oneRL + TANH( (GM_Scrit-Smod)/GM_Sd ) )
b6b11b9b2f Patr*0566           ENDIF
0c49347dc7 Alis*0567          ENDDO
                0568         ENDDO
                0569
e9fd580bc4 Jean*0570        ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
0c49347dc7 Alis*0571
                0572 C-      Large, Danabasoglu and Doney, JPO 1997
14e0496834 Jean*0573         DO j=1-OLy+1,sNy+OLy-1
                0574          DO i=1-OLx+1,sNx+OLx-1
0c49347dc7 Alis*0575
5b172de0d2 Jean*0576           IF ( SlopeSqr(i,j) .EQ. zeroRL ) THEN
b6b11b9b2f Patr*0577            taperFct(i,j) = 1. _d 0
07c79baaea Jean*0578           ELSEIF ( SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
8233d0ceb9 Jean*0579            Smod = SQRT(SlopeSqr(i,j))
5b172de0d2 Jean*0580            f1 = op5*( oneRL
                0581      &              + TANH( (GM_Scrit-Smod*convSlopeUnit)/GM_Sd ) )
                0582            Rnondim = rDepth/(Lrho(i,j)*Smod)
07c79baaea Jean*0583            IF ( Rnondim.GE.1. _d 0 ) THEN
                0584              f2 = 1. _d 0
                0585            ELSE
8233d0ceb9 Jean*0586              f2 = op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5) ) )
07c79baaea Jean*0587            ENDIF
8233d0ceb9 Jean*0588            taperFct(i,j) = f1*f2
b6b11b9b2f Patr*0589           ENDIF
0c49347dc7 Alis*0590
                0591          ENDDO
                0592         ENDDO
                0593
d8ee9652bd Gael*0594        ELSEIF (GM_taper_scheme.EQ.'stableGmAdjTap') THEN
                0595
0a637c269e Ou W*0596 #ifdef GMREDI_WITH_STABLE_ADJOINT
d8ee9652bd Gael*0597
0a637c269e Ou W*0598 C special choice for adjoint/optimization of parameters
                0599 C (~ strong clipping, reducing non linearity of kw=f(K))
d8ee9652bd Gael*0600
5a6997640b Mart*0601 C this variable could be relaced by the runtime parameter GM_maxSlope
1b7098b402 Gael*0602         slopeMax= 2. _d -3
d8ee9652bd Gael*0603
20dee61641 Mart*0604 # ifdef ALLOW_AUTODIFF_TAMC
5a6997640b Mart*0605 cCADJ STORE SlopeX(:,:) = comlev1_gmredi_slope, key=ikey, byte=isbyte
                0606 cCADJ STORE SlopeY(:,:) = comlev1_gmredi_slope, key=ikey, byte=isbyte
                0607 CADJ STORE SlopeX = loctape_gm
                0608 CADJ STORE SlopeY = loctape_gm
                0609 #endif
8233d0ceb9 Jean*0610         DO j=1-OLy+1,sNy+OLy-1
                0611          DO i=1-OLx+1,sNx+OLx-1
                0612           slopeSqTmp = SlopeX(i,j)*SlopeX(i,j)
                0613      &               + SlopeY(i,j)*SlopeY(i,j)
0a637c269e Ou W*0614
1b7098b402 Gael*0615           IF ( slopeSqTmp .GT. slopeMax**2 ) then
8233d0ceb9 Jean*0616            slopeTmp = SQRT(slopeSqTmp)
                0617            SlopeX(i,j) = SlopeX(i,j)*slopeMax/slopeTmp
                0618            SlopeY(i,j) = SlopeY(i,j)*slopeMax/slopeTmp
f5fe6ba1f3 Gael*0619           ENDIF
0a637c269e Ou W*0620          ENDDO
                0621         ENDDO
                0622
5a6997640b Mart*0623 C move the assignment of SlopeSqr to its own do-loop block from the do-loop
7c50f07931 Mart*0624 C block above to reduce TAF recomputations. The assignment of taperFct is
0a637c269e Ou W*0625 C also moved to keep the original order of operations.
8233d0ceb9 Jean*0626         DO j=1-OLy+1,sNy+OLy-1
                0627          DO i=1-OLx+1,sNx+OLx-1
d8ee9652bd Gael*0628           SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j)
8233d0ceb9 Jean*0629      &                  + SlopeY(i,j)*SlopeY(i,j)
d8ee9652bd Gael*0630           taperFct(i,j) = 1. _d 0
                0631          ENDDO
                0632         ENDDO
                0633
0a637c269e Ou W*0634 #else  /* GMREDI_WITH_STABLE_ADJOINT */
                0635
                0636         STOP 'Need to compile wth "#define GMREDI_WITH_STABLE_ADJOINT"'
                0637
d8ee9652bd Gael*0638 #endif /* GMREDI_WITH_STABLE_ADJOINT */
                0639
f42e64b3e7 Jean*0640        ELSEIF (GM_taper_scheme.NE.' ') THEN
                0641         STOP 'GMREDI_SLOPE_LIMIT: Bad GM_taper_scheme'
e9fd580bc4 Jean*0642        ENDIF
0c49347dc7 Alis*0643
bd3f833a36 Jean*0644 #endif /* GM_EXCLUDE_TAPERING */
2092dbb101 Patr*0645
e9fd580bc4 Jean*0646       ENDIF
0c49347dc7 Alis*0647
                0648 #endif /* ALLOW_GMREDI */
                0649
                0650       RETURN
                0651       END