Back to home page

MITgcm
 
 

    

File indexing completed on 2026-08-13 05:08:40 UTC

view on githubraw file Latest commit 09a9aa1d on 2026-08-10 19:46:52 UTC
89474f9a5c Mart*0001 #include "GGL90_OPTIONS.h"
7bb5a8a109 Jean*0002 #ifdef ALLOW_GENERIC_ADVDIFF
                0003 # include "GAD_OPTIONS.h"
                0004 #endif
dd9d13d532 Mart*0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
89474f9a5c Mart*0008 
                0009 CBOP
                0010 C !ROUTINE: GGL90_CALC
                0011 
                0012 C !INTERFACE: ======================================================
f688417df1 Jean*0013       SUBROUTINE GGL90_CALC(
dc6107c029 Jean*0014      I                 bi, bj, sigmaR, myTime, myIter, myThid )
                0015 
89474f9a5c Mart*0016 C !DESCRIPTION: \bv
5e48dccc42 Jean*0017 C     *==========================================================*
89474f9a5c Mart*0018 C     | SUBROUTINE GGL90_CALC                                    |
                0019 C     | o Compute all GGL90 fields defined in GGL90.h            |
5e48dccc42 Jean*0020 C     *==========================================================*
89474f9a5c Mart*0021 C     | Equation numbers refer to                                |
0320e25227 Mart*0022 C     |  Gaspar et al. (1990), JGR 95 (C9), pp 16,179            |
89474f9a5c Mart*0023 C     | Some parts of the implementation follow Blanke and       |
0320e25227 Mart*0024 C     |  Delecuse (1993), JPO, and OPA code, in particular the   |
                0025 C     |  computation of the                                      |
                0026 C     |  mixing length = max(min(lk,depth),lkmin)                |
                0027 C     | Note: Only call this S/R if Nr > 1 (no use if Nr=1)      |
5e48dccc42 Jean*0028 C     *==========================================================*
89474f9a5c Mart*0029 
                0030 C global parameters updated by ggl90_calc
5e48dccc42 Jean*0031 C     GGL90TKE     :: sub-grid turbulent kinetic energy          (m^2/s^2)
                0032 C     GGL90viscAz  :: GGL90 eddy viscosity coefficient             (m^2/s)
                0033 C     GGL90diffKzT :: GGL90 diffusion coefficient for temperature  (m^2/s)
89474f9a5c Mart*0034 C \ev
                0035 
                0036 C !USES: ============================================================
5e48dccc42 Jean*0037       IMPLICIT NONE
89474f9a5c Mart*0038 #include "SIZE.h"
                0039 #include "EEPARAMS.h"
                0040 #include "PARAMS.h"
                0041 #include "DYNVARS.h"
                0042 #include "FFIELDS.h"
                0043 #include "GRID.h"
f13fe90a48 Patr*0044 #include "GGL90.h"
f18a893d42 Mart*0045 #ifdef ALLOW_SHELFICE
                0046 # include "SHELFICE.h"
                0047 #endif
dd9d13d532 Mart*0048 #ifdef ALLOW_AUTODIFF_TAMC
                0049 # include "tamc.h"
                0050 #endif
7c50f07931 Mart*0051 
89474f9a5c Mart*0052 C !INPUT PARAMETERS: ===================================================
5e48dccc42 Jean*0053 C Routine arguments
dc6107c029 Jean*0054 C     bi, bj :: Current tile indices
                0055 C     sigmaR :: Vertical gradient of iso-neutral density
5e48dccc42 Jean*0056 C     myTime :: Current time in simulation
f688417df1 Jean*0057 C     myIter :: Current time-step number
5e48dccc42 Jean*0058 C     myThid :: My Thread Id number
89474f9a5c Mart*0059       INTEGER bi, bj
dc6107c029 Jean*0060       _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89474f9a5c Mart*0061       _RL     myTime
f688417df1 Jean*0062       INTEGER myIter
5e48dccc42 Jean*0063       INTEGER myThid
89474f9a5c Mart*0064 
                0065 #ifdef ALLOW_GGL90
                0066 
5b0716a6b3 Mart*0067 #ifdef ALLOW_DIAGNOSTICS
                0068       LOGICAL  DIAGNOSTICS_IS_ON
                0069       EXTERNAL DIAGNOSTICS_IS_ON
                0070 #endif /* ALLOW_DIAGNOSTICS */
                0071 
89474f9a5c Mart*0072 C !LOCAL VARIABLES: ====================================================
f688417df1 Jean*0073 C     iMin,iMax,jMin,jMax :: index boundaries of computation domain
0320e25227 Mart*0074 C     i, j, k          :: array computation indices
                0075 C     kSrf             :: vertical index of surface level
                0076 C     kTop             :: index of top interface (just below surf. level)
                0077 C     kBot             :: index of bottom interface (just above bottom lev.)
cdafb98dea Mart*0078 C     hFac/hFacI       :: fractional thickness of W-cell
f688417df1 Jean*0079 C     explDissFac      :: explicit Dissipation Factor (in [0-1])
                0080 C     implDissFac      :: implicit Dissipation Factor (in [0-1])
0320e25227 Mart*0081 C
                0082 C     In general, all 3D variables are defined at W-points (i.e.,
                0083 C     between k and k-1), all 2D variables are also defined at W-points
                0084 C     or at the very surface level (like uStarSquare)
                0085 C
f688417df1 Jean*0086 C     uStarSquare      :: square of friction velocity
                0087 C     verticalShear    :: (squared) vertical shear of horizontal velocity
                0088 C     Nsquare          :: squared buoyancy freqency
                0089 C     RiNumber         :: local Richardson number
                0090 C     KappaM           :: (local) viscosity parameter (eq.10)
                0091 C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
                0092 C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
                0093 C     TKEdissipation   :: dissipation of TKE
                0094 C     GGL90mixingLength:: mixing length of scheme following Banke+Delecuse
                0095 C         rMixingLength:: inverse of mixing length
                0096 C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
89474f9a5c Mart*0097       INTEGER iMin ,iMax ,jMin ,jMax
0320e25227 Mart*0098       INTEGER i, j, k
f18a893d42 Mart*0099       INTEGER kp1, km1
                0100       INTEGER kSrf, kTop, kBot
0320e25227 Mart*0101       INTEGER errCode
                0102       _RL deltaTloc
                0103       _RL explDissFac, implDissFac
                0104       _RL uStarSquare  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0105       _RL verticalShear(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0106       _RL KappaM       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0107       _RL KappaH
                0108 c     _RL Nsquare
                0109       _RL Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0110 c     _RL SQRTTKE
                0111       _RL SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0112       _RL RiNumber
87bca9545c Mart*0113 #ifdef ALLOW_GGL90_IDEMIX
0320e25227 Mart*0114       _RL IDEMIX_RiNumber
87bca9545c Mart*0115 #endif
0320e25227 Mart*0116       _RL TKEdissipation
b038e3cc4f Mart*0117       _RL tempU, tempUp, tempV, tempVp, prTemp, tmpVisc
0320e25227 Mart*0118       _RL TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0119       _RL GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0120       _RL rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0121       _RL KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0122       _RL GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
cdafb98dea Mart*0123 #ifdef ALLOW_GGL90_IDEMIX
0320e25227 Mart*0124       _RL hFacI            (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
5b0716a6b3 Mart*0125 C     IDEMIX_gTKE :: dissipation of internal wave energy is a source
                0126 C                    of TKE and mixing (output of S/R GGL90_IDEMIX)
f18a893d42 Mart*0127       _RL IDEMIX_gTKE      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
cdafb98dea Mart*0128 #endif /* ALLOW_GGL90_IDEMIX */
9293d3c672 Hajo*0129 #ifdef ALLOW_GGL90_LANGMUIR
                0130 C   uStar, vStar :: frictional velocity component
ee5f92f083 mjlo*0131       _RL depthFac
b038e3cc4f Mart*0132       _RL recip_Lasq, recip_LD
9293d3c672 Hajo*0133       _RL LCmixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0134       _RL stokesterm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0135       _RL dstokesUdR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0136       _RL dstokesVdR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
ee5f92f083 mjlo*0137       _RL uStar     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0138       _RL vStar     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
9293d3c672 Hajo*0139 #endif /* ALLOW_GGL90_LANGMUIR */
0320e25227 Mart*0140       _RL recip_hFacI      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0141       _RL hFac
f18a893d42 Mart*0142       _RS mskLoc
                0143 #ifdef ALLOW_GGL90_SMOOTH
                0144       _RS maskI            (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0145 #endif
b4ce400958 Davi*0146 C-    tri-diagonal matrix
0320e25227 Mart*0147       _RL a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0148       _RL b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0149       _RL c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0150 C     This mixed layer model is not invariant under coordinate
                0151 C     transformation to pressure coordinates, so we need these
                0152 C     factors to scale the vertical (pressure) coordinates
                0153       _RL coordFac, recip_coordFac
76f580e1f0 Mart*0154 #ifdef ALLOW_GGL90_HORIZDIFF
f688417df1 Jean*0155 C     xA, yA   :: area of lateral faces
                0156 C     dfx, dfy :: diffusive flux across lateral faces
                0157 C     gTKE     :: right hand side of diffusion equation
0320e25227 Mart*0158       _RL xA  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0159       _RL yA  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0160       _RL dfx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0161       _RL dfy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0162       _RL gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76f580e1f0 Mart*0163 #endif /* ALLOW_GGL90_HORIZDIFF */
004d5ee949 Davi*0164 #ifdef ALLOW_GGL90_SMOOTH
f6b150f7f1 Gael*0165       _RL p4, p8, p16
63bbd437b1 Jean*0166 #endif
f18a893d42 Mart*0167 #ifdef ALLOW_SHELFICE
                0168       INTEGER ki
                0169       _RL KE     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0170       _RL uFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0171       _RL vFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0172       _RL cDragU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0173       _RL cDragV (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0174       _RL stressU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0175       _RL stressV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0176       _RL kappaRX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0177 #endif
0320e25227 Mart*0178 #ifdef ALLOW_DIAGNOSTICS
5b0716a6b3 Mart*0179 # ifndef ALLOW_AUTODIFF
                0180       LOGICAL doDiagTKEmin
                0181       _RL recip_deltaT
                0182 # endif
0320e25227 Mart*0183       _RL surf_flx_tke(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0184 #endif /* ALLOW_DIAGNOSTICS */
dd9d13d532 Mart*0185 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0186 C     tkey :: tape key (depends on tiles)
                0187 C     kkey :: tape key (depends on levels and tiles)
                0188       INTEGER tkey, kkey
dd9d13d532 Mart*0189 #endif
31a3206180 Mart*0190 CEOP
63bbd437b1 Jean*0191 
                0192       PARAMETER( iMin = 2-OLx, iMax = sNx+OLx-1 )
                0193       PARAMETER( jMin = 2-OLy, jMax = sNy+OLy-1 )
                0194 #ifdef ALLOW_GGL90_SMOOTH
                0195       p4  = 0.25   _d 0
                0196       p8  = 0.125  _d 0
                0197       p16 = 0.0625 _d 0
004d5ee949 Davi*0198 #endif
89474f9a5c Mart*0199 
0320e25227 Mart*0200       IF ( usingPCoords ) THEN
                0201        kSrf = Nr
                0202        kTop = Nr
                0203       ELSE
                0204        kSrf =  1
                0205        kTop =  2
                0206       ENDIF
                0207       deltaTloc = dTtracerLev(kSrf)
                0208 
                0209       coordFac = 1. _d 0
                0210       IF ( usingPCoords) coordFac = gravity * rhoConst
                0211       recip_coordFac = 1./coordFac
                0212 
dd9d13d532 Mart*0213 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0214       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
dd9d13d532 Mart*0215 #endif /* ALLOW_AUTODIFF_TAMC */
                0216 
f688417df1 Jean*0217 C     explicit/implicit timestepping weights for dissipation
                0218       explDissFac = 0. _d 0
                0219       implDissFac = 1. _d 0 - explDissFac
89474f9a5c Mart*0220 
5b0716a6b3 Mart*0221 #ifdef ALLOW_DIAGNOSTICS
                0222 # ifndef ALLOW_AUTODIFF
                0223       doDiagTKEmin = .FALSE.
                0224 # endif
                0225       IF ( useDiagnostics ) THEN
                0226 # ifndef ALLOW_AUTODIFF
                0227        doDiagTKEmin = DIAGNOSTICS_IS_ON('GGL90Emn',myThid)
                0228 C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
                0229 C        does it only if k=1 (never the case here)
                0230        IF ( doDiagTKEmin )
                0231      &      CALL DIAGNOSTICS_COUNT('GGL90Emn',bi,bj,myThid)
                0232 # endif
                0233        DO j=1-OLy,sNy+OLy
                0234         DO i=1-OLx,sNx+OLx
                0235          surf_flx_tke(i,j) = 0.
                0236         ENDDO
                0237        ENDDO
                0238       ENDIF
                0239 #endif
                0240 
63bbd437b1 Jean*0241 C     For nonlinear free surface and especially with r*-coordinates, the
cdafb98dea Mart*0242 C     hFacs change every timestep, so we need to update them here in the
                0243 C     case of using IDEMIX.
7c50f07931 Mart*0244        DO k=1,Nr
0320e25227 Mart*0245         km1 = MAX(k-1,1)
cdafb98dea Mart*0246         DO j=1-OLy,sNy+OLy
                0247          DO i=1-OLx,sNx+OLx
63bbd437b1 Jean*0248           hFac =
5b0716a6b3 Mart*0249      &         MIN( halfRS, _hFacC(i,j,km1,bi,bj) )
                0250      &       + MIN( halfRS, _hFacC(i,j,k  ,bi,bj) )
7c50f07931 Mart*0251           recip_hFacI(i,j,k)=0. _d 0
cdafb98dea Mart*0252           IF ( hFac .NE. 0. _d 0 )
7c50f07931 Mart*0253      &         recip_hFacI(i,j,k)=1. _d 0/hFac
cdafb98dea Mart*0254 #ifdef ALLOW_GGL90_IDEMIX
                0255           hFacI(i,j,k) = hFac
                0256 #endif /* ALLOW_GGL90_IDEMIX */
                0257          ENDDO
                0258         ENDDO
                0259        ENDDO
                0260 
31f96e9372 Jean*0261 #ifdef ALLOW_GGL90_IDEMIX
                0262 C     step forward IDEMIX_E(energy) and compute tendency for TKE,
                0263 C     IDEMIX_gTKE = tau_d * IDEMIX_E**2, following Olbers and Eden (2013)
                0264       IF ( useIDEMIX) CALL GGL90_IDEMIX(
                0265      I     bi, bj, hFacI, recip_hFacI, sigmaR,
                0266      O     IDEMIX_gTKE,
                0267      I     myTime, myIter, myThid )
                0268 #endif /* ALLOW_GGL90_IDEMIX */
                0269 
89474f9a5c Mart*0270 C     Initialize local fields
0320e25227 Mart*0271       DO k=1,Nr
dc6107c029 Jean*0272        DO j=1-OLy,sNy+OLy
                0273         DO i=1-OLx,sNx+OLx
87bca9545c Mart*0274          rMixingLength(i,j,k)     = 0. _d 0
                0275          GGL90visctmp(i,j,k)      = 0. _d 0
f688417df1 Jean*0276          KappaE(i,j,k)            = 0. _d 0
                0277          TKEPrandtlNumber(i,j,k)  = 1. _d 0
                0278          GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
909cdb2275 Jean*0279 #ifndef SOLVE_DIAGONAL_LOWMEMORY
                0280          a3d(i,j,k) = 0. _d 0
                0281          b3d(i,j,k) = 1. _d 0
                0282          c3d(i,j,k) = 0. _d 0
                0283 #endif
87bca9545c Mart*0284          Nsquare(i,j,k) = 0. _d 0
                0285          SQRTTKE(i,j,k) = 0. _d 0
89474f9a5c Mart*0286         ENDDO
94c8eb5701 Jean*0287        ENDDO
89474f9a5c Mart*0288       ENDDO
dd9d13d532 Mart*0289 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0290 CADJ STORE GGL90TKE(:,:,:,bi,bj)=comlev1_bibj, key=tkey, kind=isbyte
dd9d13d532 Mart*0291 #endif
dc6107c029 Jean*0292       DO j=1-OLy,sNy+OLy
                0293        DO i=1-OLx,sNx+OLx
cdafb98dea Mart*0294         KappaM(i,j)        = 0. _d 0
0320e25227 Mart*0295         uStarSquare(i,j)   = 0. _d 0
cdafb98dea Mart*0296         verticalShear(i,j) = 0. _d 0
0320e25227 Mart*0297 c       rMixingLength(i,j,1)  = 0. _d 0
417b5b7e19 Oliv*0298 #ifdef ALLOW_AUTODIFF
b038e3cc4f Mart*0299         IF ( usingZCoords .AND. maskC(i,j,1,bi,bj).EQ.oneRS
                0300      &       .AND. GGL90TKE(i,j,1,bi,bj) .GT. zeroRL ) THEN
417b5b7e19 Oliv*0301 #endif
3e0545e2b7 Oliv*0302          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
417b5b7e19 Oliv*0303 #ifdef ALLOW_AUTODIFF
3e0545e2b7 Oliv*0304         ELSE
                0305          SQRTTKE(i,j,1) = 0. _d 0
                0306         ENDIF
417b5b7e19 Oliv*0307 #endif
87bca9545c Mart*0308 #ifdef ALLOW_GGL90_HORIZDIFF
                0309         xA(i,j)  = 0. _d 0
                0310         yA(i,j)  = 0. _d 0
                0311         dfx(i,j) = 0. _d 0
                0312         dfy(i,j) = 0. _d 0
                0313         gTKE(i,j) = 0. _d 0
                0314 #endif /* ALLOW_GGL90_HORIZDIFF */
89474f9a5c Mart*0315        ENDDO
                0316       ENDDO
                0317 
9293d3c672 Hajo*0318 #ifdef ALLOW_GGL90_LANGMUIR
                0319       IF (useLANGMUIR) THEN
                0320        recip_Lasq = 1. _d 0 / LC_num
                0321        recip_Lasq = recip_Lasq * recip_Lasq
                0322        recip_LD   = 4. _d 0 * PI / LC_lambda
                0323        DO j=1-OLy,sNy+OLy
                0324         DO i=1-OLx,sNx+OLx
                0325          stokesterm(i,j) = 0. _d 0
                0326          dstokesUdR(i,j) = 0. _d 0
                0327          dstokesVdR(i,j) = 0. _d 0
ee5f92f083 mjlo*0328          uStar(i,j) = SIGN( SQRT(ABS(surfaceForcingU(i,j,bi,bj))),
                0329      &                               surfaceForcingU(i,j,bi,bj)   )
                0330          vStar(i,j) = SIGN( SQRT(ABS(surfaceForcingV(i,j,bi,bj))),
                0331      &                               surfaceForcingV(i,j,bi,bj)   )
9293d3c672 Hajo*0332         ENDDO
                0333        ENDDO
                0334       ENDIF
                0335 #endif
                0336 
f688417df1 Jean*0337       DO k = 2, Nr
                0338        DO j=jMin,jMax
                0339         DO i=iMin,iMax
f18a893d42 Mart*0340          mskLoc = maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
417b5b7e19 Oliv*0341 #ifdef ALLOW_AUTODIFF
b038e3cc4f Mart*0342          IF ( mskLoc.EQ.oneRS
                0343      &        .AND. GGL90TKE(i,j,k,bi,bj) .GT. zeroRL ) THEN
417b5b7e19 Oliv*0344 #endif
3e0545e2b7 Oliv*0345           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
417b5b7e19 Oliv*0346 #ifdef ALLOW_AUTODIFF
3e0545e2b7 Oliv*0347          ELSE
                0348           SQRTTKE(i,j,k)=0. _d 0
                0349          ENDIF
417b5b7e19 Oliv*0350 #endif
f5bf4b6e4b Jean*0351 
89474f9a5c Mart*0352 C     buoyancy frequency
dc6107c029 Jean*0353          Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
0320e25227 Mart*0354      &                  * sigmaR(i,j,k) * coordFac
cdafb98dea Mart*0355 C     vertical shear term (dU/dz)^2+(dV/dz)^2 is computed later
                0356 C     to save some memory
b038e3cc4f Mart*0357 C     Initialise mixing length (eq. 2.35)
f688417df1 Jean*0358          GGL90mixingLength(i,j,k) = SQRTTWO *
004d5ee949 Davi*0359      &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
f18a893d42 Mart*0360      &        * mskLoc
004d5ee949 Davi*0361         ENDDO
                0362        ENDDO
                0363       ENDDO
                0364 
b038e3cc4f Mart*0365       CALL GGL90_MIXINGLENGTH(
                0366      U     GGL90mixingLength,
                0367 #ifdef ALLOW_GGL90_LANGMUIR
                0368      O     LCmixingLength,
dd9d13d532 Mart*0369 #endif
b038e3cc4f Mart*0370      O     rMixingLength,
                0371      I     iMin ,iMax ,jMin ,jMax,
                0372      I     bi, bj, myTime, myIter, myThid )
0320e25227 Mart*0373 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0374 CADJ STORE GGL90mixingLength = comlev1_bibj, key=tkey, kind=isbyte
b038e3cc4f Mart*0375 CADJ STORE rMixingLength     = comlev1_bibj, key=tkey, kind=isbyte
                0376 # ifdef ALLOW_GGL90_LANGMUIR
                0377 CADJ STORE LCmixingLength    = comlev1_bibj, key=tkey, kind=isbyte
                0378 # endif
0320e25227 Mart*0379 #endif
9293d3c672 Hajo*0380 
b038e3cc4f Mart*0381 C     start "proper" k-loop
004d5ee949 Davi*0382       DO k=2,Nr
f688417df1 Jean*0383        km1 = k-1
b038e3cc4f Mart*0384 #ifdef ALLOW_AUTODIFF_TAMC
                0385        kkey = k + (tkey-1)*Nr
                0386 #endif
f688417df1 Jean*0387 #ifdef ALLOW_GGL90_HORIZDIFF
0320e25227 Mart*0388        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
94c8eb5701 Jean*0389 C     horizontal diffusion of TKE (requires an exchange in
                0390 C      do_fields_blocking_exchanges)
76f580e1f0 Mart*0391 C     common factors
dc6107c029 Jean*0392         DO j=1-OLy,sNy+OLy
                0393          DO i=1-OLx,sNx+OLx
198cdce361 Mart*0394           xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
0320e25227 Mart*0395      &                 (MIN(.5 _d 0,_hFacW(i,j,km1,bi,bj) ) +
                0396      &                  MIN(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
198cdce361 Mart*0397           yA(i,j) = _dxG(i,j,bi,bj)*drC(k)*
0320e25227 Mart*0398      &                 (MIN(.5 _d 0,_hFacS(i,j,km1,bi,bj) ) +
                0399      &                  MIN(.5 _d 0,_hFacS(i,j,k  ,bi,bj) ) )
76f580e1f0 Mart*0400          ENDDO
94c8eb5701 Jean*0401         ENDDO
76f580e1f0 Mart*0402 C     Compute diffusive fluxes
                0403 C     ... across x-faces
dc6107c029 Jean*0404         DO j=1-OLy,sNy+OLy
                0405          dfx(1-OLx,j)=0. _d 0
                0406          DO i=1-OLx+1,sNx+OLx
76f580e1f0 Mart*0407           dfx(i,j) = -GGL90diffTKEh*xA(i,j)
                0408      &      *_recip_dxC(i,j,bi,bj)
                0409      &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
198cdce361 Mart*0410 #ifdef ISOTROPIC_COS_SCALING
76f580e1f0 Mart*0411      &      *CosFacU(j,bi,bj)
198cdce361 Mart*0412 #endif /* ISOTROPIC_COS_SCALING */
76f580e1f0 Mart*0413          ENDDO
                0414         ENDDO
                0415 C     ... across y-faces
25c8af7c05 Jean*0416         DO i=1-OLx,sNx+OLx
dc6107c029 Jean*0417          dfy(i,1-OLy)=0. _d 0
76f580e1f0 Mart*0418         ENDDO
dc6107c029 Jean*0419         DO j=1-OLy+1,sNy+OLy
                0420          DO i=1-OLx,sNx+OLx
76f580e1f0 Mart*0421           dfy(i,j) = -GGL90diffTKEh*yA(i,j)
                0422      &      *_recip_dyC(i,j,bi,bj)
                0423      &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
                0424 #ifdef ISOTROPIC_COS_SCALING
                0425      &      *CosFacV(j,bi,bj)
                0426 #endif /* ISOTROPIC_COS_SCALING */
                0427          ENDDO
94c8eb5701 Jean*0428         ENDDO
76f580e1f0 Mart*0429 C     Compute divergence of fluxes
dc6107c029 Jean*0430         DO j=1-OLy,sNy+OLy-1
                0431          DO i=1-OLx,sNx+OLx-1
31a3206180 Mart*0432           gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
cdafb98dea Mart*0433      &         *recip_hFacI(i,j,k)
198cdce361 Mart*0434      &         *((dfx(i+1,j)-dfx(i,j))
cdafb98dea Mart*0435      &         + (dfy(i,j+1)-dfy(i,j)) )
94c8eb5701 Jean*0436          ENDDO
76f580e1f0 Mart*0437         ENDDO
cdafb98dea Mart*0438 C     end if GGL90diffTKEh .eq. 0.
f688417df1 Jean*0439        ENDIF
                0440 #endif /* ALLOW_GGL90_HORIZDIFF */
                0441 
cdafb98dea Mart*0442 C     viscosity and diffusivity
9293d3c672 Hajo*0443 #ifdef ALLOW_GGL90_LANGMUIR
                0444        IF (useLANGMUIR) THEN
                0445         DO j=jMin,jMax
                0446          DO i=iMin,iMax
                0447           KappaM(i,j) = GGL90ck*LCmixingLength(i,j,k)*SQRTTKE(i,j,k)
                0448          ENDDO
                0449         ENDDO
                0450        ELSE
                0451 #endif
                0452         DO j=jMin,jMax
                0453          DO i=iMin,iMax
                0454           KappaM(i,j) = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
                0455 #ifdef ALLOW_GGL90_LANGMUIR
                0456          ENDDO
                0457         ENDDO
                0458        ENDIF
f688417df1 Jean*0459        DO j=jMin,jMax
                0460         DO i=iMin,iMax
9293d3c672 Hajo*0461 #endif
0320e25227 Mart*0462          GGL90visctmp(i,j,k) = MAX( KappaM(i,j),diffKrNrS(k)
                0463      &                            * recip_coordFac*recip_coordFac )
                0464      &        * maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
305c472a49 Jean*0465 C        note: storing GGL90visctmp like this, and using it later to compute
                0466 C              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
0320e25227 Mart*0467          KappaM(i,j) = MAX( KappaM(i,j),viscArNr(k)
                0468      &                    * recip_coordFac*recip_coordFac )
                0469      &                    * maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
cdafb98dea Mart*0470         ENDDO
                0471        ENDDO
31a3206180 Mart*0472 
63bbd437b1 Jean*0473 C     compute vertical shear (dU/dz)^2+(dV/dz)^2
                0474        IF ( calcMeanVertShear ) THEN
                0475 C     by averaging (@ grid-cell center) the 4 vertical shear compon @ U,V pos.
                0476         DO j=jMin,jMax
                0477          DO i=iMin,iMax
                0478           tempU  = ( uVel( i ,j,km1,bi,bj) - uVel( i ,j,k,bi,bj) )
                0479           tempUp = ( uVel(i+1,j,km1,bi,bj) - uVel(i+1,j,k,bi,bj) )
                0480           tempV  = ( vVel(i, j ,km1,bi,bj) - vVel(i, j ,k,bi,bj) )
                0481           tempVp = ( vVel(i,j+1,km1,bi,bj) - vVel(i,j+1,k,bi,bj) )
                0482           verticalShear(i,j) = (
31f96e9372 Jean*0483      &                 ( tempU*tempU + tempUp*tempUp )
                0484      &               + ( tempV*tempV + tempVp*tempVp )
                0485      &                         )*halfRL*recip_drC(k)*recip_drC(k)
0320e25227 Mart*0486      &                          *coordFac*coordFac
63bbd437b1 Jean*0487          ENDDO
                0488         ENDDO
                0489        ELSE
                0490 C     from the averaged flow at grid-cell center (2 compon x 2 pos.)
                0491         DO j=jMin,jMax
                0492          DO i=iMin,iMax
                0493           tempU = ( ( uVel(i,j,km1,bi,bj) + uVel(i+1,j,km1,bi,bj) )
                0494      &             -( uVel(i,j,k  ,bi,bj) + uVel(i+1,j,k  ,bi,bj) )
                0495      &            )*halfRL*recip_drC(k)
0320e25227 Mart*0496      &             *coordFac
63bbd437b1 Jean*0497           tempV = ( ( vVel(i,j,km1,bi,bj) + vVel(i,j+1,km1,bi,bj) )
                0498      &             -( vVel(i,j,k  ,bi,bj) + vVel(i,j+1,k  ,bi,bj) )
                0499      &            )*halfRL*recip_drC(k)
0320e25227 Mart*0500      &             *coordFac
63bbd437b1 Jean*0501           verticalShear(i,j) = tempU*tempU + tempV*tempV
                0502          ENDDO
                0503         ENDDO
                0504        ENDIF
b038e3cc4f Mart*0505 #ifdef ALLOW_AUTODIFF_TAMC
                0506 CADJ STORE kappaM        = comlev1_bibj_k, key = kkey, kind=isbyte
                0507 CADJ STORE verticalShear = comlev1_bibj_k, key = kkey, kind=isbyte
                0508 #endif /* ALLOW_AUTODIFF_TAMC */
63bbd437b1 Jean*0509 
9293d3c672 Hajo*0510 #ifdef ALLOW_GGL90_LANGMUIR
                0511        IF (useLANGMUIR) THEN
                0512 C       compute (dStokesU/dz) and (dStokesV/dz)
                0513         depthFac = recip_Lasq*EXP( recip_LD*rF(k) )
                0514         DO j=1-OLy,sNy+OLy
                0515          DO i=1-OLx,sNx+OLx
ee5f92f083 mjlo*0516           dstokesUdR(i,j) = recip_LD * uStar(i,j) * depthFac
                0517           dstokesVdR(i,j) = recip_LD * vStar(i,j) * depthFac
9293d3c672 Hajo*0518          ENDDO
                0519         ENDDO
                0520 
                0521         IF ( calcMeanVertShear ) THEN
                0522 C     by averaging (@ grid-cell center) the 4 vertical shear compon @
                0523 C     U,V pos.
                0524          DO j=jMin,jMax
                0525           DO i=iMin,iMax
                0526            tempU  = ( uVel( i ,j,km1,bi,bj) - uVel( i ,j,k,bi,bj) )
                0527            tempUp = ( uVel(i+1,j,km1,bi,bj) - uVel(i+1,j,k,bi,bj) )
                0528            tempV  = ( vVel(i, j ,km1,bi,bj) - vVel(i, j ,k,bi,bj) )
                0529            tempVp = ( vVel(i,j+1,km1,bi,bj) - vVel(i,j+1,k,bi,bj) )
                0530            stokesterm(i,j) = (
                0531      &              ( tempU *dstokesUdR(i,j)
                0532      &               +tempUp*dstokesUdR(i+1,j) )
                0533      &             +( tempV *dstokesVdR(i,j)
                0534      &               +tempVp*dstokesVdR(i,j+1) )
                0535      &                       )*halfRL*recip_drC(k)*coordFac*coordFac
                0536           ENDDO
                0537          ENDDO
                0538         ELSE
                0539 C     from the averaged flow at grid-cell center (2 compon x 2 pos.)
                0540          DO j=jMin,jMax
                0541           DO i=iMin,iMax
                0542            tempU = ( ( uVel(i,j,km1,bi,bj) + uVel(i+1,j,km1,bi,bj) )
                0543      &              -( uVel(i,j,k  ,bi,bj) + uVel(i+1,j,k  ,bi,bj) )
                0544      &             )*halfRL*recip_drC(k)
                0545      &              *coordFac
                0546            tempV = ( ( vVel(i,j,km1,bi,bj) + vVel(i,j+1,km1,bi,bj) )
                0547      &              -( vVel(i,j,k  ,bi,bj) + vVel(i,j+1,k  ,bi,bj) )
                0548      &             )*halfRL*recip_drC(k)
                0549      &              *coordFac
                0550            stokesterm(i,j) = halfRL*coordFac*(
                0551      &                       tempU*(dstokesUdR(i,j)+dstokesUdR(i+1,j))
                0552      &                     + tempV*(dstokesVdR(i,j)+dstokesVdR(i,j+1))
                0553      &                       )
                0554           ENDDO
                0555          ENDDO
                0556         ENDIF
                0557        ENDIF
                0558 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0559 CADJ STORE stokesterm = comlev1_bibj, key=tkey, kind=isbyte
9293d3c672 Hajo*0560 # endif
                0561 #endif /* ALLOW_GGL90_LANGMUIR */
                0562 
5b0716a6b3 Mart*0563 C     compute Prandtl number (always greater than 1)
31a3206180 Mart*0564 #ifdef ALLOW_GGL90_IDEMIX
cdafb98dea Mart*0565        IF ( useIDEMIX ) THEN
63bbd437b1 Jean*0566         DO j=jMin,jMax
                0567          DO i=iMin,iMax
                0568           RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
                0569      &         /(verticalShear(i,j)+GGL90eps)
                0570           IDEMIX_RiNumber = MAX( KappaM(i,j)*Nsquare(i,j,k), 0. _d 0)/
31f96e9372 Jean*0571      &         ( GGL90eps + IDEMIX_gTKE(i,j,k) )
5b0716a6b3 Mart*0572           prTemp          = 6.6 _d 0 * MIN( RiNumber, IDEMIX_RiNumber )
63bbd437b1 Jean*0573           TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
0320e25227 Mart*0574           TKEPrandtlNumber(i,j,k) = MAX( oneRL,TKEPrandtlNumber(i,j,k) )
63bbd437b1 Jean*0575          ENDDO
cdafb98dea Mart*0576         ENDDO
                0577        ELSE
                0578 #endif /* ALLOW_GGL90_IDEMIX */
63bbd437b1 Jean*0579         DO j=jMin,jMax
                0580          DO i=iMin,iMax
                0581           RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
                0582      &         /(verticalShear(i,j)+GGL90eps)
                0583           prTemp = 1. _d 0
                0584           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
                0585           TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
                0586          ENDDO
cdafb98dea Mart*0587         ENDDO
f2a88c9ff8 jm-c 0588 #ifdef ALLOW_GGL90_IDEMIX
cdafb98dea Mart*0589        ENDIF
f2a88c9ff8 jm-c 0590 #endif /* ALLOW_GGL90_IDEMIX */
b038e3cc4f Mart*0591 #ifdef ALLOW_AUTODIFF_TAMC
                0592 CADJ STORE TKEPrandtlNumber(:,:,k)=comlev1_bibj_k,key=kkey,kind=isbyte
                0593 #endif /* ALLOW_AUTODIFF_TAMC */
31a3206180 Mart*0594 
cdafb98dea Mart*0595        DO j=jMin,jMax
                0596         DO i=iMin,iMax
305c472a49 Jean*0597 C        diffusivity
cdafb98dea Mart*0598          KappaH = KappaM(i,j)/TKEPrandtlNumber(i,j,k)
0320e25227 Mart*0599          KappaE(i,j,k) = GGL90alpha * KappaM(i,j)
                0600      &        * maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
f688417df1 Jean*0601 
                0602 C     dissipation term
                0603          TKEdissipation = explDissFac*GGL90ceps
                0604      &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
                0605      &        *GGL90TKE(i,j,k,bi,bj)
                0606 C     partial update with sum of explicit contributions
                0607          GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
0320e25227 Mart*0608      &        + deltaTloc*(
cdafb98dea Mart*0609      &        + KappaM(i,j)*verticalShear(i,j)
f688417df1 Jean*0610      &        - KappaH*Nsquare(i,j,k)
                0611      &        - TKEdissipation
                0612      &        )
                0613         ENDDO
76f580e1f0 Mart*0614        ENDDO
f688417df1 Jean*0615 
cdafb98dea Mart*0616 #ifdef ALLOW_GGL90_IDEMIX
                0617        IF ( useIDEMIX ) THEN
                0618 C     add IDEMIX contribution to the turbulent kinetic energy
                0619         DO j=jMin,jMax
                0620          DO i=iMin,iMax
                0621           GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
31f96e9372 Jean*0622      &         + deltaTloc*IDEMIX_gTKE(i,j,k)
cdafb98dea Mart*0623          ENDDO
                0624         ENDDO
                0625        ENDIF
                0626 #endif /* ALLOW_GGL90_IDEMIX */
                0627 
9293d3c672 Hajo*0628 #ifdef ALLOW_GGL90_LANGMUIR
                0629        IF ( useLANGMUIR ) THEN
9af873c532 Hajo*0630 C     add Langmuir contribution to the turbulent kinetic energy
9293d3c672 Hajo*0631         DO j=jMin,jMax
                0632          DO i=iMin,iMax
                0633           GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
                0634      &         + deltaTloc*(KappaM(i,j)*stokesterm(i,j))
                0635          ENDDO
                0636         ENDDO
                0637        ENDIF
                0638 #endif /* ALLOW_GGL90_LANGMUIR */
                0639 
f688417df1 Jean*0640 #ifdef ALLOW_GGL90_HORIZDIFF
                0641        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
                0642 C--    Add horiz. diffusion tendency
                0643         DO j=jMin,jMax
                0644          DO i=iMin,iMax
                0645           GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
0320e25227 Mart*0646      &                          + gTKE(i,j)*deltaTloc
f688417df1 Jean*0647          ENDDO
                0648         ENDDO
                0649        ENDIF
76f580e1f0 Mart*0650 #endif /* ALLOW_GGL90_HORIZDIFF */
                0651 
f688417df1 Jean*0652 C--   end of k loop
                0653       ENDDO
0320e25227 Mart*0654       IF ( usingPCoords ) THEN
                0655 C     impose TKE(1) = 0.
                0656        DO j=jMin,jMax
                0657         DO i=iMin,iMax
                0658          GGL90TKE(i,j,1,bi,bj) = 0. _d 0
                0659         ENDDO
                0660        ENDDO
                0661       ENDIF
f688417df1 Jean*0662 
0320e25227 Mart*0663 C     ==============================================
76f580e1f0 Mart*0664 C     Implicit time step to update TKE for k=1,Nr;
0320e25227 Mart*0665 C     TKE(Nr+1)=0 by default;
                0666 C     for pressure coordinates, this translates into
                0667 C     TKE(1)  = 0, TKE(Nr+1) is the surface value
                0668 C     ==============================================
89474f9a5c Mart*0669 C     set up matrix
76f580e1f0 Mart*0670 C--   Lower diagonal
89474f9a5c Mart*0671       DO j=jMin,jMax
                0672        DO i=iMin,iMax
909cdb2275 Jean*0673          a3d(i,j,1) = 0. _d 0
89474f9a5c Mart*0674        ENDDO
                0675       ENDDO
                0676       DO k=2,Nr
5b0716a6b3 Mart*0677 #ifdef GGL90_MISSING_HFAC_BUG
                0678        IF ( .NOT.useIDEMIX ) THEN
                0679         DO j=1-OLy,sNy+OLy
                0680          DO i=1-OLx,sNx+OLx
                0681           recip_hFacI(i,j,k) = oneRS
                0682          ENDDO
                0683         ENDDO
                0684        ENDIF
                0685 #endif
37a95f6b42 Davi*0686        km1=MAX(2,k-1)
89474f9a5c Mart*0687        DO j=jMin,jMax
                0688         DO i=iMin,iMax
0320e25227 Mart*0689          IF ( usingPCoords) km1=MIN(Nr,MAX(kSurfC(i,j,bi,bj)+1,k-1))
37a95f6b42 Davi*0690 C-    We keep recip_hFacC in the diffusive flux calculation,
                0691 C-    but no hFacC in TKE volume control
                0692 C-    No need for maskC(k-1) with recip_hFacC(k-1)
0320e25227 Mart*0693          a3d(i,j,k) = -deltaTloc
004d5ee949 Davi*0694      &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
73c2f90ab3 Davi*0695      &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
5b0716a6b3 Mart*0696      &        *recip_drC(k)*maskC(i,j,k,bi,bj)*recip_hFacI(i,j,k)
0320e25227 Mart*0697      &        *coordFac*coordFac
89474f9a5c Mart*0698         ENDDO
                0699        ENDDO
                0700       ENDDO
76f580e1f0 Mart*0701 C--   Upper diagonal
89474f9a5c Mart*0702       DO j=jMin,jMax
                0703        DO i=iMin,iMax
909cdb2275 Jean*0704          c3d(i,j,1)  = 0. _d 0
89474f9a5c Mart*0705        ENDDO
                0706       ENDDO
004d5ee949 Davi*0707       DO k=2,Nr
0320e25227 Mart*0708        kp1=MIN(k+1,Nr)
89474f9a5c Mart*0709        DO j=jMin,jMax
                0710         DO i=iMin,iMax
0320e25227 Mart*0711          IF ( usingZCoords ) kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
37a95f6b42 Davi*0712 C-    We keep recip_hFacC in the diffusive flux calculation,
                0713 C-    but no hFacC in TKE volume control
                0714 C-    No need for maskC(k) with recip_hFacC(k)
0320e25227 Mart*0715          c3d(i,j,k) = -deltaTloc
5b0716a6b3 Mart*0716      &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
                0717      &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
                0718      &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)*recip_hFacI(i,j,k)
                0719      &        *coordFac*coordFac
89474f9a5c Mart*0720         ENDDO
                0721        ENDDO
                0722       ENDDO
31a3206180 Mart*0723 
                0724       IF (.NOT.GGL90_dirichlet) THEN
                0725 C      Neumann bottom boundary condition for TKE: no flux from bottom
0320e25227 Mart*0726        IF ( usingPCoords ) THEN
                0727         DO j=jMin,jMax
                0728          DO i=iMin,iMax
                0729           kBot = MIN(kSurfC(i,j,bi,bj)+1,Nr)
                0730           a3d(i,j,kBot) = 0. _d 0
                0731          ENDDO
31a3206180 Mart*0732         ENDDO
0320e25227 Mart*0733        ELSE
                0734         DO j=jMin,jMax
                0735          DO i=iMin,iMax
                0736           kBot = MAX(kLowC(i,j,bi,bj),1)
                0737           c3d(i,j,kBot) = 0. _d 0
                0738          ENDDO
                0739         ENDDO
                0740        ENDIF
31a3206180 Mart*0741       ENDIF
                0742 
76f580e1f0 Mart*0743 C--   Center diagonal
89474f9a5c Mart*0744       DO k=1,Nr
37a95f6b42 Davi*0745        km1 = MAX(k-1,1)
89474f9a5c Mart*0746        DO j=jMin,jMax
                0747         DO i=iMin,iMax
cdafb98dea Mart*0748          b3d(i,j,k) = 1. _d 0 - c3d(i,j,k) - a3d(i,j,k)
0320e25227 Mart*0749      &        + implDissFac*deltaTloc*GGL90ceps*SQRTTKE(i,j,k)
f688417df1 Jean*0750      &        * rMixingLength(i,j,k)
37a95f6b42 Davi*0751      &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
cdafb98dea Mart*0752         ENDDO
89474f9a5c Mart*0753        ENDDO
                0754       ENDDO
0320e25227 Mart*0755       IF ( usingPCoords ) THEN
                0756 C     impose TKE(1) = 0.
                0757        DO j=jMin,jMax
                0758         DO i=iMin,iMax
                0759          b3d(i,j,1) = 1. _d 0
                0760         ENDDO
                0761        ENDDO
                0762       ENDIF
89474f9a5c Mart*0763 C     end set up matrix
                0764 
                0765 C     Apply boundary condition
31f96e9372 Jean*0766       IF ( calcMeanVertShear ) THEN
                0767 C     by averaging (@ grid-cell center) the 4 components @ U,V pos.
                0768        DO j=jMin,jMax
                0769         DO i=iMin,iMax
                0770          tempU  = surfaceForcingU( i ,j,bi,bj)
                0771          tempUp = surfaceForcingU(i+1,j,bi,bj)
                0772          tempV  = surfaceForcingV(i, j ,bi,bj)
                0773          tempVp = surfaceForcingV(i,j+1,bi,bj)
b038e3cc4f Mart*0774          uStarSquare(i,j) =
31f96e9372 Jean*0775      &        ( tempU*tempU + tempUp*tempUp
                0776      &        + tempV*tempV + tempVp*tempVp
b038e3cc4f Mart*0777      &        )*halfRL
31f96e9372 Jean*0778 C Note: adding parenthesis in 4 terms sum (-> 2 group of 2) as below:
b038e3cc4f Mart*0779 c        uStarSquare(i,j) =
31f96e9372 Jean*0780 c    &        ( ( tempU*tempU + tempUp*tempUp )
                0781 c    &        + ( tempV*tempV + tempVp*tempVp )
b038e3cc4f Mart*0782 c    &        )*halfRL
31f96e9372 Jean*0783 C       seems to break restart !
                0784         ENDDO
                0785        ENDDO
                0786       ELSE
                0787        DO j=jMin,jMax
                0788         DO i=iMin,iMax
89474f9a5c Mart*0789 C     estimate friction velocity uStar from surface forcing
b038e3cc4f Mart*0790          uStarSquare(i,j) =
31f96e9372 Jean*0791      &     ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
                0792      &               + surfaceForcingU(i+1,j,  bi,bj) ) )**2
                0793      &   + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
                0794      &               + surfaceForcingV(i,  j+1,bi,bj) ) )**2
                0795         ENDDO
89474f9a5c Mart*0796        ENDDO
31f96e9372 Jean*0797       ENDIF
f18a893d42 Mart*0798 #ifdef ALLOW_SHELFICE
                0799 C     uStarSquare should not have any effect undernath iceshelves, but
                0800 C     masking uStarSquare underneath ice shelf is not necessary because
                0801 C     surfaceForcingU/V=0 in this case (see shelfice_forcing_surf.F)
                0802 C     Instead, uStarSquare is computed from the sub-glacial drag.
                0803       IF ( useSHELFICE .AND.
                0804      &     ( no_slip_shelfice .OR. SHELFICEDragLinear.NE.zeroRL
                0805      &                        .OR. SHELFICEselectDragQuadr.GE.0 )
                0806      &   ) THEN
                0807 C     First, we need to compute an early estimate of the drag
                0808 C     coefficients based ond the ocean velocity of the previous time
                0809 C     step only; because kappaRU and kappaRV are not yet available,
                0810 C     kappyRX is just set to horizontally constant viscosity
                0811 C     coefficients.
                0812        DO j=1-OLy,sNy+OLy
                0813         DO i=1-OLx,sNx+OLx
                0814          stressU(i,j) = 0. _d 0
                0815          stressV(i,j) = 0. _d 0
                0816         ENDDO
                0817        ENDDO
                0818        DO k=1,Nr+1
                0819         ki = MIN(k,Nr)
                0820         DO j=1-OLy,sNy+OLy
                0821          DO i=1-OLx,sNx+OLx
                0822           kappaRX(i,j,k) = viscArNr(ki)
                0823          ENDDO
                0824         ENDDO
                0825        ENDDO
                0826        DO k=1,Nr
                0827         DO j=1-OLy,sNy+OLy
                0828          DO i=1-OLx,sNx+OLx
                0829           KE  (i,j) = 0. _d 0
                0830           uFld(i,j) = uVel(i,j,k,bi,bj)
                0831           vFld(i,j) = vVel(i,j,k,bi,bj)
                0832          ENDDO
                0833         ENDDO
                0834         CALL SHELFICE_U_DRAG_COEFF( bi, bj, k, .FALSE.,
                0835      I       uFld, vFld, kappaRX, KE,
                0836      O       cDragU,
                0837      I       myIter, myThid )
                0838         CALL SHELFICE_V_DRAG_COEFF( bi, bj, k, .FALSE.,
                0839      I       uFld, vFld, kappaRX, KE,
                0840      O       cDragV,
                0841      I       myIter, myThid )
                0842 C-     compute explicit stress
                0843         DO j=1-OLy,sNy+OLy
                0844          DO i=1-OLx,sNx+OLx
                0845           stressU(i,j) = stressU(i,j) - cDragU(i,j)*uFld(i,j)*rUnit2mass
                0846           stressV(i,j) = stressV(i,j) - cDragV(i,j)*vFld(i,j)*rUnit2mass
                0847          ENDDO
                0848         ENDDO
                0849        ENDDO
                0850 C
                0851        DO j=jMin,jMax
09a9aa1d4d Mich*0852         DO i=iMin,iMax
f18a893d42 Mart*0853          IF ( kTopC(i,j,bi,bj) .GT. 0 ) THEN
b038e3cc4f Mart*0854           uStarSquare(i,j) =
f18a893d42 Mart*0855      &         ( stressU(i,j)*stressU(i,j)+stressU(i+1,j)*stressU(i+1,j)
                0856      &         + stressV(i,j)*stressV(i,j)+stressV(i,j+1)*stressV(i,j+1)
b038e3cc4f Mart*0857      &         )*halfRL
f18a893d42 Mart*0858          ENDIF
                0859         ENDDO
                0860        ENDDO
                0861       ENDIF
                0862 #endif
b038e3cc4f Mart*0863 #ifdef ALLOW_AUTODIFF_TAMC
                0864 CADJ STORE uStarSquare      = comlev1_bibj, key = tkey, kind = isbyte
                0865 #endif
                0866       DO j=jMin,jMax
                0867        DO i=iMin,iMax
                0868 #ifdef ALLOW_AUTODIFF
                0869          IF ( uStarSquare(i,j) .GT. zeroRL )
                0870      &        uStarSquare(i,j) = SQRT(uStarSquare(i,j))*recip_coordFac
                0871 #else
                0872          uStarSquare(i,j) = SQRT(uStarSquare(i,j))*recip_coordFac
                0873 #endif
                0874        ENDDO
                0875       ENDDO
                0876 #ifdef ALLOW_AUTODIFF_TAMC
                0877 C     avoid recomputing the above multiple times in AD routine
                0878 CADJ STORE TKEPrandtlNumber = comlev1_bibj, key = tkey, kind = isbyte
                0879 CADJ STORE GGL90visctmp     = comlev1_bibj, key = tkey, kind = isbyte
                0880 CADJ STORE kappaE           = comlev1_bibj, key = tkey, kind = isbyte
                0881 CADJ STORE a3d, b3d, c3d    = comlev1_bibj, key = tkey, kind = isbyte
                0882 #endif
f18a893d42 Mart*0883 
0320e25227 Mart*0884 C     Dirichlet surface boundary condition for TKE
                0885       IF ( usingPCoords ) THEN
                0886        DO j=jMin,jMax
                0887         DO i=iMin,iMax
f18a893d42 Mart*0888 CML#ifdef ALLOW_SHELFICE
                0889 CML         IF ( useShelfIce ) THEN
                0890 CML          kSrf = MAX(kTopC(i,j,bi,bj),1)
                0891 CML          kTop = kSrf
                0892 CML         ENDIF
                0893 CML#endif
0320e25227 Mart*0894          GGL90TKE(i,j,kSrf,bi,bj) = GGL90TKE(i,j,kSrf,bi,bj)
                0895      &        - c3d(i,j,kSrf) * maskC(i,j,kSrf,bi,bj)
                0896      &        *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare(i,j))
                0897          c3d(i,j,kSrf) = 0. _d 0
                0898         ENDDO
                0899        ENDDO
                0900       ELSE
31a3206180 Mart*0901        DO j=jMin,jMax
                0902         DO i=iMin,iMax
f18a893d42 Mart*0903 #ifdef ALLOW_SHELFICE
                0904          IF ( useShelfIce ) THEN
                0905           kSrf = MAX(1,kTopC(i,j,bi,bj))
                0906           kTop = MIN(kSrf+1,Nr)
                0907          ENDIF
                0908 #endif
0320e25227 Mart*0909          GGL90TKE(i,j,kSrf,bi,bj) = maskC(i,j,kSrf,bi,bj)
                0910      &        *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare(i,j))
                0911          GGL90TKE(i,j,kTop,bi,bj) = GGL90TKE(i,j,kTop,bi,bj)
                0912      &        - a3d(i,j,kTop)*GGL90TKE(i,j,kSrf,bi,bj)
                0913          a3d(i,j,kTop) = 0. _d 0
