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

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

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

Download in other formats: