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 @ 217

Last change on this file since 217 was 217, checked in by opalod, 19 years ago
CT : UPDATE152 : optimization: reduce the MPI communications in doing only one CALL lbc_lnk() with a 3D array instead of jpk CALL lbc_lnk() with a 2D array; or suppressing useless CALL lbc_lnk()
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 35.2 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	! local Richardson number
381	zri = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )
382	! Prandtl number
383	zpdl = 1.0
384	IF( zri >= 0.2 ) zpdl = 0.2 / zri
385	zpdl = MAX( 0.1, zpdl )
386	! coefficient (esh2)
387	zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk)
388
389	! Matrix
390	zcof = zfact1 * tmask(ji,jj,jk)
391	! lower diagonal
392	avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk ) + zmxlm(ji,jj,jk-1) ) &
393	& / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )
394	! upper diagonal
395	avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk ) ) &
396	& / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) )
397	! diagonal
398	zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
399	! right hand side in en
400	en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en (ji,jj,jk) &
401	& + rdt * zmxlm(ji,jj,jk) * zesh2
402	! save masked Prandlt number in zmxlm array
403	zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)
404	END DO
405	END DO
406	END DO
407
408	END SELECT
409
410	! 4. Matrix inversion from level 2 (tke prescribed at level 1)
411	!!--------------------------------
412	! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
413	DO jk = 3, jpkm1
414	DO jj = 2, jpjm1
415	DO ji = fs_2, fs_jpim1 ! vector opt.
416	zwd(ji,jj,jk) = zwd(ji,jj,jk) - avmv(ji,jj,jk) * avmu(ji,jj,jk-1) / zwd(ji,jj,jk-1)
417	END DO
418	END DO
419	END DO
420
421	! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
422	DO jj = 2, jpjm1
423	DO ji = fs_2, fs_jpim1 ! vector opt.
424	avmv(ji,jj,2) = en(ji,jj,2) - avmv(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke
425	END DO
426	END DO
427	DO jk = 3, jpkm1
428	DO jj = 2, jpjm1
429	DO ji = fs_2, fs_jpim1 ! vector opt.
430	avmv(ji,jj,jk) = en(ji,jj,jk) - avmv(ji,jj,jk) / zwd(ji,jj,jk-1) *avmv(ji,jj,jk-1)
431	END DO
432	END DO
433	END DO
434
435	! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
436	DO jj = 2, jpjm1
437	DO ji = fs_2, fs_jpim1 ! vector opt.
438	en(ji,jj,jpkm1) = avmv(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
439	END DO
440	END DO
441	DO jk = jpk-2, 2, -1
442	DO jj = 2, jpjm1
443	DO ji = fs_2, fs_jpim1 ! vector opt.
444	en(ji,jj,jk) = ( avmv(ji,jj,jk) - avmu(ji,jj,jk) * en(ji,jj,jk+1) ) / zwd(ji,jj,jk)
445	END DO
446	END DO
447	END DO
448
449	! Save the result in en and set minimum value of tke : emin
450	DO jk = 2, jpkm1
451	DO jj = 2, jpjm1
452	DO ji = fs_2, fs_jpim1 ! vector opt.
453	en(ji,jj,jk) = MAX( en(ji,jj,jk), emin ) * tmask(ji,jj,jk)
454	END DO
455	END DO
456	END DO
457
458	! Lateral boundary conditions on ( avt, en ) (sign unchanged)
459	CALL lbc_lnk( en , 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. )
460
461
462	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
463	! III. Before vertical eddy vicosity and diffusivity coefficients
464	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
465
466	SELECT CASE ( nave )
467
468	CASE ( 0 ) ! no horizontal average
469
470	! Vertical eddy viscosity
471
472	DO jk = 2, jpkm1 ! Horizontal slab
473	DO jj = 2, jpjm1
474	DO ji = fs_2, fs_jpim1 ! vector opt.
475	avmu(ji,jj,jk) = ( avt (ji,jj,jk) + avt (ji+1,jj ,jk) ) * umask(ji,jj,jk) &
476	& / MAX( 1., tmask(ji,jj,jk) + tmask(ji+1,jj ,jk) )
477	avmv(ji,jj,jk) = ( avt (ji,jj,jk) + avt (ji ,jj+1,jk) ) * vmask(ji,jj,jk) &
478	& / MAX( 1., tmask(ji,jj,jk) + tmask(ji ,jj+1,jk) )
479	END DO
480	END DO
481	END DO
482
483	! Lateral boundary conditions (avmu,avmv) (U- and V- points, sign unchanged)
484	CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, '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 jk = 2, jpkm1 ! Horizontal slab
496	DO jj = 2, jpjm1
497	DO ji = fs_2, fs_jpim1 ! vector opt.
498	avmu(ji,jj,jk) = ( avt(ji,jj ,jk) + avt(ji+1,jj ,jk) &
499	& +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk) &
500	& +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk)
501
502	avmv(ji,jj,jk) = ( avt(ji ,jj,jk) + avt(ji ,jj+1,jk) &
503	& +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk) &
504	& +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk)
505	END DO
506	END DO
507	END DO
508	# else
509	DO jk = 2, jpkm1 ! Horizontal slab
510	DO jj = 2, jpjm1
511	DO ji = fs_2, fs_jpim1 ! vector opt.
512	avmu(ji,jj,jk) = ( avt (ji,jj ,jk) + avt (ji+1,jj ,jk) &
513	& +.5*( avt (ji,jj-1,jk) + avt (ji+1,jj-1,jk) &
514	& +avt (ji,jj+1,jk) + avt (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk) &
515	& / MAX( 1., tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) &
516	& +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk) &
517	& +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) ) )
518
519	avmv(ji,jj,jk) = ( avt (ji ,jj,jk) + avt (ji ,jj+1,jk) &
520	& +.5*( avt (ji-1,jj,jk) + avt (ji-1,jj+1,jk) &
521	& +avt (ji+1,jj,jk) + avt (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk) &
522	& / MAX( 1., tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
523	& +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk) &
524	& +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) ) )
525	END DO
526	END DO
527	END DO
528	# endif
529
530	! Lateral boundary conditions (avmu,avmv) (sign unchanged)
531	CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. )
532
533	! Vertical eddy diffusivity
534	! ------------------------------
535	! (1 2 1)
536	! horizontal average avt = 1/16 (2 4 2)
537	! (1 2 1)
538	DO jk = 2, jpkm1 ! Horizontal slab
539	# if defined key_vectopt_memory
540	DO jj = 2, jpjm1
541	DO ji = fs_2, fs_jpim1 ! vector opt.
542	avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj ,jk) &
543	& + avmv(ji,jj,jk) + avmv(ji ,jj-1,jk) ) * etmean(ji,jj,jk)
544	END DO
545	END DO
546	# else
547	DO jj = 2, jpjm1
548	DO ji = fs_2, fs_jpim1 ! vector opt.
549	avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj ,jk) &
550	& + avmv (ji,jj,jk) + avmv (ji ,jj-1,jk) ) * tmask(ji,jj,jk) &
551	& / MAX( 1., umask(ji,jj,jk) + umask(ji-1,jj ,jk) &
552	& + vmask(ji,jj,jk) + vmask(ji ,jj-1,jk) )
553	END DO
554	END DO
555	# endif
556	END DO
557
558	END SELECT
559
560	! multiplied by the Prandtl number (npdl>1)
561	! ----------------------------------------
562	IF( npdl == 1 ) THEN
563	DO jk = 2, jpkm1 ! Horizontal slab
564	DO jj = 2, jpjm1
565	DO ji = fs_2, fs_jpim1 ! vector opt.
566	zpdl = zmxld(ji,jj,jk)
567	avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
568	END DO
569	END DO
570	END DO
571	ENDIF
572
573	! Minimum value on the eddy viscosity
574	! ----------------------------------------
575	DO jk = 2, jpkm1 ! Horizontal slab
576	DO jj = 1, jpj
577	DO ji = 1, jpi
578	avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), avmb(jk) ) * umask(ji,jj,jk)
579	avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), avmb(jk) ) * vmask(ji,jj,jk)
580	END DO
581	END DO
582	END DO
583
584	! Lateral boundary conditions on avt (sign unchanged)
585	! ------------------------------=====
586	CALL lbc_lnk( avt, 'W', 1. )
587
588	IF(l_ctl) THEN
589	WRITE(numout,*) ' tke e : ', SUM( en (1:nictl+1,1:njctl+1,:) ), ' t : ', SUM( avt (1:nictl+1,1:njctl+1,:) )
590	WRITE(numout,*) ' u : ', SUM( avmu(1:nictl+1,1:njctl+1,:) ), ' v : ', SUM( avmv(1:nictl+1,1:njctl+1,:) )
591	ENDIF
592
593	END SUBROUTINE zdf_tke
594
595	# endif
596
597	SUBROUTINE zdf_tke_init
598	!!----------------------------------------------------------------------
599	!! * ROUTINE zdf_tke_init *
600	!!
601	!! ** Purpose : Initialization of the vertical eddy diffivity and
602	!! viscosity when using a tke turbulent closure scheme
603	!!
604	!! ** Method : Read the namtke namelist and check the parameters
605	!! called at the first timestep (nit000)
606	!!
607	!! ** input : Namlist namtke
608	!!
609	!! ** Action : Increase by 1 the nstop flag is setting problem encounter
610	!!
611	!! history :
612	!! 8.5 ! 02-06 (G. Madec) original code
613	!!----------------------------------------------------------------------
614	!! * Module used
615	USE dynzdf_exp
616	USE trazdf_exp
617
618	!! * local declarations
619	!! caution vectopt_memory change the solution (last digit of the solver stat)
620	# if defined key_vectopt_memory
621	INTEGER :: ji, jj, jk, jit ! dummy loop indices
622	# else
623	INTEGER :: jk, jit ! dummy loop indices
624	# endif
625
626	NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0, &
627	ri_c, nitke, nmxl, npdl, nave, navb
628	!!----------------------------------------------------------------------
629
630	! Read Namelist namtke : Turbulente Kinetic Energy
631	! --------------------
632	REWIND ( numnam )
633	READ ( numnam, namtke )
634
635	! Compute boost associated with the Richardson critic
636	! (control values: ri_c = 0.3 ==> eboost=1.25 for npdl=1 or 2)
637	! ( ri_c = 0.222 ==> eboost=1. )
638	eboost = ri_c * ( 2. + ediss / ediff ) / 2.
639
640
641	! Parameter control and print
642	! ---------------------------
643	! Control print
644	IF(lwp) THEN
645	WRITE(numout,*)
646	WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme'
647	WRITE(numout,*) '~~~~~~~~~~~~'
648	WRITE(numout,*) ' Namelist namtke : set tke mixing parameters'
649	WRITE(numout,*) ' restart with tke from no tke ln_rstke = ', ln_rstke
650	WRITE(numout,*) ' coef. to compute avt ediff = ', ediff
651	WRITE(numout,*) ' Kolmogoroff dissipation coef. ediss = ', ediss
652	WRITE(numout,*) ' tke surface input coef. ebb = ', ebb
653	WRITE(numout,*) ' tke diffusion coef. efave = ', efave
654	WRITE(numout,*) ' minimum value of tke emin = ', emin
655	WRITE(numout,*) ' surface minimum value of tke emin0 = ', emin0
656	WRITE(numout,*) ' number of restart iter loops nitke = ', nitke
657	WRITE(numout,*) ' mixing length type nmxl = ', nmxl
658	WRITE(numout,*) ' prandl number flag npdl = ', npdl
659	WRITE(numout,*) ' horizontal average flag nave = ', nave
660	WRITE(numout,*) ' critic Richardson nb ri_c = ', ri_c
661	WRITE(numout,*) ' and its associated coeff. eboost = ', eboost
662	WRITE(numout,*) ' constant background or profile navb = ', navb
663	WRITE(numout,*)
664	ENDIF
665
666	! Check nmxl and npdl values
667	IF( nmxl < 0 .OR. nmxl > 3 ) THEN
668	IF(lwp) WRITE(numout,cform_err)
669	IF(lwp) WRITE(numout,*) ' bad flag: nmxl is < 0 or > 3 '
670	nstop = nstop + 1
671	ENDIF
672	IF ( npdl < 0 .OR. npdl > 1 ) THEN
673	IF(lwp) WRITE(numout,cform_err)
674	IF(lwp) WRITE(numout,*) ' bad flag: npdl is < 0 or > 1 '
675	nstop = nstop + 1
676	ENDIF
677
678
679	! Horizontal average : initialization of weighting arrays
680	! -------------------
681
682	SELECT CASE ( nave )
683
684	CASE ( 0 ) ! no horizontal average
685	IF(lwp) WRITE(numout,*) ' no horizontal average on avt, avmu, avmv'
686	IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !'
687	!! caution vectopt_memory change the solution (last digit of the solver stat)
688	# if defined key_vectopt_memory
689	! weighting mean arrays etmean, eumean and evmean
690	! ( 1 1 ) ( 1 )
691	! avt = 1/4 ( 1 1 ) avmu = 1/2 ( 1 1 ) avmv= 1/2 ( 1 )
692	!
693	etmean(:,:,:) = 0.e0
694	eumean(:,:,:) = 0.e0
695	evmean(:,:,:) = 0.e0
696
697	DO jk = 1, jpkm1
698	DO jj = 2, jpjm1
699	DO ji = fs_2, fs_jpim1 ! vector opt.
700	etmean(ji,jj,jk) = tmask(ji,jj,jk) &
701	& / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &
702	& + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )
703
704	eumean(ji,jj,jk) = umask(ji,jj,jk) &
705	& / MAX( 1., tmask(ji,jj,jk) + tmask(ji+1,jj ,jk) )
706
707	evmean(ji,jj,jk) = vmask(ji,jj,jk) &
708	& / MAX( 1., tmask(ji,jj,jk) + tmask(ji ,jj+1,jk) )
709	END DO
710	END DO
711	END DO
712	# endif
713
714	CASE ( 1 ) ! horizontal average
715	IF(lwp) WRITE(numout,*) ' horizontal average on avt, avmu, avmv'
716	!! caution vectopt_memory change the solution (last digit of the solver stat)
717	# if defined key_vectopt_memory
718	! weighting mean arrays etmean, eumean and evmean
719	! ( 1 1 ) ( 1/2 1/2 ) ( 1/2 1 1/2 )
720	! avt = 1/4 ( 1 1 ) avmu = 1/4 ( 1 1 ) avmv= 1/4 ( 1/2 1 1/2 )
721	! ( 1/2 1/2 )
722	etmean(:,:,:) = 0.e0
723	eumean(:,:,:) = 0.e0
724	evmean(:,:,:) = 0.e0
725
726	DO jk = 1, jpkm1
727	DO jj = 2, jpjm1
728	DO ji = fs_2, fs_jpim1 ! vector opt.
729	etmean(ji,jj,jk) = tmask(ji,jj,jk) &
730	& / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &
731	& + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )
732
733	eumean(ji,jj,jk) = umask(ji,jj,jk) &
734	& / MAX( 1., tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) &
735	& +.5 * ( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk) &
736	& +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) ) )
737
738	evmean(ji,jj,jk) = vmask(ji,jj,jk) &
739	& / MAX( 1., tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
740	& +.5 * ( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk) &
741	& +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) ) )
742	END DO
743	END DO
744	END DO
745	# endif
746
747	CASE DEFAULT
748	IF(lwp) WRITE(numout,cform_err)
749	IF(lwp) WRITE(numout,*) ' bad flag value for nave = ', nave
750	nstop = nstop + 1
751
752	END SELECT
753
754
755	! Background eddy viscosity and diffusivity profil
756	! ------------------------------------------------
757	IF( navb == 0 ) THEN
758	! Define avmb, avtb from namelist parameter
759	avmb(:) = avm0
760	avtb(:) = avt0
761	ELSE
762	! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)
763	avmb(:) = avm0
764	avtb(:) = 1.e-5 + 2.8e-8 * gdepw(:) ! m2/s
765	ENDIF
766
767	! Increase the background in the surface layers
768	avmb(1) = 10. * avmb(1) ; avtb(1) = 10. * avtb(1)
769	avmb(2) = 10. * avmb(2) ; avtb(2) = 10. * avtb(2)
770	avmb(3) = 5. * avmb(3) ; avtb(3) = 5. * avtb(3)
771	avmb(4) = 2.5 * avmb(4) ; avtb(4) = 2.5 * avtb(4)
772
773
774	! Initialization of vertical eddy coef. to the background value
775	! -------------------------------------------------------------
776	DO jk = 1, jpk
777	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
778	avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
779	avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
780	END DO
781
782
783	! Initialization of turbulent kinetic energy ( en )
784	! -------------------------------------------------
785	IF( ln_rstart ) THEN
786	! no en field in the restart file, en set by iterative loop
787	IF( ln_rstke ) THEN
788	en (:,:,:) = emin * tmask(:,:,:)
789	DO jit = 2, nitke+1
790	CALL zdf_tke( jit )
791	END DO
792	ENDIF
793	! otherwise en is already read in dtrlec called by inidtr
794	ELSE
795	! no restart: en set to emin
796	en(:,:,:) = emin * tmask(:,:,:)
797	ENDIF
798
799	END SUBROUTINE zdf_tke_init
800
801	#else
802	!!----------------------------------------------------------------------
803	!! Dummy module : NO TKE scheme
804	!!----------------------------------------------------------------------
805	LOGICAL, PUBLIC, PARAMETER :: lk_zdftke = .FALSE. !: TKE flag
806	CONTAINS
807	SUBROUTINE zdf_tke( kt ) ! Empty routine
808	WRITE(,) 'zdf_tke: You should not have seen this print! error?', kt
809	END SUBROUTINE zdf_tke
810	#endif
811
812	!!======================================================================
813	END MODULE zdftke

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

Download in other formats: