/[MITgcm]/MITgcm/pkg/ggl90/ggl90_calc.F
ViewVC logotype

Contents of /MITgcm/pkg/ggl90/ggl90_calc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.25 - (show annotations) (download)
Mon Nov 5 13:11:11 2012 UTC (11 years, 8 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint65, checkpoint65b, checkpoint65a
Changes since 1.24: +1 -4 lines
fix previous modification: remove spurious (and incorrect) masking

1 C $Header: /u/gcmpack/MITgcm/pkg/ggl90/ggl90_calc.F,v 1.24 2012/11/05 08:11:45 mlosch Exp $
2 C $Name: $
3
4 #include "GGL90_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: GGL90_CALC
8
9 C !INTERFACE: ======================================================
10 SUBROUTINE GGL90_CALC(
11 I bi, bj, sigmaR, myTime, myIter, myThid )
12
13
14 C !DESCRIPTION: \bv
15 C *==========================================================*
16 C | SUBROUTINE GGL90_CALC |
17 C | o Compute all GGL90 fields defined in GGL90.h |
18 C *==========================================================*
19 C | Equation numbers refer to |
20 C | Gaspar et al. (1990), JGR 95 (C9), pp 16,179 |
21 C | Some parts of the implementation follow Blanke and |
22 C | Delecuse (1993), JPO, and OPA code, in particular the |
23 C | computation of the |
24 C | mixing length = max(min(lk,depth),lkmin) |
25 C *==========================================================*
26
27 C global parameters updated by ggl90_calc
28 C GGL90TKE :: sub-grid turbulent kinetic energy (m^2/s^2)
29 C GGL90viscAz :: GGL90 eddy viscosity coefficient (m^2/s)
30 C GGL90diffKzT :: GGL90 diffusion coefficient for temperature (m^2/s)
31 C \ev
32
33 C !USES: ============================================================
34 IMPLICIT NONE
35 #include "SIZE.h"
36 #include "EEPARAMS.h"
37 #include "PARAMS.h"
38 #include "DYNVARS.h"
39 #include "GGL90.h"
40 #include "FFIELDS.h"
41 #include "GRID.h"
42
43 C !INPUT PARAMETERS: ===================================================
44 C Routine arguments
45 C bi, bj :: Current tile indices
46 C sigmaR :: Vertical gradient of iso-neutral density
47 C myTime :: Current time in simulation
48 C myIter :: Current time-step number
49 C myThid :: My Thread Id number
50 INTEGER bi, bj
51 _RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
52 _RL myTime
53 INTEGER myIter
54 INTEGER myThid
55 CEOP
56
57 #ifdef ALLOW_GGL90
58
59 C !LOCAL VARIABLES: ====================================================
60 C Local constants
61 C iMin,iMax,jMin,jMax :: index boundaries of computation domain
62 C i, j, k, kp1,km1 :: array computation indices
63 C kSurf, kBottom :: vertical indices of domain boundaries
64 C explDissFac :: explicit Dissipation Factor (in [0-1])
65 C implDissFac :: implicit Dissipation Factor (in [0-1])
66 C uStarSquare :: square of friction velocity
67 C verticalShear :: (squared) vertical shear of horizontal velocity
68 C Nsquare :: squared buoyancy freqency
69 C RiNumber :: local Richardson number
70 C KappaM :: (local) viscosity parameter (eq.10)
71 C KappaH :: (local) diffusivity parameter for temperature (eq.11)
72 C KappaE :: (local) diffusivity parameter for TKE (eq.15)
73 C TKEdissipation :: dissipation of TKE
74 C GGL90mixingLength:: mixing length of scheme following Banke+Delecuse
75 C rMixingLength:: inverse of mixing length
76 C totalDepth :: thickness of water column (inverse of recip_Rcol)
77 C TKEPrandtlNumber :: here, an empirical function of the Richardson number
78 INTEGER iMin ,iMax ,jMin ,jMax
79 INTEGER i, j, k, kp1, km1, kSurf, kBottom
80 _RL explDissFac, implDissFac
81 _RL uStarSquare
82 _RL verticalShear
83 _RL KappaM, KappaH
84 c _RL Nsquare
85 _RL Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86 _RL deltaTggl90
87 c _RL SQRTTKE
88 _RL SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89 _RL RiNumber
90 _RL TKEdissipation
91 _RL tempU, tempV, prTemp
92 _RL MaxLength, tmpmlx, tmpVisc
93 _RL TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
94 _RL GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
95 _RL rMixingLength (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96 _RL mxLength_Dn (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97 _RL KappaE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
98 _RL totalDepth (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99 _RL GGL90visctmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100 C- tri-diagonal matrix
101 _RL a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102 _RL b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103 _RL c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104 INTEGER errCode
105 #ifdef ALLOW_GGL90_HORIZDIFF
106 C hFac :: fractional thickness of W-cell
107 C xA, yA :: area of lateral faces
108 C dfx, dfy :: diffusive flux across lateral faces
109 C gTKE :: right hand side of diffusion equation
110 _RL hFac
111 _RL xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112 _RL yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113 _RL dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114 _RL dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115 _RL gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116 #endif /* ALLOW_GGL90_HORIZDIFF */
117 #ifdef ALLOW_GGL90_SMOOTH
118 _RL p4, p8, p16
119 p4=0.25 _d 0
120 p8=0.125 _d 0
121 p16=0.0625 _d 0
122 #endif
123 iMin = 2-OLx
124 iMax = sNx+OLx-1
125 jMin = 2-OLy
126 jMax = sNy+OLy-1
127
128 C set separate time step (should be deltaTtracer)
129 deltaTggl90 = dTtracerLev(1)
130
131 kSurf = 1
132 C explicit/implicit timestepping weights for dissipation
133 explDissFac = 0. _d 0
134 implDissFac = 1. _d 0 - explDissFac
135
136 C Initialize local fields
137 DO k = 1, Nr
138 DO j=1-OLy,sNy+OLy
139 DO i=1-OLx,sNx+OLx
140 KappaE(i,j,k) = 0. _d 0
141 TKEPrandtlNumber(i,j,k) = 1. _d 0
142 GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
143 GGL90visctmp(i,j,k) = 0. _d 0
144 #ifndef SOLVE_DIAGONAL_LOWMEMORY
145 a3d(i,j,k) = 0. _d 0
146 b3d(i,j,k) = 1. _d 0
147 c3d(i,j,k) = 0. _d 0
148 #endif
149 ENDDO
150 ENDDO
151 ENDDO
152 DO j=1-OLy,sNy+OLy
153 DO i=1-OLx,sNx+OLx
154 totalDepth(i,j) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
155 rMixingLength(i,j,1) = 0. _d 0
156 mxLength_Dn(i,j,1) = GGL90mixingLengthMin
157 SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
158 ENDDO
159 ENDDO
160
161 C start k-loop
162 DO k = 2, Nr
163 c km1 = k-1
164 c kp1 = MIN(Nr,k+1)
165 DO j=jMin,jMax
166 DO i=iMin,iMax
167 SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
168
169 C buoyancy frequency
170 Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
171 & * sigmaR(i,j,k)
172 cC vertical shear term (dU/dz)^2+(dV/dz)^2
173 c tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
174 c & -( uVel(i,j,k ,bi,bj)+uVel(i+1,j,k ,bi,bj)) )
175 c & *recip_drC(k)
176 c tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
177 c & -( vVel(i,j,k ,bi,bj)+vVel(i,j+1,k ,bi,bj)) )
178 c & *recip_drC(k)
179 c verticalShear = tempU*tempU + tempV*tempV
180 c RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
181 cC compute Prandtl number (always greater than 0)
182 c prTemp = 1. _d 0
183 c IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
184 c TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
185 C mixing length
186 GGL90mixingLength(i,j,k) = SQRTTWO *
187 & SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
188 ENDDO
189 ENDDO
190 ENDDO
191
192 C- ensure mixing between first and second level
193 IF (mxlSurfFlag) THEN
194 DO j=jMin,jMax
195 DO i=iMin,iMax
196 GGL90mixingLength(i,j,2)=drF(1)
197 ENDDO
198 ENDDO
199 ENDIF
200
201 C- Impose upper and lower bound for mixing length
202 IF ( mxlMaxFlag .EQ. 0 ) THEN
203
204 DO k=2,Nr
205 DO j=jMin,jMax
206 DO i=iMin,iMax
207 MaxLength=totalDepth(i,j)
208 GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
209 & MaxLength)
210 ENDDO
211 ENDDO
212 ENDDO
213
214 DO k=2,Nr
215 DO j=jMin,jMax
216 DO i=iMin,iMax
217 GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
218 & GGL90mixingLengthMin)
219 rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
220 ENDDO
221 ENDDO
222 ENDDO
223
224 ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
225
226 DO k=2,Nr
227 DO j=jMin,jMax
228 DO i=iMin,iMax
229 MaxLength=MIN(Ro_surf(i,j,bi,bj)-rF(k),rF(k)-R_low(i,j,bi,bj))
230 c MaxLength=MAX(MaxLength,20. _d 0)
231 GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
232 & MaxLength)
233 ENDDO
234 ENDDO
235 ENDDO
236
237 DO k=2,Nr
238 DO j=jMin,jMax
239 DO i=iMin,iMax
240 GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
241 & GGL90mixingLengthMin)
242 rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
243 ENDDO
244 ENDDO
245 ENDDO
246
247 ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
248
249 DO k=2,Nr
250 DO j=jMin,jMax
251 DO i=iMin,iMax
252 GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
253 & GGL90mixingLength(i,j,k-1)+drF(k-1))
254 ENDDO
255 ENDDO
256 ENDDO
257 DO j=jMin,jMax
258 DO i=iMin,iMax
259 GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
260 & GGL90mixingLengthMin+drF(Nr))
261 ENDDO
262 ENDDO
263 DO k=Nr-1,2,-1
264 DO j=jMin,jMax
265 DO i=iMin,iMax
266 GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
267 & GGL90mixingLength(i,j,k+1)+drF(k))
268 ENDDO
269 ENDDO
270 ENDDO
271
272 DO k=2,Nr
273 DO j=jMin,jMax
274 DO i=iMin,iMax
275 GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
276 & GGL90mixingLengthMin)
277 rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
278 ENDDO
279 ENDDO
280 ENDDO
281
282 ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
283
284 DO k=2,Nr
285 DO j=jMin,jMax
286 DO i=iMin,iMax
287 mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
288 & mxLength_Dn(i,j,k-1)+drF(k-1))
289 ENDDO
290 ENDDO
291 ENDDO
292 DO j=jMin,jMax
293 DO i=iMin,iMax
294 GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
295 & GGL90mixingLengthMin+drF(Nr))
296 ENDDO
297 ENDDO
298 DO k=Nr-1,2,-1
299 DO j=jMin,jMax
300 DO i=iMin,iMax
301 GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
302 & GGL90mixingLength(i,j,k+1)+drF(k))
303 ENDDO
304 ENDDO
305 ENDDO
306
307 DO k=2,Nr
308 DO j=jMin,jMax
309 DO i=iMin,iMax
310 GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
311 & mxLength_Dn(i,j,k))
312 tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
313 tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
314 rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
315 ENDDO
316 ENDDO
317 ENDDO
318
319 ELSE
320 STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
321 ENDIF
322
323 C- Impose minimum mixing length (to avoid division by zero)
324 c DO k=2,Nr
325 c DO j=jMin,jMax
326 c DO i=iMin,iMax
327 c GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
328 c & GGL90mixingLengthMin)
329 c rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
330 c ENDDO
331 c ENDDO
332 c ENDDO
333
334 DO k=2,Nr
335 km1 = k-1
336
337 #ifdef ALLOW_GGL90_HORIZDIFF
338 IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
339 C horizontal diffusion of TKE (requires an exchange in
340 C do_fields_blocking_exchanges)
341 C common factors
342 DO j=1-OLy,sNy+OLy
343 DO i=1-OLx,sNx+OLx
344 xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
345 & (min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
346 & min(.5 _d 0,_hFacW(i,j,k ,bi,bj) ) )
347 yA(i,j) = _dxG(i,j,bi,bj)*drC(k)*
348 & (min(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
349 & min(.5 _d 0,_hFacS(i,j,k ,bi,bj) ) )
350 ENDDO
351 ENDDO
352 C Compute diffusive fluxes
353 C ... across x-faces
354 DO j=1-OLy,sNy+OLy
355 dfx(1-OLx,j)=0. _d 0
356 DO i=1-OLx+1,sNx+OLx
357 dfx(i,j) = -GGL90diffTKEh*xA(i,j)
358 & *_recip_dxC(i,j,bi,bj)
359 & *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
360 #ifdef ISOTROPIC_COS_SCALING
361 & *CosFacU(j,bi,bj)
362 #endif /* ISOTROPIC_COS_SCALING */
363 ENDDO
364 ENDDO
365 C ... across y-faces
366 DO i=1-OLx,sNx+OLx
367 dfy(i,1-OLy)=0. _d 0
368 ENDDO
369 DO j=1-OLy+1,sNy+OLy
370 DO i=1-OLx,sNx+OLx
371 dfy(i,j) = -GGL90diffTKEh*yA(i,j)
372 & *_recip_dyC(i,j,bi,bj)
373 & *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
374 #ifdef ISOTROPIC_COS_SCALING
375 & *CosFacV(j,bi,bj)
376 #endif /* ISOTROPIC_COS_SCALING */
377 ENDDO
378 ENDDO
379 C Compute divergence of fluxes
380 DO j=1-OLy,sNy+OLy-1
381 DO i=1-OLx,sNx+OLx-1
382 hFac = min(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +
383 & min(.5 _d 0,_hFacC(i,j,k ,bi,bj) )
384 gTKE(i,j) = 0.0
385 if ( hFac .ne. 0.0 )
386 & gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac
387 & *((dfx(i+1,j)-dfx(i,j))
388 & +(dfy(i,j+1)-dfy(i,j)) )
389 ENDDO
390 ENDDO
391 C end if GGL90diffTKEh .eq. 0.
392 ENDIF
393 #endif /* ALLOW_GGL90_HORIZDIFF */
394
395 DO j=jMin,jMax
396 DO i=iMin,iMax
397 C vertical shear term (dU/dz)^2+(dV/dz)^2
398 tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
399 & -( uVel(i,j,k ,bi,bj)+uVel(i+1,j,k ,bi,bj)) )
400 & *recip_drC(k)
401 tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
402 & -( vVel(i,j,k ,bi,bj)+vVel(i,j+1,k ,bi,bj)) )
403 & *recip_drC(k)
404 verticalShear = tempU*tempU + tempV*tempV
405 RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
406 C compute Prandtl number (always greater than 0)
407 prTemp = 1. _d 0
408 IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
409 TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
410 c TKEPrandtlNumber(i,j,k) = 1. _d 0
411
412 C viscosity and diffusivity
413 KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
414 GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
415 & * maskC(i,j,k,bi,bj)
416 c note: storing GGL90visctmp like this, and using it later to compute
417 c GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
418 KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
419 KappaH = KappaM/TKEPrandtlNumber(i,j,k)
420 KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
421
422 C dissipation term
423 TKEdissipation = explDissFac*GGL90ceps
424 & *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
425 & *GGL90TKE(i,j,k,bi,bj)
426 C partial update with sum of explicit contributions
427 GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
428 & + deltaTggl90*(
429 & + KappaM*verticalShear
430 & - KappaH*Nsquare(i,j,k)
431 & - TKEdissipation
432 & )
433 ENDDO
434 ENDDO
435
436 #ifdef ALLOW_GGL90_HORIZDIFF
437 IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
438 C-- Add horiz. diffusion tendency
439 DO j=jMin,jMax
440 DO i=iMin,iMax
441 GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
442 & + gTKE(i,j)*deltaTggl90
443 ENDDO
444 ENDDO
445 ENDIF
446 #endif /* ALLOW_GGL90_HORIZDIFF */
447
448 C-- end of k loop
449 ENDDO
450
451 C ============================================
452 C Implicit time step to update TKE for k=1,Nr;
453 C TKE(Nr+1)=0 by default
454 C ============================================
455 C set up matrix
456 C-- Lower diagonal
457 DO j=jMin,jMax
458 DO i=iMin,iMax
459 a3d(i,j,1) = 0. _d 0
460 ENDDO
461 ENDDO
462 DO k=2,Nr
463 km1=MAX(2,k-1)
464 DO j=jMin,jMax
465 DO i=iMin,iMax
466 C- We keep recip_hFacC in the diffusive flux calculation,
467 C- but no hFacC in TKE volume control
468 C- No need for maskC(k-1) with recip_hFacC(k-1)
469 a3d(i,j,k) = -deltaTggl90
470 & *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
471 & *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
472 & *recip_drC(k)*maskC(i,j,k,bi,bj)
473 ENDDO
474 ENDDO
475 ENDDO
476 C-- Upper diagonal
477 DO j=jMin,jMax
478 DO i=iMin,iMax
479 c3d(i,j,1) = 0. _d 0
480 ENDDO
481 ENDDO
482 DO k=2,Nr
483 DO j=jMin,jMax
484 DO i=iMin,iMax
485 kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
486 C- We keep recip_hFacC in the diffusive flux calculation,
487 C- but no hFacC in TKE volume control
488 C- No need for maskC(k) with recip_hFacC(k)
489 c3d(i,j,k) = -deltaTggl90
490 & *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
491 & *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
492 & *recip_drC(k)*maskC(i,j,k-1,bi,bj)
493 ENDDO
494 ENDDO
495 ENDDO
496 C-- Center diagonal
497 DO k=1,Nr
498 km1 = MAX(k-1,1)
499 DO j=jMin,jMax
500 DO i=iMin,iMax
501 b3d(i,j,k) = 1. _d 0 - c3d(i,j,k) - a3d(i,j,k)
502 & + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
503 & * rMixingLength(i,j,k)
504 & * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
505 ENDDO
506 ENDDO
507 ENDDO
508 C end set up matrix
509
510 C Apply boundary condition
511 kp1 = MIN(Nr,kSurf+1)
512 DO j=jMin,jMax
513 DO i=iMin,iMax
514 C estimate friction velocity uStar from surface forcing
515 uStarSquare = SQRT(
516 & ( .5 _d 0*( surfaceForcingU(i, j, bi,bj)
517 & + surfaceForcingU(i+1,j, bi,bj) ) )**2
518 & + ( .5 _d 0*( surfaceForcingV(i, j, bi,bj)
519 & + surfaceForcingV(i, j+1,bi,bj) ) )**2
520 & )
521 C Dirichlet surface boundary condition for TKE
522 GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
523 & *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
524 GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
525 & - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
526 a3d(i,j,kp1) = 0. _d 0
527 C Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
528 kBottom = MAX(kLowC(i,j,bi,bj),1)
529 GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
530 & - GGL90TKEbottom*c3d(i,j,kBottom)
531 c3d(i,j,kBottom) = 0. _d 0
532 ENDDO
533 ENDDO
534
535 C solve tri-diagonal system
536 CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
537 I a3d, b3d, c3d,
538 U GGL90TKE,
539 O errCode,
540 I bi, bj, myThid )
541
542 DO k=1,Nr
543 DO j=jMin,jMax
544 DO i=iMin,iMax
545 C impose minimum TKE to avoid numerical undershoots below zero
546 GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
547 & *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
548 ENDDO
549 ENDDO
550 ENDDO
551
552 C end of time step
553 C ===============================
554
555 DO k=2,Nr
556 DO j=1,sNy
557 DO i=1,sNx
558 #ifdef ALLOW_GGL90_SMOOTH
559 tmpVisc=
560 & (
561 & p4 * GGL90visctmp(i ,j ,k) * mskCor(i ,j ,bi,bj)
562 & +p8 *( GGL90visctmp(i-1,j ,k) * mskCor(i-1,j ,bi,bj)
563 & + GGL90visctmp(i ,j-1,k) * mskCor(i ,j-1,bi,bj)
564 & + GGL90visctmp(i+1,j ,k) * mskCor(i+1,j ,bi,bj)
565 & + GGL90visctmp(i ,j+1,k) * mskCor(i ,j+1,bi,bj))
566 & +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
567 & + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
568 & + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
569 & + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
570 & )
571 & /(p4
572 & +p8 *( maskC(i-1,j ,k,bi,bj) * mskCor(i-1,j ,bi,bj)
573 & + maskC(i ,j-1,k,bi,bj) * mskCor(i ,j-1,bi,bj)
574 & + maskC(i+1,j ,k,bi,bj) * mskCor(i+1,j ,bi,bj)
575 & + maskC(i ,j+1,k,bi,bj) * mskCor(i ,j+1,bi,bj))
576 & +p16*( maskC(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)
577 & + maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)
578 & + maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
579 & + maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))
580 & )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
581 #else
582 tmpVisc = GGL90visctmp(i,j,k)
583 #endif
584 tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
585 GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
586 ENDDO
587 ENDDO
588 ENDDO
589
590 DO k=2,Nr
591 DO j=1,sNy
592 DO i=1,sNx+1
593 #ifdef ALLOW_GGL90_SMOOTH
594 tmpVisc =
595 & (
596 & p4 *(GGL90visctmp(i ,j ,k) * mskCor(i ,j ,bi,bj)
597 & +GGL90visctmp(i-1,j ,k) * mskCor(i-1,j ,bi,bj))
598 & +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
599 & +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
600 & +GGL90visctmp(i ,j-1,k) * mskCor(i ,j-1,bi,bj)
601 & +GGL90visctmp(i ,j+1,k) * mskCor(i ,j+1,bi,bj))
602 & )
603 & /(p4 * 2. _d 0
604 & +p8 *( maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
605 & + maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
606 & + maskC(i ,j-1,k,bi,bj) * mskCor(i ,j-1,bi,bj)
607 & + maskC(i ,j+1,k,bi,bj) * mskCor(i ,j+1,bi,bj))
608 & )
609 & *maskC(i ,j,k,bi,bj)*mskCor(i ,j,bi,bj)
610 & *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
611 #else
612 tmpVisc = _maskW(i,j,k,bi,bj) *
613 & (.5 _d 0*(GGL90visctmp(i,j,k)
614 & +GGL90visctmp(i-1,j,k))
615 & )
616 #endif
617 tmpVisc = MIN( tmpVisc , GGL90viscMax )
618 GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
619 ENDDO
620 ENDDO
621 ENDDO
622
623 DO k=2,Nr
624 DO j=1,sNy+1
625 DO i=1,sNx
626 #ifdef ALLOW_GGL90_SMOOTH
627 tmpVisc =
628 & (
629 & p4 *(GGL90visctmp(i ,j ,k) * mskCor(i ,j ,bi,bj)
630 & +GGL90visctmp(i ,j-1,k) * mskCor(i ,j-1,bi,bj))
631 & +p8 *(GGL90visctmp(i-1,j ,k) * mskCor(i-1,j ,bi,bj)
632 & +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
633 & +GGL90visctmp(i+1,j ,k) * mskCor(i+1,j ,bi,bj)
634 & +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))
635 & )
636 & /(p4 * 2. _d 0
637 & +p8 *( maskC(i-1,j ,k,bi,bj) * mskCor(i-1,j ,bi,bj)
638 & + maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
639 & + maskC(i+1,j ,k,bi,bj) * mskCor(i+1,j ,bi,bj)
640 & + maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))
641 & )
642 & *maskC(i,j ,k,bi,bj)*mskCor(i,j ,bi,bj)
643 & *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
644 #else
645 tmpVisc = _maskS(i,j,k,bi,bj) *
646 & (.5 _d 0*(GGL90visctmp(i,j,k)
647 & +GGL90visctmp(i,j-1,k))
648 & )
649
650 #endif
651 tmpVisc = MIN( tmpVisc , GGL90viscMax )
652 GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
653 ENDDO
654 ENDDO
655 ENDDO
656
657 #ifdef ALLOW_DIAGNOSTICS
658 IF ( useDiagnostics ) THEN
659 CALL DIAGNOSTICS_FILL( GGL90TKE ,'GGL90TKE',
660 & 0,Nr, 1, bi, bj, myThid )
661 CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
662 & 0,Nr, 1, bi, bj, myThid )
663 CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
664 & 0,Nr, 1, bi, bj, myThid )
665 CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
666 & 0,Nr, 1, bi, bj, myThid )
667 CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
668 & 0,Nr, 2, bi, bj, myThid )
669 CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
670 & 0,Nr, 2, bi, bj, myThid )
671 ENDIF
672 #endif
673
674 #endif /* ALLOW_GGL90 */
675
676 RETURN
677 END

  ViewVC Help
Powered by ViewVC 1.1.22