31a3206180 Mart*0914         ENDDO
                0915        ENDDO
                0916       ENDIF
                0917 
0320e25227 Mart*0918       IF (GGL90_dirichlet) THEN
                0919 C      Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
                0920        IF ( usingPCoords ) THEN
                0921         DO j=jMin,jMax
                0922          DO i=iMin,iMax
                0923           kBot = MIN(kSurfC(i,j,bi,bj)+1,Nr)
                0924           GGL90TKE(i,j,kBot,bi,bj) = GGL90TKE(i,j,kBot,bi,bj)
                0925      &                             - GGL90TKEbottom*a3d(i,j,kBot)
                0926           a3d(i,j,kBot) = 0. _d 0
                0927          ENDDO
                0928         ENDDO
                0929        ELSE
                0930         DO j=jMin,jMax
                0931          DO i=iMin,iMax
                0932           kBot = MAX(kLowC(i,j,bi,bj),1)
                0933           GGL90TKE(i,j,kBot,bi,bj) = GGL90TKE(i,j,kBot,bi,bj)
                0934      &                             - GGL90TKEbottom*c3d(i,j,kBot)
                0935           c3d(i,j,kBot) = 0. _d 0
                0936          ENDDO
                0937         ENDDO
                0938        ENDIF
                0939       ENDIF
                0940 
dd9d13d532 Mart*0941 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0942 CADJ STORE GGL90TKE(:,:,:,bi,bj)=comlev1_bibj, key=tkey, kind=isbyte
dd9d13d532 Mart*0943 #endif
f688417df1 Jean*0944 C     solve tri-diagonal system
fa8d87e2db Jean*0945       errCode = -1
37a95f6b42 Davi*0946       CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
909cdb2275 Jean*0947      I                        a3d, b3d, c3d,
8a58850ca8 Jean*0948      U                        GGL90TKE(1-OLx,1-OLy,1,bi,bj),
37a95f6b42 Davi*0949      O                        errCode,
f688417df1 Jean*0950      I                        bi, bj, myThid )
f5bf4b6e4b Jean*0951 
dd9d13d532 Mart*0952 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0953 CADJ STORE GGL90TKE(:,:,:,bi,bj)=comlev1_bibj, key=tkey, kind=isbyte
dd9d13d532 Mart*0954 #endif
0320e25227 Mart*0955       DO k=2,Nr
5b0716a6b3 Mart*0956 #if ( defined ALLOW_DIAGNOSTICS && !defined ALLOW_AUTODIFF )
                0957 C     This diagnostics code causes extra recomputations so we skip it.
                0958        IF ( doDiagTKEmin ) THEN
                0959         DO j=1,sNy
                0960          DO i=1,sNx
                0961           surf_flx_tke(i,j) = GGL90TKE(i,j,k,bi,bj)
                0962      &         * maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
                0963          ENDDO
                0964         ENDDO
                0965        ENDIF
                0966 #endif
f688417df1 Jean*0967        DO j=jMin,jMax
                0968         DO i=iMin,iMax
0320e25227 Mart*0969 C     impose minimum TKE to avoid numerical undershoots below zero;
                0970 C     level k=1 is either prescribed surface boundary condition (z-coords) or
                0971 C     bottom boundary conditions, which by definition is zero
                0972          GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
f688417df1 Jean*0973      &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
89474f9a5c Mart*0974         ENDDO
                0975        ENDDO
5b0716a6b3 Mart*0976 #if ( defined ALLOW_DIAGNOSTICS && !defined ALLOW_AUTODIFF )
                0977        IF ( doDiagTKEmin ) THEN
                0978         recip_deltaT = 1. _d 0 / deltaTloc
                0979         DO j=1,sNy
                0980          DO i=1,sNx
                0981           surf_flx_tke(i,j) = (GGL90TKE(i,j,k,bi,bj)-surf_flx_tke(i,j))
                0982      &         *recip_deltaT
                0983          ENDDO
                0984         ENDDO
                0985         CALL DIAGNOSTICS_FILL( surf_flx_tke ,'GGL90Emn',
                0986      &                         k, 1, 2, bi, bj, myThid )
                0987        ENDIF
                0988 #endif
94c8eb5701 Jean*0989       ENDDO
004d5ee949 Davi*0990 
76f580e1f0 Mart*0991 C     end of time step
                0992 C     ===============================
