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