pkg/ggl90/ggl90_calc.F

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