004d5ee949 Davi*0993 
f688417df1 Jean*0994       DO k=2,Nr
f18a893d42 Mart*0995 #ifdef ALLOW_GGL90_SMOOTH
                0996        DO j=1-OLy,sNy+OLy
                0997         DO i=1-OLx,sNx+OLx
                0998          maskI(i,j) = maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
                0999      &        *mskCor(i,j,bi,bj)
                1000          GGL90visctmp(i,j,k) = GGL90visctmp(i,j,k)*mskCor(i,j,bi,bj)
                1001         ENDDO
                1002        ENDDO
                1003 #endif
f688417df1 Jean*1004        DO j=1,sNy
                1005         DO i=1,sNx
f6b150f7f1 Gael*1006 #ifdef ALLOW_GGL90_SMOOTH
63bbd437b1 Jean*1007          tmpVisc = (
f18a893d42 Mart*1008      &     p4 *    GGL90visctmp(i  ,j  ,k)
                1009      &    +p8 *( ( GGL90visctmp(i-1,j  ,k) + GGL90visctmp(i+1,j  ,k) )
                1010      &         + ( GGL90visctmp(i  ,j-1,k) + GGL90visctmp(i  ,j+1,k) ) )
                1011      &    +p16*( ( GGL90visctmp(i+1,j+1,k) + GGL90visctmp(i-1,j-1,k) )
                1012      &         + ( GGL90visctmp(i+1,j-1,k) + GGL90visctmp(i-1,j+1,k) ) )
63bbd437b1 Jean*1013      &             )/(
                1014      &     p4
f18a893d42 Mart*1015      &    +p8 *(( maskI(i-1,j  ) + maskI(i+1,j  ) )
                1016      &         +( maskI(i  ,j-1) + maskI(i  ,j+1) ) )
                1017      &    +p16*(( maskI(i+1,j+1) + maskI(i-1,j-1) )
                1018      &         +( maskI(i+1,j-1) + maskI(i-1,j+1) ) )
                1019      &               )*maskI(i,j)
