New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

zdftke.F90 in trunk/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 96

Last change on this file since 96 was 49, checked in by opalod, 20 years ago
CT : UPDATE023 : Addition of new diagnostics controled with logical key l_ctl
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 35.1 KB

Line
1	MODULE zdftke
2	!!======================================================================
3	!! * MODULE zdftke *
4	!! Ocean physics: vertical mixing coefficient compute from the tke
5	!! turbulent closure parameterization
6	!!=====================================================================
7	#if defined key_zdftke \|\| defined key_esopa
8	!!----------------------------------------------------------------------
9	!! 'key_zdftke' TKE scheme
10	!!----------------------------------------------------------------------
11	!! zdf_tke : update momentum and tracer Kz from a tke scheme
12	!! zdf_tke_init : initialization, namelist read, and parameters control
13	!!----------------------------------------------------------------------
14	!! * Modules used
15	USE oce ! ocean dynamics and active tracers
16	USE dom_oce ! ocean space and time domain
17	USE zdf_oce ! ocean vertical physics
18	USE in_out_manager ! I/O manager
19	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
20	USE phycst ! physical constants
21	USE taumod ! surface stress
22
23	IMPLICIT NONE
24	PRIVATE
25
26	!! * Routine accessibility
27	PUBLIC zdf_tke ! routine called by step.F90
28
29	!! * Share Module variables
30	LOGICAL, PUBLIC, PARAMETER :: lk_zdftke = .TRUE. !: TKE vertical mixing flag
31	LOGICAL, PUBLIC :: & !!: tke namelist (namtke)
32	ln_rstke = .FALSE. !: =T restart with tke from a run without tke with
33	! ! a none zero initial value for en
34	REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & !:
35	en !: now turbulent kinetic energy
36
37	!! * Module variables
38	INTEGER :: & !!! tke namelist (namtke)
39	nitke = 50 , & ! number of restart iterative loops
40	nmxl = 2 , & ! = 0/1/2/3 flag for the type of mixing length used
41	npdl = 1 , & ! = 0/1/2 flag on prandtl number on vert. eddy coeff.
42	nave = 1 , & ! = 0/1 flag for horizontal average on avt, avmu, avmv
43	navb = 0 ! = 0/1 flag for constant or profile background avt
44	REAL(wp) :: & !!! tke namlist (namtke)
45	ediff = 0.1_wp , & ! coeff. for vertical eddy coef.; avt=ediffmxlsqrt(e)
46	ediss = 0.7_wp , & ! coef. of the Kolmogoroff dissipation
47	ebb = 3.75_wp , & ! coef. of the surface input of tke
48	efave = 1._wp , & ! coef. for the tke vert. diff. coeff.; avtke=efave*avm
49	emin = 0.7071e-6_wp , & ! minimum value of tke (m2/s2)
50	emin0 = 1.e-4_wp , & ! surface minimum value of tke (m2/s2)
51	ri_c = 2._wp / 9._wp ! critic Richardson number
52	REAL(wp) :: &
53	eboost ! multiplicative coeff of the shear product.
54
55	!! caution vectopt_memory change the solution (last digit of the solver stat)
56	# if defined key_vectopt_memory
57	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
58	etmean, & ! coefficient used for horizontal smoothing
59	eumean, & ! at t-, u- and v-points
60	evmean !
61	# endif
62
63	!! * Substitutions
64	# include "domzgr_substitute.h90"
65	# include "vectopt_loop_substitute.h90"
66	!!----------------------------------------------------------------------
67	!! OPA 9.0 , LODYC-IPSL (2003)
68	!!----------------------------------------------------------------------
69
70	CONTAINS
71
72	# if defined key_autotasking
73	!!----------------------------------------------------------------------
74	!! 'key_autotasking' : j-k-i loop (j-slab)
75	!!----------------------------------------------------------------------
76	# include "zdftke_atsk.h90"
77
78	# else
79	!!----------------------------------------------------------------------
80	!! default option : k-j-i loop
81	!!----------------------------------------------------------------------
82
83	SUBROUTINE zdf_tke ( kt )
84	!!----------------------------------------------------------------------
85	!! * ROUTINE zdf_tke *
86	!!
87	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
88	!! coefficients using a 1.5 turbulent closure scheme.
89	!!
90	!! ** Method : The time evolution of the turbulent kinetic energy
91	!! (tke) is computed from a prognostic equation :
92	!! d(en)/dt = eboost eav (d(u)/dz)**2 ! shear production
93	!! + d( efave eav d(en)/dz )/dz ! diffusion of tke
94	!! + grav/rau0 pdl eav d(rau)/dz ! stratif. destruc.
95	!! - ediss / emxl en**(2/3) ! dissipation
96	!! with the boundary conditions:
97	!! surface: en = max( emin0,ebb sqrt(taux^2 + tauy^2) )
98	!! bottom : en = emin
99	!! -1- The dissipation and mixing turbulent lengh scales are computed
100	!! from the usual diagnostic buoyancy length scale:
101	!! mxl= 1/(sqrt(en)/N) WHERE N is the brunt-vaisala frequency
102	!! Four cases :
103	!! nmxl=0 : mxl bounded by the distance to surface and bottom.
104	!! zmxld = zmxlm = mxl
105	!! nmxl=1 : mxl bounded by the vertical scale factor.
106	!! zmxld = zmxlm = mxl
107	!! nmxl=2 : mxl bounded such that the vertical derivative of mxl
108	!! is less than 1 (\|d/dz(xml)\|<1).
109	!! zmxld = zmxlm = mxl
110	!! nmxl=3 : lup = mxl bounded using \|d/dz(xml)\|<1 from the surface
111	!! to the bottom
112	!! ldown = mxl bounded using \|d/dz(xml)\|<1 from the bottom
113	!! to the surface
114	!! zmxld = sqrt (lup*ldown) ; zmxlm = min(lup,ldown)
115	!! -2- Compute the now Turbulent kinetic energy. The time differencing
116	!! is implicit for vertical diffusion term, linearized for kolmo-
117	!! goroff dissipation term, and explicit forward for both buoyancy
118	!! and dynamic production terms. Thus a tridiagonal linear system is
119	!! solved.
120	!! Note that - the shear production is multiplied by eboost in order
121	!! to set the critic richardson number to ri_c (namelist parameter)
122	!! - the destruction by stratification term is multiplied
123	!! by the Prandtl number (defined by an empirical funtion of the local
124	!! Richardson number) if npdl=1 (namelist parameter)
125	!! coefficient (zesh2):
126	!! -3- Compute the now vertical eddy vicosity and diffusivity
127	!! coefficients from en (before the time stepping) and zmxlm:
128	!! avm = max( avtb, ediffzmxlmen^1/2 )
129	!! avt = max( avmb, pdl*avm ) (pdl=1 if npdl=0)
130	!! eav = max( avmb, avm )
131	!! avt and avm are horizontally averaged to avoid numerical insta-
132	!! bilities.
133	!! N.B. The computation is done from jk=2 to jpkm1 except for
134	!! en. Surface value of avt avmu avmv are set once a time to
135	!! their background value in routine zdf_tke_init.
136	!!
137	!! ** Action : compute en (now turbulent kinetic energy)
138	!! update avt, avmu, avmv (before vertical eddy coef.)
139	!!
140	!! References :
141	!! Gaspar et al., jgr, 95, 1990,
142	!! Blanke and Delecluse, jpo, 1991
143	!! History :
144	!! 6.0 ! 91-03 (b. blanke) Original code
145	!! 7.0 ! 91-11 (G. Madec) bug fix
146	!! 7.1 ! 92-10 (G. Madec) new mixing length and eav
147	!! 7.2 ! 93-03 (M. Guyon) symetrical conditions
148	!! 7.3 ! 94-08 (G. Madec, M. Imbard) npdl flag
149	!! 7.5 ! 96-01 (G. Madec) s-coordinates
150	!! 8.0 ! 97-07 (G. Madec) lbc
151	!! 8.1 ! 99-01 (E. Stretta) new option for the mixing length
152	!! 8.5 ! 02-08 (G. Madec) ri_c and Free form, F90
153	!!----------------------------------------------------------------------
154	!! * Modules used
155	USE oce , zwd => ua, & ! use ua as workspace
156	& zmxlm => ta, & ! use ta as workspace
157	& zmxld => sa ! use sa as workspace
158
159	!! * arguments
160	INTEGER, INTENT( in ) :: kt ! ocean time step
161
162	!! * local declarations
163	INTEGER :: ji, jj, jk ! dummy loop arguments
164	REAL(wp) :: &
165	zmlmin, zbbrau, & ! temporary scalars
166	zfact1, zfact2, zfact3, & !
167	zrn2, zesurf, & !
168	ztx2, zty2, zav, & !
169	zcoef, zcof, zsh2, & !
170	zdku, zdkv, zpdl, zri, & !
171	zsqen, zesh2, & !
172	zemxl, zemlm, zemlp
173	!!--------------------------------------------------------------------
174
175	! Initialization (first time-step only)
176	! --------------
177	IF( kt == nit000 ) CALL zdf_tke_init
178
179	! Local constant initialization
180	zmlmin = 1.e-8
181	zbbrau = .5 * ebb / rau0
182	zfact1 = -.5 * rdt * efave
183	zfact2 = 1.5 * rdt * ediss
184	zfact3 = 0.5 * rdt * ediss
185
186
187	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
188	! I. Mixing length
189	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
190
191
192	! Buoyancy length scale: l=sqrt(2e/n*2)
193	! ---------------------
194	zmxlm(:,:, 1 ) = zmlmin ! surface set to the minimum value
195	zmxlm(:,:,jpk) = zmlmin ! bottom set to the minimum value
196	!CDIR NOVERRCHK
197	DO jk = 2, jpkm1
198	!CDIR NOVERRCHK
199	DO jj = 2, jpjm1
200	!CDIR NOVERRCHK
201	DO ji = fs_2, fs_jpim1 ! vector opt.
202	zrn2 = MAX( rn2(ji,jj,jk), rsmall )
203	zmxlm(ji,jj,jk) = MAX( SQRT( 2. * en(ji,jj,jk) / zrn2 ), zmlmin )
204	END DO
205	END DO
206	END DO
207
208
209	! Physical limits for the mixing length
210	! -------------------------------------
211	zmxld(:,:, 1 ) = zmlmin ! surface set to the minimum value
212	zmxld(:,:,jpk) = zmlmin ! bottom set to the minimum value
213
214	SELECT CASE ( nmxl )
215
216	CASE ( 0 ) ! bounded by the distance to surface and bottom
217	DO jk = 2, jpkm1
218	DO jj = 2, jpjm1
219	DO ji = fs_2, fs_jpim1 ! vector opt.
220	zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk), &
221	& fsdepw(ji,jj,mbathy(ji,jj)) - fsdepw(ji,jj,jk) )
222	zmxlm(ji,jj,jk) = zemxl
223	zmxld(ji,jj,jk) = zemxl
224	END DO
225	END DO
226	END DO
227
228	CASE ( 1 ) ! bounded by the vertical scale factor
229	DO jk = 2, jpkm1
230	DO jj = 2, jpjm1
231	DO ji = fs_2, fs_jpim1 ! vector opt.
232	zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
233	zmxlm(ji,jj,jk) = zemxl
234	zmxld(ji,jj,jk) = zemxl
235	END DO
236	END DO
237	END DO
238
239	CASE ( 2 ) ! \|dk[xml]\| bounded by e3t :
240	DO jk = 2, jpkm1 ! from the surface to the bottom :
241	DO jj = 2, jpjm1
242	DO ji = fs_2, fs_jpim1 ! vector opt.
243	zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
244	END DO
245	END DO
246	END DO
247	DO jk = jpkm1, 2, -1 ! from the bottom to the surface :
248	DO jj = 2, jpjm1
249	DO ji = fs_2, fs_jpim1 ! vector opt.
250	zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
251	zmxlm(ji,jj,jk) = zemxl
252	zmxld(ji,jj,jk) = zemxl
253	END DO
254	END DO
255	END DO
256
257	CASE ( 3 ) ! lup and ldown, \|dk[xml]\| bounded by e3t :
258	DO jk = 2, jpkm1 ! from the surface to the bottom : lup
259	DO jj = 2, jpjm1
260	DO ji = fs_2, fs_jpim1 ! vector opt.
261	zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
262	END DO
263	END DO
264	END DO
265	DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown
266	DO jj = 2, jpjm1
267	DO ji = fs_2, fs_jpim1 ! vector opt.
268	zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
269	END DO
270	END DO
271	END DO
272	!CDIR NOVERRCHK
273	DO jk = 2, jpkm1
274	!CDIR NOVERRCHK
275	DO jj = 2, jpjm1
276	!CDIR NOVERRCHK
277	DO ji = fs_2, fs_jpim1 ! vector opt.
278	zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) )
279	zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
280	zmxlm(ji,jj,jk) = zemlm
281	zmxld(ji,jj,jk) = zemlp
282	END DO
283	END DO
284	END DO
285
286	END SELECT
287
288
289	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
290	! II Tubulent kinetic energy time stepping
291	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
292
293
294	! 1. Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt
295	! ---------------------------------------------------------------------
296	!CDIR NOVERRCHK
297	DO jk = 2, jpkm1
298	!CDIR NOVERRCHK
299	DO jj = 2, jpjm1
300	!CDIR NOVERRCHK
301	DO ji = fs_2, fs_jpim1 ! vector opt.
302	zsqen = SQRT( en(ji,jj,jk) )
303	zav = ediff * zmxlm(ji,jj,jk) * zsqen
304	avt (ji,jj,jk) = MAX( zav, avtb(jk) ) * tmask(ji,jj,jk)
305	zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)
306	zmxld(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
307	END DO
308	END DO
309	END DO
310
311	! 2. Surface boundary condition on tke and its eddy viscosity (zmxlm)
312	! -------------------------------------------------
313	! en(1) = ebb sqrt(taux^2+tauy^2) / rau0 (min value emin0)
314	! zmxlm(1) = avmb(1) and zmxlm(jpk) = 0.
315	!CDIR NOVERRCHK
316	DO jj = 2, jpjm1
317	!CDIR NOVERRCHK
318	DO ji = fs_2, fs_jpim1 ! vector opt.
319	ztx2 = taux(ji-1,jj ) + taux(ji,jj)
320	zty2 = tauy(ji ,jj-1) + tauy(ji,jj)
321	zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
322	en (ji,jj,1) = MAX( zesurf, emin0 ) * tmask(ji,jj,1)
323	zmxlm(ji,jj,1 ) = avmb(1) * tmask(ji,jj,1)
324	zmxlm(ji,jj,jpk) = 0.e0
325	END DO
326	END DO
327
328	! 3. Now Turbulent kinetic energy (output in en)
329	! -------------------------------
330	! Resolution of a tridiagonal linear system by a "methode de chasse"
331	! computation from level 2 to jpkm1 (e(1) already computed and
332	! e(jpk)=0 ).
333
334	SELECT CASE ( npdl )
335
336	CASE ( 0 ) ! No Prandtl number
337	DO jk = 2, jpkm1
338	DO jj = 2, jpjm1
339	DO ji = fs_2, fs_jpim1 ! vector opt.
340	! zesh2 = eboost * (du/dz)^2 - N^2
341	zcoef = 0.5 / fse3w(ji,jj,jk)
342	! shear
343	zdku = zcoef * ( ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1) &
344	& - ub(ji-1, jj ,jk ) - ub(ji,jj,jk ) )
345	zdkv = zcoef * ( vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1) &
346	& - vb( ji ,jj-1,jk ) - vb(ji,jj,jk ) )
347	! coefficient (zesh2)
348	zesh2 = eboost * ( zdkuzdku + zdkvzdkv ) - rn2(ji,jj,jk)
349
350	! Matrix
351	zcof = zfact1 * tmask(ji,jj,jk)
352	! lower diagonal
353	avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk ) + zmxlm(ji,jj,jk-1) ) &
354	& / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )
355	! upper diagonal
356	avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk ) ) &
357	& / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) )
358	! diagonal
359	zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
360	! right hand side in en
361	en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en (ji,jj,jk) &
362	& + rdt * zmxlm(ji,jj,jk) * zesh2
363	END DO
364	END DO
365	END DO
366
367	CASE ( 1 ) ! Prandtl number
368	DO jk = 2, jpkm1
369	DO jj = 2, jpjm1
370	DO ji = fs_2, fs_jpim1 ! vector opt.
371	! zesh2 = eboost * (du/dz)^2 - pdl * N^2
372	zcoef = 0.5 / fse3w(ji,jj,jk)
373	! shear
374	zdku = zcoef * ( ub(ji-1,jj ,jk-1) + ub(ji,jj,jk-1) &
375	& - ub(ji-1,jj ,jk ) - ub(ji,jj,jk ) )
376	zdkv = zcoef * ( vb(ji ,jj-1,jk-1) + vb(ji,jj,jk-1) &
377	& - vb(ji ,jj-1,jk ) - vb(ji,jj,jk ) )
378	! square of vertical shear
379	zsh2 = zdku * zdku + zdkv * zdkv
380	! Prandtl number
381	zri = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )
382	zpdl = 1.0
383	IF( zri >= 0.2 ) zpdl = 0.2 / zri
384	zpdl = MAX( 0.1, zpdl )
385	! coefficient (esh2)
386	zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk)
387
388	! Matrix
389	zcof = zfact1 * tmask(ji,jj,jk)
390	! lower diagonal
391	avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk ) + zmxlm(ji,jj,jk-1) ) &
392	& / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )
393	! upper diagonal
394	avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk ) ) &
395	& / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) )
396	! diagonal
397	zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
398	! right hand side in en
399	en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en (ji,jj,jk) &
400	& + rdt * zmxlm(ji,jj,jk) * zesh2
401	! save masked Prandlt number in zmxlm array
402	zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)
403	END DO
404	END DO
405	END DO
406
407	END SELECT
408
409	! 4. Matrix inversion from level 2 (tke prescribed at level 1)
410	!!--------------------------------
411	! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
412	DO jk = 3, jpkm1
413	DO jj = 2, jpjm1
414	DO ji = fs_2, fs_jpim1 ! vector opt.
415	zwd(ji,jj,jk) = zwd(ji,jj,jk) - avmv(ji,jj,jk) * avmu(ji,jj,jk-1) / zwd(ji,jj,jk-1)
416	END DO
417	END DO
418	END DO
419
420	! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
421	DO jj = 2, jpjm1
422	DO ji = fs_2, fs_jpim1 ! vector opt.
423	avmv(ji,jj,2) = en(ji,jj,2) - avmv(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke
424	END DO
425	END DO
426	DO jk = 3, jpkm1
427	DO jj = 2, jpjm1
428	DO ji = fs_2, fs_jpim1 ! vector opt.
429	avmv(ji,jj,jk) = en(ji,jj,jk) - avmv(ji,jj,jk) / zwd(ji,jj,jk-1) *avmv(ji,jj,jk-1)
430	END DO
431	END DO
432	END DO
433
434	! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
435	DO jj = 2, jpjm1
436	DO ji = fs_2, fs_jpim1 ! vector opt.
437	en(ji,jj,jpkm1) = avmv(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
438	END DO
439	END DO
440	DO jk = jpk-2, 2, -1
441	DO jj = 2, jpjm1
442	DO ji = fs_2, fs_jpim1 ! vector opt.
443	en(ji,jj,jk) = ( avmv(ji,jj,jk) - avmu(ji,jj,jk) * en(ji,jj,jk+1) ) / zwd(ji,jj,jk)
444	END DO
445	END DO
446	END DO
447
448	! Save the result in en and set minimum value of tke : emin
449	DO jk = 2, jpkm1
450	DO jj = 2, jpjm1
451	DO ji = fs_2, fs_jpim1 ! vector opt.
452	en(ji,jj,jk) = MAX( en(ji,jj,jk), emin ) * tmask(ji,jj,jk)
453	END DO
454	END DO
455	END DO
456
457	! Lateral boundary conditions on ( avt, en ) (sign unchanged)
458	CALL lbc_lnk( en , 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. )
459
460
461	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
462	! III. Before vertical eddy vicosity and diffusivity coefficients
463	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
464
465	! ! ===============
466	DO jk = 2, jpkm1 ! Horizontal slab
467	! ! ===============
468	SELECT CASE ( nave )
469
470	CASE ( 0 ) ! no horizontal average
471
472	! Vertical eddy viscosity
473
474	DO jj = 2, jpjm1
475	DO ji = fs_2, fs_jpim1 ! vector opt.
476	avmu(ji,jj,jk) = ( avt (ji,jj,jk) + avt (ji+1,jj ,jk) ) * umask(ji,jj,jk) &
477	& / MAX( 1., tmask(ji,jj,jk) + tmask(ji+1,jj ,jk) )
478	avmv(ji,jj,jk) = ( avt (ji,jj,jk) + avt (ji ,jj+1,jk) ) * vmask(ji,jj,jk) &
479	& / MAX( 1., tmask(ji,jj,jk) + tmask(ji ,jj+1,jk) )
480	END DO
481	END DO
482
483	! Lateral boundary conditions (avmu,avmv) (U- and V- points, sign unchanged)
484	CALL lbc_lnk( avmu(:,:,jk), 'U', 1. ) ; CALL lbc_lnk( avmv(:,:,jk), 'V', 1. )
485
486	CASE ( 1 ) ! horizontal average
487
488	! ( 1/2 1/2 )
489	! Eddy viscosity: horizontal average: avmu = 1/4 ( 1 1 )
490	! ( 1/2 1 1/2 ) ( 1/2 1/2 )
491	! avmv = 1/4 ( 1/2 1 1/2 )
492
493	!! caution vectopt_memory change the solution (last digit of the solver stat)
494	# if defined key_vectopt_memory
495	DO jj = 2, jpjm1
496	DO ji = fs_2, fs_jpim1 ! vector opt.
497	avmu(ji,jj,jk) = ( avt(ji,jj ,jk) + avt(ji+1,jj ,jk) &
498	& +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk) &
499	& +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk)
500
501	avmv(ji,jj,jk) = ( avt(ji ,jj,jk) + avt(ji ,jj+1,jk) &
502	& +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk) &
503	& +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk)
504	END DO
505	END DO
506	# else
507	DO jj = 2, jpjm1
508	DO ji = fs_2, fs_jpim1 ! vector opt.
509	avmu(ji,jj,jk) = ( avt (ji,jj ,jk) + avt (ji+1,jj ,jk) &
510	& +.5*( avt (ji,jj-1,jk) + avt (ji+1,jj-1,jk) &
511	& +avt (ji,jj+1,jk) + avt (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk) &
512	& / MAX( 1., tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) &
513	& +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk) &
514	& +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) ) )
515
516	avmv(ji,jj,jk) = ( avt (ji ,jj,jk) + avt (ji ,jj+1,jk) &
517	& +.5*( avt (ji-1,jj,jk) + avt (ji-1,jj+1,jk) &
518	& +avt (ji+1,jj,jk) + avt (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk) &
519	& / MAX( 1., tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
520	& +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk) &
521	& +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) ) )
522	END DO
523	END DO
524	# endif
525
526	! Lateral boundary conditions (avmu,avmv) (sign unchanged)
527	CALL lbc_lnk( avmu(:,:,jk), 'U', 1. ) ; CALL lbc_lnk( avmv(:,:,jk), 'V', 1. )
528
529	! Vertical eddy diffusivity
530	! ------------------------------
531	! (1 2 1)
532	! horizontal average avt = 1/16 (2 4 2)
533	! (1 2 1)
534	# if defined key_vectopt_memory
535	DO jj = 2, jpjm1
536	DO ji = fs_2, fs_jpim1 ! vector opt.
537	avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj ,jk) &
538	& + avmv(ji,jj,jk) + avmv(ji ,jj-1,jk) ) * etmean(ji,jj,jk)
539	END DO
540	END DO
541	# else
542	DO jj = 2, jpjm1
543	DO ji = fs_2, fs_jpim1 ! vector opt.
544	avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj ,jk) &
545	& + avmv (ji,jj,jk) + avmv (ji ,jj-1,jk) ) * tmask(ji,jj,jk) &
546	& / MAX( 1., umask(ji,jj,jk) + umask(ji-1,jj ,jk) &
547	& + vmask(ji,jj,jk) + vmask(ji ,jj-1,jk) )
548	END DO
549	END DO
550	# endif
551	END SELECT
552
553	! multiplied by the Prandtl number (npdl>1)
554	! ----------------------------------------
555	IF( npdl == 1 ) THEN
556	DO jj = 2, jpjm1
557	DO ji = fs_2, fs_jpim1 ! vector opt.
558	zpdl = zmxld(ji,jj,jk)
559	avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
560	END DO
561	END DO
562	ENDIF
563
564	! Minimum value on the eddy viscosity
565	! ----------------------------------------
566	DO jj = 1, jpj
567	DO ji = 1, jpi
568	avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), avmb(jk) ) * umask(ji,jj,jk)
569	avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), avmb(jk) ) * vmask(ji,jj,jk)
570	END DO
571	END DO
572	! ! ===============
573	END DO ! End of slab
574	! ! ===============
575
576	! Lateral boundary conditions on avt (sign unchanged)
577	! ------------------------------=====
578	CALL lbc_lnk( avt, 'W', 1. )
579
580	IF( l_ctl ) THEN
581	WRITE(numout,*) ' tke e : ', SUM( en (:,:,:) ), ' t : ', SUM( avt (:,:,:) )
582	WRITE(numout,*) ' u : ', SUM( avmu(:,:,:) ), ' v : ', SUM( avmv(:,:,:) )
583	ENDIF
584
585	END SUBROUTINE zdf_tke
586
587	# endif
588
589	SUBROUTINE zdf_tke_init
590	!!----------------------------------------------------------------------
591	!! * ROUTINE zdf_tke_init *
592	!!
593	!! ** Purpose : Initialization of the vertical eddy diffivity and
594	!! viscosity when using a tke turbulent closure scheme
595	!!
596	!! ** Method : Read the namtke namelist and check the parameters
597	!! called at the first timestep (nit000)
598	!!
599	!! ** input : Namlist namtke
600	!!
601	!! ** Action : Increase by 1 the nstop flag is setting problem encounter
602	!!
603	!! history :
604	!! 8.5 ! 02-06 (G. Madec) original code
605	!!----------------------------------------------------------------------
606	!! * Module used
607	USE dynzdf_exp
608	USE trazdf_exp
609
610	!! * local declarations
611	!! caution vectopt_memory change the solution (last digit of the solver stat)
612	# if defined key_vectopt_memory
613	INTEGER :: ji, jj, jk, jit ! dummy loop indices
614	# else
615	INTEGER :: jk, jit ! dummy loop indices
616	# endif
617
618	NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0, &
619	ri_c, nitke, nmxl, npdl, nave, navb
620	!!----------------------------------------------------------------------
621
622	! Read Namelist namtke : Turbulente Kinetic Energy
623	! --------------------
624	REWIND ( numnam )
625	READ ( numnam, namtke )
626
627	! Compute boost associated with the Richardson critic
628	! (control values: ri_c = 0.3 ==> eboost=1.25 for npdl=1 or 2)
629	! ( ri_c = 0.222 ==> eboost=1. )
630	eboost = ri_c * ( 2. + ediss / ediff ) / 2.
631
632
633	! Parameter control and print
634	! ---------------------------
635	! Control print
636	IF(lwp) THEN
637	WRITE(numout,*)
638	WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme'
639	WRITE(numout,*) '~~~~~~~~~~~~'
640	WRITE(numout,*) ' Namelist namtke : set tke mixing parameters'
641	WRITE(numout,*) ' restart with tke from no tke ln_rstke = ', ln_rstke
642	WRITE(numout,*) ' coef. to compute avt ediff = ', ediff
643	WRITE(numout,*) ' Kolmogoroff dissipation coef. ediss = ', ediss
644	WRITE(numout,*) ' tke surface input coef. ebb = ', ebb
645	WRITE(numout,*) ' tke diffusion coef. efave = ', efave
646	WRITE(numout,*) ' minimum value of tke emin = ', emin
647	WRITE(numout,*) ' surface minimum value of tke emin0 = ', emin0
648	WRITE(numout,*) ' number of restart iter loops nitke = ', nitke
649	WRITE(numout,*) ' mixing length type nmxl = ', nmxl
650	WRITE(numout,*) ' prandl number flag npdl = ', npdl
651	WRITE(numout,*) ' horizontal average flag nave = ', nave
652	WRITE(numout,*) ' critic Richardson nb ri_c = ', ri_c
653	WRITE(numout,*) ' and its associated coeff. eboost = ', eboost
654	WRITE(numout,*) ' constant background or profile navb = ', navb
655	WRITE(numout,*)
656	ENDIF
657
658	! Check nmxl and npdl values
659	IF( nmxl < 0 .OR. nmxl > 3 ) THEN
660	IF(lwp) WRITE(numout,cform_err)
661	IF(lwp) WRITE(numout,*) ' bad flag: nmxl is < 0 or > 3 '
662	nstop = nstop + 1
663	ENDIF
664	IF ( npdl < 0 .OR. npdl > 1 ) THEN
665	IF(lwp) WRITE(numout,cform_err)
666	IF(lwp) WRITE(numout,*) ' bad flag: npdl is < 0 or > 1 '
667	nstop = nstop + 1
668	ENDIF
669
670
671	! Horizontal average : initialization of weighting arrays
672	! -------------------
673
674	SELECT CASE ( nave )
675
676	CASE ( 0 ) ! no horizontal average
677	IF(lwp) WRITE(numout,*) ' no horizontal average on avt, avmu, avmv'
678	IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !'
679	!! caution vectopt_memory change the solution (last digit of the solver stat)
680	# if defined key_vectopt_memory
681	! weighting mean arrays etmean, eumean and evmean
682	! ( 1 1 ) ( 1 )
683	! avt = 1/4 ( 1 1 ) avmu = 1/2 ( 1 1 ) avmv= 1/2 ( 1 )
684	!
685	etmean(:,:,:) = 0.e0
686	eumean(:,:,:) = 0.e0
687	evmean(:,:,:) = 0.e0
688
689	DO jk = 1, jpkm1
690	DO jj = 2, jpjm1
691	DO ji = fs_2, fs_jpim1 ! vector opt.
692	etmean(ji,jj,jk) = tmask(ji,jj,jk) &
693	& / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &
694	& + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )
695
696	eumean(ji,jj,jk) = umask(ji,jj,jk) &
697	& / MAX( 1., tmask(ji,jj,jk) + tmask(ji+1,jj ,jk) )
698
699	evmean(ji,jj,jk) = vmask(ji,jj,jk) &
700	& / MAX( 1., tmask(ji,jj,jk) + tmask(ji ,jj+1,jk) )
701	END DO
702	END DO
703	END DO
704	# endif
705
706	CASE ( 1 ) ! horizontal average
707	IF(lwp) WRITE(numout,*) ' horizontal average on avt, avmu, avmv'
708	!! caution vectopt_memory change the solution (last digit of the solver stat)
709	# if defined key_vectopt_memory
710	! weighting mean arrays etmean, eumean and evmean
711	! ( 1 1 ) ( 1/2 1/2 ) ( 1/2 1 1/2 )
712	! avt = 1/4 ( 1 1 ) avmu = 1/4 ( 1 1 ) avmv= 1/4 ( 1/2 1 1/2 )
713	! ( 1/2 1/2 )
714	etmean(:,:,:) = 0.e0
715	eumean(:,:,:) = 0.e0
716	evmean(:,:,:) = 0.e0
717
718	DO jk = 1, jpkm1
719	DO jj = 2, jpjm1
720	DO ji = fs_2, fs_jpim1 ! vector opt.
721	etmean(ji,jj,jk) = tmask(ji,jj,jk) &
722	& / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &
723	& + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )
724
725	eumean(ji,jj,jk) = umask(ji,jj,jk) &
726	& / MAX( 1., tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) &
727	& +.5 * ( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk) &
728	& +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) ) )
729
730	evmean(ji,jj,jk) = vmask(ji,jj,jk) &
731	& / MAX( 1., tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
732	& +.5 * ( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk) &
733	& +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) ) )
734	END DO
735	END DO
736	END DO
737	# endif
738
739	CASE DEFAULT
740	IF(lwp) WRITE(numout,cform_err)
741	IF(lwp) WRITE(numout,*) ' bad flag value for nave = ', nave
742	nstop = nstop + 1
743
744	END SELECT
745
746
747	! Background eddy viscosity and diffusivity profil
748	! ------------------------------------------------
749	IF( navb == 0 ) THEN
750	! Define avmb, avtb from namelist parameter
751	avmb(:) = avm0
752	avtb(:) = avt0
753	ELSE
754	! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)
755	avmb(:) = avm0
756	avtb(:) = 1.e-5 + 2.8e-8 * gdepw(:) ! m2/s
757	ENDIF
758
759	! Increase the background in the surface layers
760	avmb(1) = 10. * avmb(1) ; avtb(1) = 10. * avtb(1)
761	avmb(2) = 10. * avmb(2) ; avtb(2) = 10. * avtb(2)
762	avmb(3) = 5. * avmb(3) ; avtb(3) = 5. * avtb(3)
763	avmb(4) = 2.5 * avmb(4) ; avtb(4) = 2.5 * avtb(4)
764
765
766	! Initialization of vertical eddy coef. to the background value
767	! -------------------------------------------------------------
768	DO jk = 1, jpk
769	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
770	avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
771	avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
772	END DO
773
774
775	! Initialization of turbulent kinetic energy ( en )
776	! -------------------------------------------------
777	IF( ln_rstart ) THEN
778	! no en field in the restart file, en set by iterative loop
779	IF( ln_rstke ) THEN
780	en (:,:,:) = emin * tmask(:,:,:)
781	DO jit = 2, nitke+1
782	CALL zdf_tke( jit )
783	END DO
784	ENDIF
785	! otherwise en is already read in dtrlec called by inidtr
786	ELSE
787	! no restart: en set to emin
788	en(:,:,:) = emin * tmask(:,:,:)
789	ENDIF
790
791	END SUBROUTINE zdf_tke_init
792
793	#else
794	!!----------------------------------------------------------------------
795	!! Dummy module : NO TKE scheme
796	!!----------------------------------------------------------------------
797	LOGICAL, PUBLIC, PARAMETER :: lk_zdftke = .FALSE. !: TKE flag
798	CONTAINS
799	SUBROUTINE zdf_tke( kt ) ! Empty routine
800	WRITE(,) 'zdf_tke: You should not have seen this print! error?', kt
801	END SUBROUTINE zdf_tke
802	#endif
803
804	!!======================================================================
805	END MODULE zdftke

Note: See TracBrowser for help on using the repository browser.

Download in other formats: