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
0007
0008
0009
f688417df1 Jean*0010 SUBROUTINE GGL90_CALC(
dc6107c029 Jean*0011 I bi, bj, sigmaR, myTime, myIter, myThid )
0012
89474f9a5c Mart*0013
5e48dccc42 Jean*0014
89474f9a5c Mart*0015
0016
5e48dccc42 Jean*0017
89474f9a5c Mart*0018
0320e25227 Mart*0019
89474f9a5c Mart*0020
0320e25227 Mart*0021
0022
0023
0024
5e48dccc42 Jean*0025
89474f9a5c Mart*0026
0027
5e48dccc42 Jean*0028
0029
0030
89474f9a5c Mart*0031
0032
0033
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
5e48dccc42 Jean*0050
dc6107c029 Jean*0051
0052
5e48dccc42 Jean*0053
f688417df1 Jean*0054
5e48dccc42 Jean*0055
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
f688417df1 Jean*0070
0320e25227 Mart*0071
0072
0073
0074
cdafb98dea Mart*0075
f688417df1 Jean*0076
0077
0320e25227 Mart*0078
0079
0080
0081
0082
f688417df1 Jean*0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
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
0106 _RL Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0107
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
0123
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
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
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
0146
0147
0148 _RL coordFac, recip_coordFac
76f580e1f0 Mart*0149 #ifdef ALLOW_GGL90_HORIZDIFF
f688417df1 Jean*0150
0151
0152
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
0182
0183 INTEGER tkey, kkey
dd9d13d532 Mart*0184 #endif
31a3206180 Mart*0185
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
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
0224
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
cdafb98dea Mart*0237
0238
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
0258
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
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
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
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
dc6107c029 Jean*0344 Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
0320e25227 Mart*0345 & * sigmaR(i,j,k) * coordFac
cdafb98dea Mart*0346
0347
b038e3cc4f Mart*0348
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
b038e3cc4f Mart*0366
0367 # ifdef ALLOW_GGL90_LANGMUIR
0368
0369 # endif
0320e25227 Mart*0370 #endif
9293d3c672 Hajo*0371
b038e3cc4f Mart*0372
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
0381
76f580e1f0 Mart*0382
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
0394
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
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
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
f688417df1 Jean*0430 ENDIF
0431 #endif /* ALLOW_GGL90_HORIZDIFF */
0432
cdafb98dea Mart*0433
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
0457
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
0465 IF ( calcMeanVertShear ) THEN
0466
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
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
0498
0499 #endif /* ALLOW_AUTODIFF_TAMC */
63bbd437b1 Jean*0500
9293d3c672 Hajo*0501 #ifdef ALLOW_GGL90_LANGMUIR
0502 IF (useLANGMUIR) THEN
0503
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
0518
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
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
9293d3c672 Hajo*0555 # endif
0556 #endif /* ALLOW_GGL90_LANGMUIR */
0557
5b0716a6b3 Mart*0558
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
0588 #endif /* ALLOW_AUTODIFF_TAMC */
31a3206180 Mart*0589
cdafb98dea Mart*0590 DO j=jMin,jMax
0591 DO i=iMin,iMax
305c472a49 Jean*0592
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
0598 TKEdissipation = explDissFac*GGL90ceps
0599 & *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
0600 & *GGL90TKE(i,j,k,bi,bj)
0601
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
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
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
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
0648 ENDDO
0320e25227 Mart*0649 IF ( usingPCoords ) THEN
0650
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
76f580e1f0 Mart*0659
0320e25227 Mart*0660
0661
0662
0663
89474f9a5c Mart*0664
76f580e1f0 Mart*0665
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
0686
0687
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
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
0708
0709
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
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
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
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
0759
0760
31f96e9372 Jean*0761 IF ( calcMeanVertShear ) THEN
0762
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
b038e3cc4f Mart*0774
31f96e9372 Jean*0775
0776
b038e3cc4f Mart*0777
31f96e9372 Jean*0778
0779 ENDDO
0780 ENDDO
0781 ELSE
0782 DO j=jMin,jMax
0783 DO i=iMin,iMax
89474f9a5c Mart*0784
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
0795
0796
0797
0798 IF ( useSHELFICE .AND.
0799 & ( no_slip_shelfice .OR. SHELFICEDragLinear.NE.zeroRL
0800 & .OR. SHELFICEselectDragQuadr.GE.0 )
0801 & ) THEN
0802
0803
0804
0805
0806
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
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
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
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
0873
0874
0875
0876
0877 #endif
f18a893d42 Mart*0878
0320e25227 Mart*0879
0880 IF ( usingPCoords ) THEN
0881 DO j=jMin,jMax
0882 DO i=iMin,iMax
f18a893d42 Mart*0883
0884
0885
0886
0887
0888
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
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
dd9d13d532 Mart*0938 #endif
f688417df1 Jean*0939
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
dd9d13d532 Mart*0949 #endif
0320e25227 Mart*0950 DO k=2,Nr
5b0716a6b3 Mart*0951 #if ( defined ALLOW_DIAGNOSTICS &&
0952
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
0965
0966
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 &&
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
0987
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
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
5b0716a6b3 Mart*1111 IF ( usingPCoords ) THEN
1112 DO j=jMin,jMax
1113 DO i=iMin,iMax
f18a893d42 Mart*1114
1115
1116
1117
1118
1119
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
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
73c2f90ab3 Davi*1168 ENDIF
31a3206180 Mart*1169 #endif /* ALLOW_DIAGNOSTICS */
89474f9a5c Mart*1170
1171 #endif /* ALLOW_GGL90 */
1172
1173 RETURN
1174 END