f6b150f7f1 Gael*1020 #else
f688417df1 Jean*1021          tmpVisc = GGL90visctmp(i,j,k)
f6b150f7f1 Gael*1022 #endif
                1023          tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
0320e25227 Mart*1024      &        * coordFac*coordFac
78524d1402 Jean*1025          GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrS(k) )
f6b150f7f1 Gael*1026         ENDDO
                1027        ENDDO
                1028 
f688417df1 Jean*1029        DO j=1,sNy
                1030         DO i=1,sNx+1
f6b150f7f1 Gael*1031 #ifdef ALLOW_GGL90_SMOOTH
63bbd437b1 Jean*1032          tmpVisc = (
f18a893d42 Mart*1033      &     p4 *(   GGL90visctmp(i-1,j  ,k) + GGL90visctmp(i,j  ,k) )
                1034      &    +p8 *( ( GGL90visctmp(i-1,j-1,k) + GGL90visctmp(i,j-1,k) )
                1035      &         + ( GGL90visctmp(i-1,j+1,k) + GGL90visctmp(i,j+1,k) ) )
63bbd437b1 Jean*1036      &             )/(
                1037      &     p4 * 2. _d 0
f18a893d42 Mart*1038      &    +p8 *(( maskI(i-1,j-1) + maskI(i,j-1) )
                1039      &         +( maskI(i-1,j+1) + maskI(i,j+1) ) )
                1040      &               )*maskI(i-1,j) * maskI(i,j)
