/[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.24 - (hide annotations) (download)
Mon Nov 5 08:11:45 2012 UTC (11 years, 7 months ago) by mlosch
Branch: MAIN
Changes since 1.23: +22 -11 lines 
fix computation of control volumes for the case of
ALLOW_GGL90_HORIZDIFF (probably not used very often): move control
volumes up half a grid cell from C-centered to W-centeredupdating
(pointed out by C. Eden)
1 mlosch 1.24 C $Header: /u/gcmpack/MITgcm/pkg/ggl90/ggl90_calc.F,v 1.23 2012/08/08 22:22:42 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 jmc 1.21 I bi, bj, sigmaR, myTime, myIter, myThid )
12    
13 mlosch 1.1
14     C !DESCRIPTION: \bv
15 jmc 1.12 C *==========================================================*
16 mlosch 1.1 C | SUBROUTINE GGL90_CALC |
17     C | o Compute all GGL90 fields defined in GGL90.h |
18 jmc 1.12 C *==========================================================*
19 mlosch 1.1 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 jmc 1.12 C *==========================================================*
26 mlosch 1.1
27     C global parameters updated by ggl90_calc
28 jmc 1.12 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 mlosch 1.1 C \ev
32    
33     C !USES: ============================================================
34 jmc 1.12 IMPLICIT NONE
35 mlosch 1.1 #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 jmc 1.12 C Routine arguments
45 jmc 1.21 C bi, bj :: Current tile indices
46     C sigmaR :: Vertical gradient of iso-neutral density
47 jmc 1.12 C myTime :: Current time in simulation
48 jmc 1.19 C myIter :: Current time-step number
49 jmc 1.12 C myThid :: My Thread Id number
50 mlosch 1.1 INTEGER bi, bj
51 jmc 1.21 _RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
52 jmc 1.12 _RL myTime
53 jmc 1.19 INTEGER myIter
54 mlosch 1.1 INTEGER myThid
55 jmc 1.12 CEOP
56 mlosch 1.1
57     #ifdef ALLOW_GGL90
58    
59     C !LOCAL VARIABLES: ====================================================
60 jmc 1.12 C Local constants
61 jmc 1.19 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 mlosch 1.1 INTEGER iMin ,iMax ,jMin ,jMax
79 jmc 1.19 INTEGER i, j, k, kp1, km1, kSurf, kBottom
80     _RL explDissFac, implDissFac
81 mlosch 1.1 _RL uStarSquare
82     _RL verticalShear
83     _RL KappaM, KappaH
84 dfer 1.11 c _RL Nsquare
85     _RL Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86 mlosch 1.1 _RL deltaTggl90
87 dfer 1.11 c _RL SQRTTKE
88     _RL SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89 mlosch 1.1 _RL RiNumber
90     _RL TKEdissipation
91     _RL tempU, tempV, prTemp
92 gforget 1.16 _RL MaxLength, tmpmlx, tmpVisc
93 mlosch 1.1 _RL TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
94     _RL GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
95 dfer 1.13 _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 mlosch 1.1 _RL KappaE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
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 jmc 1.20 _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 dfer 1.15 INTEGER errCode
105 mlosch 1.2 #ifdef ALLOW_GGL90_HORIZDIFF
106 mlosch 1.24 C hFac :: fractional thickness of W-cell
107 jmc 1.19 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 mlosch 1.24 _RL hFac
111 mlosch 1.2 _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 jmc 1.19 _RL gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116 mlosch 1.2 #endif /* ALLOW_GGL90_HORIZDIFF */
117 dfer 1.11 #ifdef ALLOW_GGL90_SMOOTH
118 gforget 1.16 _RL p4, p8, p16
119 dfer 1.11 p4=0.25 _d 0
120     p8=0.125 _d 0
121     p16=0.0625 _d 0
122     #endif
123 mlosch 1.1 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 jmc 1.4 deltaTggl90 = dTtracerLev(1)
130 jmc 1.12
131 mlosch 1.1 kSurf = 1
132 jmc 1.19 C explicit/implicit timestepping weights for dissipation
133     explDissFac = 0. _d 0
134     implDissFac = 1. _d 0 - explDissFac
135 mlosch 1.1
136     C Initialize local fields
137 jmc 1.19 DO k = 1, Nr
138 jmc 1.21 DO j=1-OLy,sNy+OLy
139     DO i=1-OLx,sNx+OLx
140 jmc 1.19 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 jmc 1.20 #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 mlosch 1.1 ENDDO
150 jmc 1.8 ENDDO
151 mlosch 1.1 ENDDO
152 jmc 1.21 DO j=1-OLy,sNy+OLy
153     DO i=1-OLx,sNx+OLx
154 jmc 1.19 totalDepth(i,j) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
155 jmc 1.14 rMixingLength(i,j,1) = 0. _d 0
156 jmc 1.19 mxLength_Dn(i,j,1) = GGL90mixingLengthMin
157 jmc 1.14 SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
158 mlosch 1.1 ENDDO
159     ENDDO
160    
161     C start k-loop
162 jmc 1.19 DO k = 2, Nr
163 jmc 1.22 c km1 = k-1
164 jmc 1.19 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 jmc 1.14
169 mlosch 1.1 C buoyancy frequency
170 jmc 1.21 Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
171     & * sigmaR(i,j,k)
172 dfer 1.11 cC vertical shear term (dU/dz)^2+(dV/dz)^2
173 jmc 1.19 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 dfer 1.11 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 jmc 1.19 c TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
185 dfer 1.11 C mixing length
186 jmc 1.19 GGL90mixingLength(i,j,k) = SQRTTWO *
187 dfer 1.11 & SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
188     ENDDO
189     ENDDO
190     ENDDO
191    
192 gforget 1.23 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 dfer 1.13 C- Impose upper and lower bound for mixing length
202 dfer 1.11 IF ( mxlMaxFlag .EQ. 0 ) THEN
203 jmc 1.21
204 dfer 1.11 DO k=2,Nr
205 jmc 1.19 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 dfer 1.13 & MaxLength)
210 dfer 1.11 ENDDO
211     ENDDO
212     ENDDO
213 dfer 1.13
214     DO k=2,Nr
215 jmc 1.19 DO j=jMin,jMax
216     DO i=iMin,iMax
217     GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
218 dfer 1.13 & GGL90mixingLengthMin)
219 jmc 1.19 rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
220 dfer 1.13 ENDDO
221     ENDDO
222     ENDDO
223    
224 dfer 1.11 ELSEIF ( mxlMaxFlag .EQ. 1 ) 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=MIN(Ro_surf(i,j,bi,bj)-rF(k),rF(k)-R_low(i,j,bi,bj))
230 dfer 1.11 c MaxLength=MAX(MaxLength,20. _d 0)
231 jmc 1.19 GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
232 dfer 1.13 & MaxLength)
233     ENDDO
234     ENDDO
235     ENDDO
236    
237     DO k=2,Nr
238 jmc 1.19 DO j=jMin,jMax
239     DO i=iMin,iMax
240     GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
241 dfer 1.13 & GGL90mixingLengthMin)
242 jmc 1.19 rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
243 dfer 1.11 ENDDO
244     ENDDO
245     ENDDO
246 dfer 1.13
247 dfer 1.11 ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
248 jmc 1.21
249 dfer 1.11 DO k=2,Nr
250 jmc 1.19 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 dfer 1.11 ENDDO
255     ENDDO
256     ENDDO
257 jmc 1.19 DO j=jMin,jMax
258     DO i=iMin,iMax
259     GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
260 dfer 1.11 & GGL90mixingLengthMin+drF(Nr))
261     ENDDO
262     ENDDO
263     DO k=Nr-1,2,-1
264 jmc 1.19 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 jmc 1.12 ENDDO
269     ENDDO
270 dfer 1.11 ENDDO
271    
272 dfer 1.13 DO k=2,Nr
273 jmc 1.19 DO j=jMin,jMax
274     DO i=iMin,iMax
275     GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
276 dfer 1.13 & GGL90mixingLengthMin)
277 jmc 1.19 rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
278 dfer 1.13 ENDDO
279     ENDDO
280     ENDDO
281    
282     ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
283 jmc 1.21
284 dfer 1.13 DO k=2,Nr
285 jmc 1.19 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 dfer 1.13 ENDDO
290     ENDDO
291     ENDDO
292 jmc 1.19 DO j=jMin,jMax
293     DO i=iMin,iMax
294     GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
295 dfer 1.13 & GGL90mixingLengthMin+drF(Nr))
296     ENDDO
297     ENDDO
298     DO k=Nr-1,2,-1
299 jmc 1.19 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 dfer 1.13 ENDDO
304 dfer 1.11 ENDDO
305     ENDDO
306 dfer 1.13
307     DO k=2,Nr
308 jmc 1.19 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 dfer 1.13 tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
314 jmc 1.19 rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
315 dfer 1.13 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 jmc 1.19 c DO j=jMin,jMax
326     c DO i=iMin,iMax
327     c GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
328 dfer 1.13 c & GGL90mixingLengthMin)
329 jmc 1.19 c rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
330 dfer 1.13 c ENDDO
331     c ENDDO
332     c ENDDO
333 dfer 1.11
334     DO k=2,Nr
335 jmc 1.19 km1 = k-1
336 dfer 1.11
337 jmc 1.19 #ifdef ALLOW_GGL90_HORIZDIFF
338 mlosch 1.24 IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
339 jmc 1.8 C horizontal diffusion of TKE (requires an exchange in
340     C do_fields_blocking_exchanges)
341 mlosch 1.2 C common factors
342 jmc 1.21 DO j=1-OLy,sNy+OLy
343     DO i=1-OLx,sNx+OLx
344 mlosch 1.24 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 mlosch 1.2 ENDDO
351 jmc 1.8 ENDDO
352 mlosch 1.2 C Compute diffusive fluxes
353     C ... across x-faces
354 jmc 1.21 DO j=1-OLy,sNy+OLy
355     dfx(1-OLx,j)=0. _d 0
356     DO i=1-OLx+1,sNx+OLx
357 mlosch 1.2 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 mlosch 1.24 #ifdef ISOTROPIC_COS_SCALING
361 mlosch 1.2 & *CosFacU(j,bi,bj)
362 mlosch 1.24 #endif /* ISOTROPIC_COS_SCALING */
363     & *maskW(i,j,k,bi,bj) ! paranoia setting
364 mlosch 1.2 ENDDO
365     ENDDO
366     C ... across y-faces
367 jmc 1.21 DO i=1-OLx,sNx+OLx
368     dfy(i,1-OLy)=0. _d 0
369 mlosch 1.2 ENDDO
370 jmc 1.21 DO j=1-OLy+1,sNy+OLy
371     DO i=1-OLx,sNx+OLx
372 mlosch 1.2 dfy(i,j) = -GGL90diffTKEh*yA(i,j)
373     & *_recip_dyC(i,j,bi,bj)
374     & *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
375     #ifdef ISOTROPIC_COS_SCALING
376     & *CosFacV(j,bi,bj)
377     #endif /* ISOTROPIC_COS_SCALING */
378 mlosch 1.24 & *maskS(i,j,k,bi,bj) ! paranoia setting
379 mlosch 1.2 ENDDO
380 jmc 1.8 ENDDO
381 mlosch 1.2 C Compute divergence of fluxes
382 jmc 1.21 DO j=1-OLy,sNy+OLy-1
383     DO i=1-OLx,sNx+OLx-1
384 mlosch 1.24 hFac = min(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +
385     & min(.5 _d 0,_hFacC(i,j,k ,bi,bj) )
386     gTKE(i,j) = 0.0
387     if ( hFac .ne. 0.0 )
388     & gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac
389     & *((dfx(i+1,j)-dfx(i,j))
390     & +(dfy(i,j+1)-dfy(i,j)) )
391     & *maskC(i,j,k,bi,bj) ! paranoia setting
392 jmc 1.8 ENDDO
393 mlosch 1.2 ENDDO
394 jmc 1.19 C end if GGL90diffTKEh .eq. 0.
395     ENDIF
396     #endif /* ALLOW_GGL90_HORIZDIFF */
397    
398     DO j=jMin,jMax
399     DO i=iMin,iMax
400     C vertical shear term (dU/dz)^2+(dV/dz)^2
401     tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
402     & -( uVel(i,j,k ,bi,bj)+uVel(i+1,j,k ,bi,bj)) )
403     & *recip_drC(k)
404     tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
405     & -( vVel(i,j,k ,bi,bj)+vVel(i,j+1,k ,bi,bj)) )
406     & *recip_drC(k)
407     verticalShear = tempU*tempU + tempV*tempV
408     RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
409     C compute Prandtl number (always greater than 0)
410     prTemp = 1. _d 0
411     IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
412     TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
413     c TKEPrandtlNumber(i,j,k) = 1. _d 0
414    
415     C viscosity and diffusivity
416     KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
417     GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
418     & * maskC(i,j,k,bi,bj)
419     c note: storing GGL90visctmp like this, and using it later to compute
420     c GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
421     KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
422     KappaH = KappaM/TKEPrandtlNumber(i,j,k)
423     KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
424    
425     C dissipation term
426     TKEdissipation = explDissFac*GGL90ceps
427     & *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
428     & *GGL90TKE(i,j,k,bi,bj)
429     C partial update with sum of explicit contributions
430     GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
431     & + deltaTggl90*(
432     & + KappaM*verticalShear
433     & - KappaH*Nsquare(i,j,k)
434     & - TKEdissipation
435     & )
436     ENDDO
437 mlosch 1.2 ENDDO
438 jmc 1.19
439     #ifdef ALLOW_GGL90_HORIZDIFF
440     IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
441     C-- Add horiz. diffusion tendency
442     DO j=jMin,jMax
443     DO i=iMin,iMax
444     GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
445     & + gTKE(i,j)*deltaTggl90
446     ENDDO
447     ENDDO
448     ENDIF
449 mlosch 1.2 #endif /* ALLOW_GGL90_HORIZDIFF */
450    
451 jmc 1.19 C-- end of k loop
452     ENDDO
453    
454 mlosch 1.2 C ============================================
455     C Implicit time step to update TKE for k=1,Nr;
456     C TKE(Nr+1)=0 by default
457     C ============================================
458 mlosch 1.1 C set up matrix
459 mlosch 1.2 C-- Lower diagonal
460 mlosch 1.1 DO j=jMin,jMax
461     DO i=iMin,iMax
462 jmc 1.20 a3d(i,j,1) = 0. _d 0
463 mlosch 1.1 ENDDO
464     ENDDO
465     DO k=2,Nr
466 dfer 1.15 km1=MAX(2,k-1)
467 mlosch 1.1 DO j=jMin,jMax
468     DO i=iMin,iMax
469 dfer 1.15 C- We keep recip_hFacC in the diffusive flux calculation,
470     C- but no hFacC in TKE volume control
471     C- No need for maskC(k-1) with recip_hFacC(k-1)
472 jmc 1.20 a3d(i,j,k) = -deltaTggl90
473 dfer 1.11 & *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
474 dfer 1.10 & *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
475 dfer 1.15 & *recip_drC(k)*maskC(i,j,k,bi,bj)
476 mlosch 1.1 ENDDO
477     ENDDO
478     ENDDO
479 mlosch 1.2 C-- Upper diagonal
480 mlosch 1.1 DO j=jMin,jMax
481     DO i=iMin,iMax
482 jmc 1.20 c3d(i,j,1) = 0. _d 0
483 mlosch 1.1 ENDDO
484     ENDDO
485 dfer 1.11 DO k=2,Nr
486 mlosch 1.1 DO j=jMin,jMax
487     DO i=iMin,iMax
488 dfer 1.15 kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
489     C- We keep recip_hFacC in the diffusive flux calculation,
490     C- but no hFacC in TKE volume control
491     C- No need for maskC(k) with recip_hFacC(k)
492 jmc 1.20 c3d(i,j,k) = -deltaTggl90
493 dfer 1.11 & *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
494     & *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
495 dfer 1.15 & *recip_drC(k)*maskC(i,j,k-1,bi,bj)
496 mlosch 1.1 ENDDO
497     ENDDO
498     ENDDO
499 mlosch 1.2 C-- Center diagonal
500 mlosch 1.1 DO k=1,Nr
501 dfer 1.15 km1 = MAX(k-1,1)
502 mlosch 1.1 DO j=jMin,jMax
503     DO i=iMin,iMax
504 jmc 1.20 b3d(i,j,k) = 1. _d 0 - c3d(i,j,k) - a3d(i,j,k)
505 jmc 1.19 & + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
506     & * rMixingLength(i,j,k)
507 dfer 1.15 & * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
508 mlosch 1.1 ENDDO
509     ENDDO
510     ENDDO
511     C end set up matrix
512    
513     C Apply boundary condition
514 dfer 1.15 kp1 = MIN(Nr,kSurf+1)
515 jmc 1.19 DO j=jMin,jMax
516     DO i=iMin,iMax
517 mlosch 1.1 C estimate friction velocity uStar from surface forcing
518 jmc 1.8 uStarSquare = SQRT(
519 jmc 1.19 & ( .5 _d 0*( surfaceForcingU(i, j, bi,bj)
520     & + surfaceForcingU(i+1,j, bi,bj) ) )**2
521     & + ( .5 _d 0*( surfaceForcingV(i, j, bi,bj)
522     & + surfaceForcingV(i, j+1,bi,bj) ) )**2
523 mlosch 1.1 & )
524     C Dirichlet surface boundary condition for TKE
525 jmc 1.19 GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
526     & *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
527     GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
528 jmc 1.20 & - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
529     a3d(i,j,kp1) = 0. _d 0
530 mlosch 1.1 C Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
531 jmc 1.19 kBottom = MAX(kLowC(i,j,bi,bj),1)
532     GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
533 jmc 1.20 & - GGL90TKEbottom*c3d(i,j,kBottom)
534     c3d(i,j,kBottom) = 0. _d 0
535 mlosch 1.1 ENDDO
536 jmc 1.8 ENDDO
537 jmc 1.14
538 jmc 1.19 C solve tri-diagonal system
539 dfer 1.15 CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
540 jmc 1.20 I a3d, b3d, c3d,
541 jmc 1.19 U GGL90TKE,
542 dfer 1.15 O errCode,
543 jmc 1.19 I bi, bj, myThid )
544 jmc 1.14
545 jmc 1.19 DO k=1,Nr
546     DO j=jMin,jMax
547     DO i=iMin,iMax
548 mlosch 1.1 C impose minimum TKE to avoid numerical undershoots below zero
549 jmc 1.19 GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
550     & *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
551 mlosch 1.1 ENDDO
552     ENDDO
553 jmc 1.8 ENDDO
554 dfer 1.11
555 mlosch 1.2 C end of time step
556     C ===============================
557 dfer 1.11
558 jmc 1.19 DO k=2,Nr
559     DO j=1,sNy
560     DO i=1,sNx
561 gforget 1.16 #ifdef ALLOW_GGL90_SMOOTH
562     tmpVisc=
563 dfer 1.11 & (
564 gforget 1.16 & p4 * GGL90visctmp(i ,j ,k) * mskCor(i ,j ,bi,bj)
565     & +p8 *( GGL90visctmp(i-1,j ,k) * mskCor(i-1,j ,bi,bj)
566     & + GGL90visctmp(i ,j-1,k) * mskCor(i ,j-1,bi,bj)
567     & + GGL90visctmp(i+1,j ,k) * mskCor(i+1,j ,bi,bj)
568     & + GGL90visctmp(i ,j+1,k) * mskCor(i ,j+1,bi,bj))
569     & +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
570     & + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
571     & + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
572     & + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
573 dfer 1.11 & )
574     & /(p4
575     & +p8 *( maskC(i-1,j ,k,bi,bj) * mskCor(i-1,j ,bi,bj)
576     & + maskC(i ,j-1,k,bi,bj) * mskCor(i ,j-1,bi,bj)
577     & + maskC(i+1,j ,k,bi,bj) * mskCor(i+1,j ,bi,bj)
578     & + maskC(i ,j+1,k,bi,bj) * mskCor(i ,j+1,bi,bj))
579     & +p16*( maskC(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)
580     & + maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)
581     & + maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
582     & + maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))
583     & )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
584 gforget 1.16 #else
585 jmc 1.19 tmpVisc = GGL90visctmp(i,j,k)
586 gforget 1.16 #endif
587     tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
588 jmc 1.19 GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
589 dfer 1.11 ENDDO
590     ENDDO
591     ENDDO
592 gforget 1.16
593 jmc 1.19 DO k=2,Nr
594     DO j=1,sNy
595     DO i=1,sNx+1
596 gforget 1.16 #ifdef ALLOW_GGL90_SMOOTH
597 jmc 1.19 tmpVisc =
598 gforget 1.16 & (
599     & p4 *(GGL90visctmp(i ,j ,k) * mskCor(i ,j ,bi,bj)
600     & +GGL90visctmp(i-1,j ,k) * mskCor(i-1,j ,bi,bj))
601     & +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
602     & +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
603     & +GGL90visctmp(i ,j-1,k) * mskCor(i ,j-1,bi,bj)
604     & +GGL90visctmp(i ,j+1,k) * mskCor(i ,j+1,bi,bj))
605     & )
606     & /(p4 * 2. _d 0
607     & +p8 *( maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
608     & + maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
609     & + maskC(i ,j-1,k,bi,bj) * mskCor(i ,j-1,bi,bj)
610     & + maskC(i ,j+1,k,bi,bj) * mskCor(i ,j+1,bi,bj))
611     & )
612     & *maskC(i ,j,k,bi,bj)*mskCor(i ,j,bi,bj)
613     & *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
614     #else
615     tmpVisc = _maskW(i,j,k,bi,bj) *
616     & (.5 _d 0*(GGL90visctmp(i,j,k)
617     & +GGL90visctmp(i-1,j,k))
618     & )
619 dfer 1.11 #endif
620 jmc 1.19 tmpVisc = MIN( tmpVisc , GGL90viscMax )
621     GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
622 gforget 1.16 ENDDO
623     ENDDO
624     ENDDO
625    
626 jmc 1.19 DO k=2,Nr
627     DO j=1,sNy+1
628     DO i=1,sNx
629 gforget 1.16 #ifdef ALLOW_GGL90_SMOOTH
630     tmpVisc =
631     & (
632     & p4 *(GGL90visctmp(i ,j ,k) * mskCor(i ,j ,bi,bj)
633     & +GGL90visctmp(i ,j-1,k) * mskCor(i ,j-1,bi,bj))
634     & +p8 *(GGL90visctmp(i-1,j ,k) * mskCor(i-1,j ,bi,bj)
635     & +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
636     & +GGL90visctmp(i+1,j ,k) * mskCor(i+1,j ,bi,bj)
637     & +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))
638     & )
639     & /(p4 * 2. _d 0
640     & +p8 *( maskC(i-1,j ,k,bi,bj) * mskCor(i-1,j ,bi,bj)
641     & + maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
642     & + maskC(i+1,j ,k,bi,bj) * mskCor(i+1,j ,bi,bj)
643     & + maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))
644     & )
645     & *maskC(i,j ,k,bi,bj)*mskCor(i,j ,bi,bj)
646     & *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
647     #else
648     tmpVisc = _maskS(i,j,k,bi,bj) *
649     & (.5 _d 0*(GGL90visctmp(i,j,k)
650     & +GGL90visctmp(i,j-1,k))
651     & )
652    
653     #endif
654     tmpVisc = MIN( tmpVisc , GGL90viscMax )
655 jmc 1.19 GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
656 gforget 1.16 ENDDO
657     ENDDO
658     ENDDO
659 dfer 1.10
660     #ifdef ALLOW_DIAGNOSTICS
661     IF ( useDiagnostics ) THEN
662     CALL DIAGNOSTICS_FILL( GGL90TKE ,'GGL90TKE',
663     & 0,Nr, 1, bi, bj, myThid )
664 gforget 1.16 CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
665     & 0,Nr, 1, bi, bj, myThid )
666     CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
667 dfer 1.10 & 0,Nr, 1, bi, bj, myThid )
668     CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
669     & 0,Nr, 1, bi, bj, myThid )
670     CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
671     & 0,Nr, 2, bi, bj, myThid )
672     CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
673     & 0,Nr, 2, bi, bj, myThid )
674     ENDIF
675     #endif
676 mlosch 1.1
677     #endif /* ALLOW_GGL90 */
678    
679     RETURN
680     END
  ViewVC Help
Powered by ViewVC 1.1.22