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

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

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

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