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