f6b150f7f1 Gael*1041 #else
f18a893d42 Mart*1042          tmpVisc = _maskW(i,j,k-1,bi,bj) * _maskW(i,j,k,bi,bj) * halfRL
63bbd437b1 Jean*1043      &          *( GGL90visctmp(i-1,j,k)
0320e25227 Mart*1044      &           + GGL90visctmp(i,  j,k) )
f6b150f7f1 Gael*1045 #endif
63bbd437b1 Jean*1046          tmpVisc = MIN( tmpVisc , GGL90viscMax )
0320e25227 Mart*1047      &        * coordFac*coordFac
63bbd437b1 Jean*1048          GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
004d5ee949 Davi*1049         ENDDO
                1050        ENDDO
f6b150f7f1 Gael*1051 
f688417df1 Jean*1052        DO j=1,sNy+1
                1053         DO i=1,sNx
f6b150f7f1 Gael*1054 #ifdef ALLOW_GGL90_SMOOTH
63bbd437b1 Jean*1055          tmpVisc = (
f18a893d42 Mart*1056      &     p4 *(   GGL90visctmp(i  ,j-1,k) + GGL90visctmp(i  ,j,k) )
                1057      &    +p8 *( ( GGL90visctmp(i-1,j-1,k) + GGL90visctmp(i-1,j,k) )
                1058      &         + ( GGL90visctmp(i+1,j-1,k) + GGL90visctmp(i+1,j,k) ) )
63bbd437b1 Jean*1059      &             )/(
                1060      &     p4 * 2. _d 0
f18a893d42 Mart*1061      &    +p8 *(( maskI(i-1,j-1) + maskI(i-1,j) )
                1062      &         +( maskI(i+1,j-1) + maskI(i+1,j) ) )
                1063      &               )*maskI(i,j-1) * maskI(i,j)
f6b150f7f1 Gael*1064 #else
f18a893d42 Mart*1065          tmpVisc = _maskS(i,j,k-1,bi,bj) * _maskS(i,j,k,bi,bj) * halfRL
63bbd437b1 Jean*1066      &          *( GGL90visctmp(i,j-1,k)
0320e25227 Mart*1067      &           + GGL90visctmp(i,j,  k) )
004d5ee949 Davi*1068 #endif
63bbd437b1 Jean*1069          tmpVisc = MIN( tmpVisc , GGL90viscMax )
0320e25227 Mart*1070      &        * coordFac*coordFac
63bbd437b1 Jean*1071          GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
f6b150f7f1 Gael*1072         ENDDO
                1073        ENDDO
                1074       ENDDO
73c2f90ab3 Davi*1075 
                1076 #ifdef ALLOW_DIAGNOSTICS
                1077       IF ( useDiagnostics ) THEN
305c472a49 Jean*1078         CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
                1079      &                         0,Nr, 1, bi, bj, myThid )
                1080         CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
                1081      &                         0,Nr, 1, bi, bj, myThid )
                1082         CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
                1083      &                         0,Nr, 1, bi, bj, myThid )
                1084         CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
                1085      &                         0,Nr, 1, bi, bj, myThid )
                1086         CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
                1087      &                         0,Nr, 2, bi, bj, myThid )
9293d3c672 Hajo*1088 #ifdef ALLOW_GGL90_LANGMUIR
                1089         IF (useLANGMUIR) THEN
                1090          CALL DIAGNOSTICS_FILL( LCmixingLength,'GGL90Lmx',
                1091      &                          0,Nr, 2, bi, bj, myThid )
                1092         ELSE
                1093 #else
                1094         IF (.TRUE.) THEN
                1095 #endif /* ALLOW_GGL90_LANGMUIR */
                1096          CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
                1097      &                          0,Nr, 2, bi, bj, myThid )
                1098         ENDIF
305c472a49 Jean*1099 
5b0716a6b3 Mart*1100 C     avoid extra 3D diagnostics field and abuse unused field
                1101         IF ( DIAGNOSTICS_IS_ON('GGL90KN2',myThid) ) THEN
                1102          DO k=1,Nr
                1103           DO j=1,sNy
                1104            DO i=1,sNx
                1105             TKEPrandtlNumber(i,j,k) =
                1106      &           GGL90diffKr(i,j,k,bi,bj) * Nsquare(i,j,k)
                1107            ENDDO
                1108           ENDDO
                1109          ENDDO
                1110          CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90KN2',
                1111      &                          0, Nr, 2, bi, bj, myThid )
                1112         ENDIF
                1113 
                1114         IF ( DIAGNOSTICS_IS_ON('GGL90flx',myThid) ) THEN
305c472a49 Jean*1115 C     diagnose surface flux of TKE
5b0716a6b3 Mart*1116          IF ( usingPCoords ) THEN
                1117           DO j=jMin,jMax
                1118            DO i=iMin,iMax
f18a893d42 Mart*1119 CML#ifdef ALLOW_SHELFICE
                1120 CML           IF ( useShelfIce ) THEN
                1121 CML            kSrf = MAX(kTopC(i,j,bi,bj),1)
                1122 CML            kTop = kSrf
                1123 CML           ENDIF
                1124 CML#endif
5b0716a6b3 Mart*1125             surf_flx_tke(i,j) =
                1126      &           (MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare(i,j))
                1127      &           - GGL90TKE(i,j,kSrf,bi,bj) )
                1128      &           *recip_drF(kSrf)*recip_hFacC(i,j,kSrf,bi,bj)
                1129      &           *KappaE(i,j,kSrf)
                1130      &           *coordFac
                1131            ENDDO
0320e25227 Mart*1132           ENDDO
5b0716a6b3 Mart*1133          ELSE
                1134           DO j=jMin,jMax
                1135            DO i=iMin,iMax
f18a893d42 Mart*1136 #ifdef ALLOW_SHELFICE
5b0716a6b3 Mart*1137             IF ( useShelfIce ) THEN
                1138              kSrf = MAX(1,kTopC(i,j,bi,bj))
                1139              kTop = MIN(kSrf+1,Nr)
                1140             ENDIF
f18a893d42 Mart*1141 #endif
5b0716a6b3 Mart*1142             surf_flx_tke(i,j) =(GGL90TKE(i,j,kSrf,bi,bj)-
                1143      &                          GGL90TKE(i,j,kTop,bi,bj))
0320e25227 Mart*1144      &         *recip_drF(kSrf)*recip_hFacC(i,j,kSrf,bi,bj)
                1145      &         *KappaE(i,j,kTop)
5b0716a6b3 Mart*1146            ENDDO
0320e25227 Mart*1147           ENDDO
5b0716a6b3 Mart*1148          ENDIF
                1149          CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90flx',
                1150      &                          0, 1, 2, bi, bj, myThid )
0320e25227 Mart*1151         ENDIF
31a3206180 Mart*1152 
5b0716a6b3 Mart*1153         IF ( DIAGNOSTICS_IS_ON('GGL90tau',myThid) ) THEN
                1154          k=kSrf
                1155          DO j=jMin,jMax
                1156           DO i=iMin,iMax
f18a893d42 Mart*1157 #ifdef ALLOW_SHELFICE
5b0716a6b3 Mart*1158            IF ( useShelfIce ) k = MAX(1,kTopC(i,j,bi,bj))
f18a893d42 Mart*1159 #endif
305c472a49 Jean*1160 C     diagnose work done by the wind
5b0716a6b3 Mart*1161            surf_flx_tke(i,j) =
305c472a49 Jean*1162      &      halfRL*( surfaceForcingU(i,  j,bi,bj)*uVel(i  ,j,k,bi,bj)
                1163      &              +surfaceForcingU(i+1,j,bi,bj)*uVel(i+1,j,k,bi,bj))
                1164      &    + halfRL*( surfaceForcingV(i,j,  bi,bj)*vVel(i,j  ,k,bi,bj)
                1165      &              +surfaceForcingV(i,j+1,bi,bj)*vVel(i,j+1,k,bi,bj))
5b0716a6b3 Mart*1166            surf_flx_tke(i,j) = surf_flx_tke(i,j) *recip_coordFac
                1167           ENDDO
305c472a49 Jean*1168          ENDDO
5b0716a6b3 Mart*1169          CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90tau',
                1170      &                          0, 1, 2, bi, bj, myThid )
                1171         ENDIF
                1172 C     endif useDiagnostics
73c2f90ab3 Davi*1173       ENDIF
31a3206180 Mart*1174 #endif /* ALLOW_DIAGNOSTICS */
89474f9a5c Mart*1175 
                1176 #endif /* ALLOW_GGL90 */
                1177 
                1178       RETURN
                1179       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!