Back to home page

MITgcm
 
 

    

File indexing completed on 2023-12-14 06:10:51 UTC

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