8 |
|
|
9 |
C !INTERFACE: ====================================================== |
C !INTERFACE: ====================================================== |
10 |
SUBROUTINE GGL90_CALC( |
SUBROUTINE GGL90_CALC( |
11 |
I bi, bj, myTime, myIter, myThid ) |
I bi, bj, sigmaR, myTime, myIter, myThid ) |
12 |
|
|
13 |
|
|
14 |
C !DESCRIPTION: \bv |
C !DESCRIPTION: \bv |
15 |
C *==========================================================* |
C *==========================================================* |
42 |
|
|
43 |
C !INPUT PARAMETERS: =================================================== |
C !INPUT PARAMETERS: =================================================== |
44 |
C Routine arguments |
C Routine arguments |
45 |
C bi, bj :: array indices on which to apply calculations |
C bi, bj :: Current tile indices |
46 |
|
C sigmaR :: Vertical gradient of iso-neutral density |
47 |
C myTime :: Current time in simulation |
C myTime :: Current time in simulation |
48 |
C myIter :: Current time-step number |
C myIter :: Current time-step number |
49 |
C myThid :: My Thread Id number |
C myThid :: My Thread Id number |
50 |
INTEGER bi, bj |
INTEGER bi, bj |
51 |
|
_RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
52 |
_RL myTime |
_RL myTime |
53 |
INTEGER myIter |
INTEGER myIter |
54 |
INTEGER myThid |
INTEGER myThid |
75 |
C rMixingLength:: inverse of mixing length |
C rMixingLength:: inverse of mixing length |
76 |
C totalDepth :: thickness of water column (inverse of recip_Rcol) |
C totalDepth :: thickness of water column (inverse of recip_Rcol) |
77 |
C TKEPrandtlNumber :: here, an empirical function of the Richardson number |
C TKEPrandtlNumber :: here, an empirical function of the Richardson number |
|
C rhoK, rhoKm1 :: density at layer k and km1 (relative to k) |
|
78 |
INTEGER iMin ,iMax ,jMin ,jMax |
INTEGER iMin ,iMax ,jMin ,jMax |
79 |
INTEGER i, j, k, kp1, km1, kSurf, kBottom |
INTEGER i, j, k, kp1, km1, kSurf, kBottom |
80 |
_RL explDissFac, implDissFac |
_RL explDissFac, implDissFac |
95 |
_RL rMixingLength (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
_RL rMixingLength (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
96 |
_RL mxLength_Dn (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
_RL mxLength_Dn (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
97 |
_RL KappaE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
_RL KappaE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
|
_RL rhoK (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
|
|
_RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
|
98 |
_RL totalDepth (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL totalDepth (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
99 |
_RL GGL90visctmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
_RL GGL90visctmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
100 |
C- tri-diagonal matrix |
C- tri-diagonal matrix |
133 |
|
|
134 |
C Initialize local fields |
C Initialize local fields |
135 |
DO k = 1, Nr |
DO k = 1, Nr |
136 |
DO j=1-Oly,sNy+Oly |
DO j=1-OLy,sNy+OLy |
137 |
DO i=1-Olx,sNx+Olx |
DO i=1-OLx,sNx+OLx |
138 |
KappaE(i,j,k) = 0. _d 0 |
KappaE(i,j,k) = 0. _d 0 |
139 |
TKEPrandtlNumber(i,j,k) = 1. _d 0 |
TKEPrandtlNumber(i,j,k) = 1. _d 0 |
140 |
GGL90mixingLength(i,j,k) = GGL90mixingLengthMin |
GGL90mixingLength(i,j,k) = GGL90mixingLengthMin |
147 |
ENDDO |
ENDDO |
148 |
ENDDO |
ENDDO |
149 |
ENDDO |
ENDDO |
150 |
DO j=1-Oly,sNy+Oly |
DO j=1-OLy,sNy+OLy |
151 |
DO i=1-Olx,sNx+Olx |
DO i=1-OLx,sNx+OLx |
|
rhoK(i,j) = 0. _d 0 |
|
|
rhoKm1(i,j) = 0. _d 0 |
|
152 |
totalDepth(i,j) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) |
totalDepth(i,j) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) |
153 |
rMixingLength(i,j,1) = 0. _d 0 |
rMixingLength(i,j,1) = 0. _d 0 |
154 |
mxLength_Dn(i,j,1) = GGL90mixingLengthMin |
mxLength_Dn(i,j,1) = GGL90mixingLengthMin |
158 |
|
|
159 |
C start k-loop |
C start k-loop |
160 |
DO k = 2, Nr |
DO k = 2, Nr |
161 |
km1 = k-1 |
c km1 = k-1 |
162 |
c kp1 = MIN(Nr,k+1) |
c kp1 = MIN(Nr,k+1) |
|
CALL FIND_RHO_2D( |
|
|
I iMin, iMax, jMin, jMax, k, |
|
|
I theta(1-OLx,1-OLy,km1,bi,bj), salt(1-OLx,1-OLy,km1,bi,bj), |
|
|
O rhoKm1, |
|
|
I km1, bi, bj, myThid ) |
|
|
|
|
|
CALL FIND_RHO_2D( |
|
|
I iMin, iMax, jMin, jMax, k, |
|
|
I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj), |
|
|
O rhoK, |
|
|
I k, bi, bj, myThid ) |
|
163 |
DO j=jMin,jMax |
DO j=jMin,jMax |
164 |
DO i=iMin,iMax |
DO i=iMin,iMax |
165 |
SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) ) |
SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) ) |
166 |
|
|
167 |
C buoyancy frequency |
C buoyancy frequency |
168 |
Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(k) |
Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst |
169 |
& * ( rhoKm1(i,j) - rhoK(i,j) )*maskC(i,j,k,bi,bj) |
& * sigmaR(i,j,k) |
170 |
cC vertical shear term (dU/dz)^2+(dV/dz)^2 |
cC vertical shear term (dU/dz)^2+(dV/dz)^2 |
171 |
c tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj) |
c tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj) |
172 |
c & -( uVel(i,j,k ,bi,bj)+uVel(i+1,j,k ,bi,bj)) ) |
c & -( uVel(i,j,k ,bi,bj)+uVel(i+1,j,k ,bi,bj)) ) |
187 |
ENDDO |
ENDDO |
188 |
ENDDO |
ENDDO |
189 |
|
|
190 |
|
C- ensure mixing between first and second level |
191 |
|
IF (mxlSurfFlag) THEN |
192 |
|
DO j=jMin,jMax |
193 |
|
DO i=iMin,iMax |
194 |
|
GGL90mixingLength(i,j,2)=drF(1) |
195 |
|
ENDDO |
196 |
|
ENDDO |
197 |
|
ENDIF |
198 |
|
|
199 |
C- Impose upper and lower bound for mixing length |
C- Impose upper and lower bound for mixing length |
200 |
IF ( mxlMaxFlag .EQ. 0 ) THEN |
IF ( mxlMaxFlag .EQ. 0 ) THEN |
201 |
C- |
|
202 |
DO k=2,Nr |
DO k=2,Nr |
203 |
DO j=jMin,jMax |
DO j=jMin,jMax |
204 |
DO i=iMin,iMax |
DO i=iMin,iMax |
220 |
ENDDO |
ENDDO |
221 |
|
|
222 |
ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN |
ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN |
223 |
C- |
|
224 |
DO k=2,Nr |
DO k=2,Nr |
225 |
DO j=jMin,jMax |
DO j=jMin,jMax |
226 |
DO i=iMin,iMax |
DO i=iMin,iMax |
243 |
ENDDO |
ENDDO |
244 |
|
|
245 |
ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN |
ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN |
246 |
C- |
|
|
cgf ensure mixing between first and second level |
|
|
c DO j=jMin,jMax |
|
|
c DO i=iMin,iMax |
|
|
c GGL90mixingLength(i,j,2)=drF(1) |
|
|
c ENDDO |
|
|
c ENDDO |
|
|
cgf |
|
247 |
DO k=2,Nr |
DO k=2,Nr |
248 |
DO j=jMin,jMax |
DO j=jMin,jMax |
249 |
DO i=iMin,iMax |
DO i=iMin,iMax |
278 |
ENDDO |
ENDDO |
279 |
|
|
280 |
ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN |
ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN |
281 |
C- |
|
282 |
DO k=2,Nr |
DO k=2,Nr |
283 |
DO j=jMin,jMax |
DO j=jMin,jMax |
284 |
DO i=iMin,iMax |
DO i=iMin,iMax |
337 |
C horizontal diffusion of TKE (requires an exchange in |
C horizontal diffusion of TKE (requires an exchange in |
338 |
C do_fields_blocking_exchanges) |
C do_fields_blocking_exchanges) |
339 |
C common factors |
C common factors |
340 |
DO j=1-Oly,sNy+Oly |
DO j=1-OLy,sNy+OLy |
341 |
DO i=1-Olx,sNx+Olx |
DO i=1-OLx,sNx+OLx |
342 |
xA(i,j) = _dyG(i,j,bi,bj) |
xA(i,j) = _dyG(i,j,bi,bj) |
343 |
& *drF(k)*_hFacW(i,j,k,bi,bj) |
& *drF(k)*_hFacW(i,j,k,bi,bj) |
344 |
yA(i,j) = _dxG(i,j,bi,bj) |
yA(i,j) = _dxG(i,j,bi,bj) |
347 |
ENDDO |
ENDDO |
348 |
C Compute diffusive fluxes |
C Compute diffusive fluxes |
349 |
C ... across x-faces |
C ... across x-faces |
350 |
DO j=1-Oly,sNy+Oly |
DO j=1-OLy,sNy+OLy |
351 |
dfx(1-Olx,j)=0. _d 0 |
dfx(1-OLx,j)=0. _d 0 |
352 |
DO i=1-Olx+1,sNx+Olx |
DO i=1-OLx+1,sNx+OLx |
353 |
dfx(i,j) = -GGL90diffTKEh*xA(i,j) |
dfx(i,j) = -GGL90diffTKEh*xA(i,j) |
354 |
& *_recip_dxC(i,j,bi,bj) |
& *_recip_dxC(i,j,bi,bj) |
355 |
& *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj)) |
& *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj)) |
357 |
ENDDO |
ENDDO |
358 |
ENDDO |
ENDDO |
359 |
C ... across y-faces |
C ... across y-faces |
360 |
DO i=1-Olx,sNx+Olx |
DO i=1-OLx,sNx+OLx |
361 |
dfy(i,1-Oly)=0. _d 0 |
dfy(i,1-OLy)=0. _d 0 |
362 |
ENDDO |
ENDDO |
363 |
DO j=1-Oly+1,sNy+Oly |
DO j=1-OLy+1,sNy+OLy |
364 |
DO i=1-Olx,sNx+Olx |
DO i=1-OLx,sNx+OLx |
365 |
dfy(i,j) = -GGL90diffTKEh*yA(i,j) |
dfy(i,j) = -GGL90diffTKEh*yA(i,j) |
366 |
& *_recip_dyC(i,j,bi,bj) |
& *_recip_dyC(i,j,bi,bj) |
367 |
& *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj)) |
& *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj)) |
371 |
ENDDO |
ENDDO |
372 |
ENDDO |
ENDDO |
373 |
C Compute divergence of fluxes |
C Compute divergence of fluxes |
374 |
DO j=1-Oly,sNy+Oly-1 |
DO j=1-OLy,sNy+OLy-1 |
375 |
DO i=1-Olx,sNx+Olx-1 |
DO i=1-OLx,sNx+OLx-1 |
376 |
gTKE(i,j) = |
gTKE(i,j) = |
377 |
& -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj) |
& -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj) |
378 |
& *( (dfx(i+1,j)-dfx(i,j)) |
& *( (dfx(i+1,j)-dfx(i,j)) |