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