Context Navigation

source: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/zdftke.F90 @ 3293

Last change on this file since 3293 was 3293, checked in by cetlod, 7 years ago
CM6.0.11 : add new NEMO routines & updated namelists and xml
File size: 49.9 KB

Line
1	MODULE zdftke
2	!!======================================================================
3	!! * MODULE zdftke *
4	!! Ocean physics: vertical mixing coefficient computed from the tke
5	!! turbulent closure parameterization
6	!!=====================================================================
7	!! History : OPA ! 1991-03 (b. blanke) Original code
8	!! 7.0 ! 1991-11 (G. Madec) bug fix
9	!! 7.1 ! 1992-10 (G. Madec) new mixing length and eav
10	!! 7.2 ! 1993-03 (M. Guyon) symetrical conditions
11	!! 7.3 ! 1994-08 (G. Madec, M. Imbard) nn_pdl flag
12	!! 7.5 ! 1996-01 (G. Madec) s-coordinates
13	!! 8.0 ! 1997-07 (G. Madec) lbc
14	!! 8.1 ! 1999-01 (E. Stretta) new option for the mixing length
15	!! NEMO 1.0 ! 2002-06 (G. Madec) add tke_init routine
16	!! - ! 2004-10 (C. Ethe ) 1D configuration
17	!! 2.0 ! 2006-07 (S. Masson) distributed restart using iom
18	!! 3.0 ! 2008-05 (C. Ethe, G.Madec) : update TKE physics:
19	!! ! - tke penetration (wind steering)
20	!! ! - suface condition for tke & mixing length
21	!! ! - Langmuir cells
22	!! - ! 2008-05 (J.-M. Molines, G. Madec) 2D form of avtb
23	!! - ! 2008-06 (G. Madec) style + DOCTOR name for namelist parameters
24	!! - ! 2008-12 (G. Reffray) stable discretization of the production term
25	!! 3.2 ! 2009-06 (G. Madec, S. Masson) TKE restart compatible with key_cpl
26	!! ! + cleaning of the parameters + bugs correction
27	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
28	!! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability
29	!!----------------------------------------------------------------------
30	#if defined key_zdftke \|\| defined key_esopa
31	!!----------------------------------------------------------------------
32	!! 'key_zdftke' TKE vertical physics
33	!!----------------------------------------------------------------------
34	!! zdf_tke : update momentum and tracer Kz from a tke scheme
35	!! tke_tke : tke time stepping: update tke at now time step (en)
36	!! tke_avn : compute mixing length scale and deduce avm and avt
37	!! zdf_tke_init : initialization, namelist read, and parameters control
38	!! tke_rst : read/write tke restart in ocean restart file
39	!!----------------------------------------------------------------------
40	USE oce ! ocean: dynamics and active tracers variables
41	USE phycst ! physical constants
42	USE dom_oce ! domain: ocean
43	USE domvvl ! domain: variable volume layer
44	USE sbc_oce ! surface boundary condition: ocean
45	USE zdf_oce ! vertical physics: ocean variables
46	USE zdfmxl ! vertical physics: mixed layer
47	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
48	USE prtctl ! Print control
49	USE in_out_manager ! I/O manager
50	USE iom ! I/O manager library
51	USE lib_mpp ! MPP library
52	USE wrk_nemo ! work arrays
53	USE timing ! Timing
54	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
55	#if defined key_agrif
56	USE agrif_opa_interp
57	USE agrif_opa_update
58	#endif
59	!CE
60	#if defined key_lim3
61	USE ice, ONLY: htm_i, ht_i
62	#endif
63	#if defined key_lim2 \|\| defined key_cice
64	USE sbc_ice, ONLY: ht_i
65	#endif
66	!CE
67
68	IMPLICIT NONE
69	PRIVATE
70
71	PUBLIC zdf_tke ! routine called in step module
72	PUBLIC zdf_tke_init ! routine called in opa module
73	PUBLIC tke_rst ! routine called in step module
74
75	LOGICAL , PUBLIC, PARAMETER :: lk_zdftke = .TRUE. !: TKE vertical mixing flag
76
77	! !! Namelist namzdf_tke
78	LOGICAL :: ln_mxl0 ! mixing length scale surface value as function of wind stress or not
79	!CE
80	INTEGER :: nn_mxl0 ! type of scaling under sea-ice (=0/1/2/3)
81	REAL(wp) :: rn_hice ! ice thickness value when scaling under sea-ice
82	!CE
83	INTEGER :: nn_mxl ! type of mixing length (=0/1/2/3)
84	REAL(wp) :: rn_mxl0 ! surface min value of mixing length (kappaz_o=0.40.1 m) [m]
85	INTEGER :: nn_pdl ! Prandtl number or not (ratio avt/avm) (=0/1)
86	REAL(wp) :: rn_ediff ! coefficient for avt: avt=rn_ediffmxlsqrt(e)
87	REAL(wp) :: rn_ediss ! coefficient of the Kolmogoroff dissipation
88	REAL(wp) :: rn_ebb ! coefficient of the surface input of tke
89	REAL(wp) :: rn_emin ! minimum value of tke [m2/s2]
90	REAL(wp) :: rn_emin0 ! surface minimum value of tke [m2/s2]
91	REAL(wp) :: rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
92	INTEGER :: nn_etau ! type of depth penetration of surface tke (=0/1/2/3)
93	INTEGER :: nn_htau ! type of tke profile of penetration (=0/1)
94	REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean
95	LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not
96	REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells
97
98	REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
99	REAL(wp) :: rmxl_min ! minimum mixing length value (deduced from rn_ediff and rn_emin values) [m]
100	REAL(wp) :: rhftau_add = 1.e-3_wp ! add offset applied to HF part of taum (nn_etau=3)
101	REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3)
102
103	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau)
104	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation
105	#if defined key_c1d
106	! !! 1D cfg only ('key_c1d')
107	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_dis, e_mix !: dissipation and mixing turbulent lengh scales
108	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers
109	#endif
110
111	!! * Substitutions
112	# include "domzgr_substitute.h90"
113	# include "vectopt_loop_substitute.h90"
114	!!----------------------------------------------------------------------
115	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
116	!! $Id$
117	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
118	!!----------------------------------------------------------------------
119	CONTAINS
120
121	INTEGER FUNCTION zdf_tke_alloc()
122	!!----------------------------------------------------------------------
123	!! * FUNCTION zdf_tke_alloc *
124	!!----------------------------------------------------------------------
125	ALLOCATE( &
126	#if defined key_c1d
127	& e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) , &
128	& e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) , &
129	#endif
130	& htau (jpi,jpj) , dissl(jpi,jpj,jpk) , STAT= zdf_tke_alloc )
131	!
132	IF( lk_mpp ) CALL mpp_sum ( zdf_tke_alloc )
133	IF( zdf_tke_alloc /= 0 ) CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
134	!
135	END FUNCTION zdf_tke_alloc
136
137
138	SUBROUTINE zdf_tke( kt )
139	!!----------------------------------------------------------------------
140	!! * ROUTINE zdf_tke *
141	!!
142	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
143	!! coefficients using a turbulent closure scheme (TKE).
144	!!
145	!! ** Method : The time evolution of the turbulent kinetic energy (tke)
146	!! is computed from a prognostic equation :
147	!! d(en)/dt = avm (d(u)/dz)**2 ! shear production
148	!! + d( avm d(en)/dz )/dz ! diffusion of tke
149	!! + avt N^2 ! stratif. destruc.
150	!! - rn_ediss / emxl en**(2/3) ! Kolmogoroff dissipation
151	!! with the boundary conditions:
152	!! surface: en = max( rn_emin0, rn_ebb * taum )
153	!! bottom : en = rn_emin
154	!! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
155	!!
156	!! The now Turbulent kinetic energy is computed using the following
157	!! time stepping: implicit for vertical diffusion term, linearized semi
158	!! implicit for kolmogoroff dissipation term, and explicit forward for
159	!! both buoyancy and shear production terms. Therefore a tridiagonal
160	!! linear system is solved. Note that buoyancy and shear terms are
161	!! discretized in a energy conserving form (Bruchard 2002).
162	!!
163	!! The dissipative and mixing length scale are computed from en and
164	!! the stratification (see tke_avn)
165	!!
166	!! The now vertical eddy vicosity and diffusivity coefficients are
167	!! given by:
168	!! avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
169	!! avt = max( avmb, pdl * avm )
170	!! eav = max( avmb, avm )
171	!! where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
172	!! given by an empirical funtion of the localRichardson number if nn_pdl=1
173	!!
174	!! ** Action : compute en (now turbulent kinetic energy)
175	!! update avt, avmu, avmv (before vertical eddy coef.)
176	!!
177	!! References : Gaspar et al., JGR, 1990,
178	!! Blanke and Delecluse, JPO, 1991
179	!! Mellor and Blumberg, JPO 2004
180	!! Axell, JGR, 2002
181	!! Bruchard OM 2002
182	!!----------------------------------------------------------------------
183	INTEGER, INTENT(in) :: kt ! ocean time step
184	!!----------------------------------------------------------------------
185	!
186	IF( kt /= nit000 ) THEN ! restore before value to compute tke
187	avt (:,:,:) = avt_k (:,:,:)
188	avm (:,:,:) = avm_k (:,:,:)
189	avmu(:,:,:) = avmu_k(:,:,:)
190	avmv(:,:,:) = avmv_k(:,:,:)
191	ENDIF
192	!
193	CALL tke_tke ! now tke (en)
194	!
195	CALL tke_avn ! now avt, avm, avmu, avmv
196	!
197	avt_k (:,:,:) = avt (:,:,:)
198	avm_k (:,:,:) = avm (:,:,:)
199	avmu_k(:,:,:) = avmu(:,:,:)
200	avmv_k(:,:,:) = avmv(:,:,:)
201	!
202	#if defined key_agrif
203	! Update child grid f => parent grid
204	IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only
205	#endif
206	!
207	END SUBROUTINE zdf_tke
208
209
210	SUBROUTINE tke_tke
211	!!----------------------------------------------------------------------
212	!! * ROUTINE tke_tke *
213	!!
214	!! ** Purpose : Compute the now Turbulente Kinetic Energy (TKE)
215	!!
216	!! ** Method : - TKE surface boundary condition
217	!! - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
218	!! - source term due to shear (saved in avmu, avmv arrays)
219	!! - Now TKE : resolution of the TKE equation by inverting
220	!! a tridiagonal linear system by a "methode de chasse"
221	!! - increase TKE due to surface and internal wave breaking
222	!!
223	!! ** Action : - en : now turbulent kinetic energy)
224	!! - avmu, avmv : production of TKE by shear at u and v-points
225	!! (= Kz dz[Ub] * dz[Un] )
226	!! ---------------------------------------------------------------------
227	INTEGER :: ji, jj, jk ! dummy loop arguments
228	!!bfr INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! temporary scalar
229	!!bfr INTEGER :: ikbt, ikbumm1, ikbvmm1 ! temporary scalar
230	REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3
231	REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient
232	REAL(wp) :: zbbrau, zesh2 ! temporary scalars
233	REAL(wp) :: zfact1, zfact2, zfact3 ! - -
234	REAL(wp) :: ztx2 , zty2 , zcof ! - -
235	REAL(wp) :: ztau , zdif ! - -
236	REAL(wp) :: zus , zwlc , zind ! - -
237	REAL(wp) :: zzd_up, zzd_lw ! - -
238	!!bfr REAL(wp) :: zebot ! - -
239	INTEGER , POINTER, DIMENSION(:,: ) :: imlc
240	REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc
241	REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw
242	!!--------------------------------------------------------------------
243	!
244	IF( nn_timing == 1 ) CALL timing_start('tke_tke')
245	!
246	CALL wrk_alloc( jpi,jpj, imlc ) ! integer
247	CALL wrk_alloc( jpi,jpj, zhlc )
248	CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )
249	!
250	zbbrau = rn_ebb / rau0 ! Local constant initialisation
251	zfact1 = -.5_wp * rdt
252	zfact2 = 1.5_wp * rdt * rn_ediss
253	zfact3 = 0.5_wp * rn_ediss
254	!
255	!
256	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
257	! ! Surface boundary condition on tke
258	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
259	IF ( ln_isfcav ) THEN
260	DO jj = 2, jpjm1 ! en(mikt(ji,jj)) = rn_emin
261	DO ji = fs_2, fs_jpim1 ! vector opt.
262	en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1)
263	END DO
264	END DO
265	END IF
266	DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0)
267	DO ji = fs_2, fs_jpim1 ! vector opt.
268	en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
269	END DO
270	END DO
271
272	!!bfr - start commented area
273	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
274	! ! Bottom boundary condition on tke
275	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
276	!
277	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
278	! Tests to date have found the bottom boundary condition on tke to have very little effect.
279	! The condition is coded here for completion but commented out until there is proof that the
280	! computational cost is justified
281	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
282	! en(bot) = (rn_ebb0/rau0)0.5sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
283	!CDIR NOVERRCHK
284	!! DO jj = 2, jpjm1
285	!CDIR NOVERRCHK
286	!! DO ji = fs_2, fs_jpim1 ! vector opt.
287	!! ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
288	!! bfrua(ji ,jj) * ub(ji ,jj,mbku(ji ,jj) )
289	!! zty2 = bfrva(ji,jj ) * vb(ji,jj ,mbkv(ji,jj )) + &
290	!! bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) )
291	!! zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.
292	!! en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)
293	!! END DO
294	!! END DO
295	!!bfr - end commented area
296	!
297	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
298	IF( ln_lc ) THEN ! Langmuir circulation source term added to tke (Axell JGR 2002)
299	! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
300	!
301	! !* total energy produce by LC : cumulative sum over jk
302	zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * fsdepw(:,:,1) * fse3w(:,:,1)
303	DO jk = 2, jpk
304	zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
305	END DO
306	! !* finite Langmuir Circulation depth
307	zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
308	imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land)
309	DO jk = jpkm1, 2, -1
310	DO jj = 1, jpj ! Last w-level at which zpelc>=0.5usus
311	DO ji = 1, jpi ! with us=0.016*wind(starting from jpk-1)
312	zus = zcof * taum(ji,jj)
313	IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk
314	END DO
315	END DO
316	END DO
317	! ! finite LC depth
318	DO jj = 1, jpj
319	DO ji = 1, jpi
320	zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
321	END DO
322	END DO
323	zcof = 0.016 / SQRT( zrhoa * zcdrag )
324	!CDIR NOVERRCHK
325	DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en
326	!CDIR NOVERRCHK
327	DO jj = 2, jpjm1
328	!CDIR NOVERRCHK
329	DO ji = fs_2, fs_jpim1 ! vector opt.
330	zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift
331	! ! vertical velocity due to LC
332	zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
333	zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
334	! ! TKE Langmuir circulation source term
335	en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.fr_i(ji,jj) ) ( zwlc * zwlc * zwlc ) / &
336	& zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
337	END DO
338	END DO
339	END DO
340	!
341	ENDIF
342	!
343	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
344	! ! Now Turbulent kinetic energy (output in en)
345	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
346	! ! Resolution of a tridiagonal linear system by a "methode de chasse"
347	! ! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ).
348	! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
349	!
350	DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form)
351	DO jj = 1, jpj ! here avmu, avmv used as workspace
352	DO ji = 1, jpi
353	avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) &
354	& * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) &
355	& / ( fse3uw_n(ji,jj,jk) &
356	& * fse3uw_b(ji,jj,jk) )
357	avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) &
358	& * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) &
359	& / ( fse3vw_n(ji,jj,jk) &
360	& * fse3vw_b(ji,jj,jk) )
361	END DO
362	END DO
363	END DO
364	!
365	DO jk = 2, jpkm1 !* Matrix and right hand side in en
366	DO jj = 2, jpjm1
367	DO ji = fs_2, fs_jpim1 ! vector opt.
368	zcof = zfact1 * tmask(ji,jj,jk)
369	# if defined key_zdftmx_new
370	! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability)
371	zzd_up = zcof * ( MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp ) ) & ! upper diagonal
372	& / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk ) )
373	zzd_lw = zcof * ( MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp ) ) & ! lower diagonal
374	& / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )
375	# else
376	zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal
377	& / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk ) )
378	zzd_lw = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & ! lower diagonal
379	& / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )
380	# endif
381	! ! shear prod. at w-point weightened by mask
382	zesh2 = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &
383	& + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
384	!
385	zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw)
386	zd_lw(ji,jj,jk) = zzd_lw
387	zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
388	!
389	! ! right hand side in en
390	en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) &
391	& + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) &
392	& * wmask(ji,jj,jk)
393	END DO
394	END DO
395	END DO
396	! !* Matrix inversion from level 2 (tke prescribed at level 1)
397	DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
398	DO jj = 2, jpjm1
399	DO ji = fs_2, fs_jpim1 ! vector opt.
400	zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
401	END DO
402	END DO
403	END DO
404	!
405	! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
406	DO jj = 2, jpjm1
407	DO ji = fs_2, fs_jpim1 ! vector opt.
408	zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke
409	END DO
410	END DO
411	DO jk = 3, jpkm1
412	DO jj = 2, jpjm1
413	DO ji = fs_2, fs_jpim1 ! vector opt.
414	zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
415	END DO
416	END DO
417	END DO
418	!
419	! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
420	DO jj = 2, jpjm1
421	DO ji = fs_2, fs_jpim1 ! vector opt.
422	en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
423	END DO
424	END DO
425	DO jk = jpk-2, 2, -1
426	DO jj = 2, jpjm1
427	DO ji = fs_2, fs_jpim1 ! vector opt.
428	en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
429	END DO
430	END DO
431	END DO
432	DO jk = 2, jpkm1 ! set the minimum value of tke
433	DO jj = 2, jpjm1
434	DO ji = fs_2, fs_jpim1 ! vector opt.
435	en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
436	END DO
437	END DO
438	END DO
439
440	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
441	! ! TKE due to surface and internal wave breaking
442	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
443	IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction)
444	DO jk = 2, jpkm1
445	DO jj = 2, jpjm1
446	DO ji = fs_2, fs_jpim1 ! vector opt.
447	en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) &
448	& * MAX(0.,1._wp - 2.fr_i(ji,jj) ) wmask(ji,jj,jk) * tmask(ji,jj,1)
449	END DO
450	END DO
451	END DO
452	ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction)
453	DO jj = 2, jpjm1
454	DO ji = fs_2, fs_jpim1 ! vector opt.
455	jk = nmln(ji,jj)
456	en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) &
457	& * MAX(0.,1._wp - 2.fr_i(ji,jj) ) wmask(ji,jj,jk) * tmask(ji,jj,1)
458	END DO
459	END DO
460	ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability)
461	!CDIR NOVERRCHK
462	DO jk = 2, jpkm1
463	!CDIR NOVERRCHK
464	DO jj = 2, jpjm1
465	!CDIR NOVERRCHK
466	DO ji = fs_2, fs_jpim1 ! vector opt.
467	ztx2 = utau(ji-1,jj ) + utau(ji,jj)
468	zty2 = vtau(ji ,jj-1) + vtau(ji,jj)
469	ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress
470	zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean
471	zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications...
472	en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) &
473	& * MAX(0.,1._wp - 2.fr_i(ji,jj) ) wmask(ji,jj,jk) * tmask(ji,jj,1)
474	END DO
475	END DO
476	END DO
477	ENDIF
478	CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged)
479	!
480	CALL wrk_dealloc( jpi,jpj, imlc ) ! integer
481	CALL wrk_dealloc( jpi,jpj, zhlc )
482	CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )
483	!
484	IF( nn_timing == 1 ) CALL timing_stop('tke_tke')
485	!
486	END SUBROUTINE tke_tke
487
488
489	SUBROUTINE tke_avn
490	!!----------------------------------------------------------------------
491	!! * ROUTINE tke_avn *
492	!!
493	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
494	!!
495	!! ** Method : At this stage, en, the now TKE, is known (computed in
496	!! the tke_tke routine). First, the now mixing lenth is
497	!! computed from en and the strafification (N^2), then the mixings
498	!! coefficients are computed.
499	!! - Mixing length : a first evaluation of the mixing lengh
500	!! scales is:
501	!! mxl = sqrt(2*en) / N
502	!! where N is the brunt-vaisala frequency, with a minimum value set
503	!! to rmxl_min (rn_mxl0) in the interior (surface) ocean.
504	!! The mixing and dissipative length scale are bound as follow :
505	!! nn_mxl=0 : mxl bounded by the distance to surface and bottom.
506	!! zmxld = zmxlm = mxl
507	!! nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
508	!! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
509	!! less than 1 (\|d/dz(mxl)\|<1) and zmxld = zmxlm = mxl
510	!! nn_mxl=3 : mxl is bounded from the surface to the bottom usings
511	!! \|d/dz(xml)\|<1 to obtain lup, and from the bottom to
512	!! the surface to obtain ldown. the resulting length
513	!! scales are:
514	!! zmxld = sqrt( lup * ldown )
515	!! zmxlm = min ( lup , ldown )
516	!! - Vertical eddy viscosity and diffusivity:
517	!! avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
518	!! avt = max( avmb, pdlr * avm )
519	!! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
520	!!
521	!! ** Action : - avt : now vertical eddy diffusivity (w-point)
522	!! - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
523	!!----------------------------------------------------------------------
524	INTEGER :: ji, jj, jk ! dummy loop indices
525	REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars
526	REAL(wp) :: zdku, zpdlr, zri, zsqen ! - -
527	REAL(wp) :: zdkv, zemxl, zemlm, zemlp ! - -
528	REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
529	!!--------------------------------------------------------------------
530	!
531	IF( nn_timing == 1 ) CALL timing_start('tke_avn')
532
533	CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )
534
535	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
536	! ! Mixing length
537	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
538	!
539	! !* Buoyancy length scale: l=sqrt(2e/n*2)
540	!
541	! initialisation of interior minimum value (avoid a 2d loop with mikt)
542	zmxlm(:,:,:) = rmxl_min
543	zmxld(:,:,:) = rmxl_min
544	!
545	IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn2.e5taum/(rau0*g)
546	!
547	zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
548	!
549	SELECT CASE( nn_mxl0 ) ! Type of scaling under sea-ice
550	!
551	CASE( 0 ) ! No scaling under sea-ice
552	zmxlm(:,:,1) = zraug * taum(:,:) * tmask(:,:,1)
553	!
554	CASE( 1 ) ! scaling with constant sea-ice thickness
555	zmxlm(:,:,1) = rn_hice
556	!
557	CASE( 2 ) ! scaling with mean sea-ice thickness
558	#if defined key_lim3
559	zmxlm(:,:,1) = ( 1. - fr_i(:,:) ) * zraug * taum(:,:) * tmask(:,:,1) + fr_i(:,:) * htm_i(:,:) * 2.
560	#elif ( defined key_lim2 \|\| defined key_cice )
561	zmxlm(:,:,1) = ( 1. - fr_i(:,:) ) * zraug * taum(:,:) * tmask(:,:,1) + fr_i(:,:) * MAXVAL( ht_i, 3 )
562	#else
563	zmxlm(:,:,1) = zraug * taum(:,:) * tmask(:,:,1)
564	#endif
565	!
566	CASE( 3 ) ! scaling with max sea-ice thickness
567	#if defined key_lim3 \|\| defined key_lim2 \|\| defined key_cice
568	zmxlm(:,:,1) = ( 1. - fr_i(:,:) ) * zraug * taum(:,:) * tmask(:,:,1) + fr_i(:,:) * MAXVAL( ht_i, 3 )
569	#else
570	zmxlm(:,:,1) = zraug * taum(:,:) * tmask(:,:,1)
571	#endif
572	!
573	END SELECT
574	!
575	DO jj = 2, jpjm1
576	DO ji = fs_2, fs_jpim1
577	zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
578	END DO
579	END DO
580	!
581	ELSE
582	zmxlm(:,:,1) = rn_mxl0
583	ENDIF
584	!
585	IF( iom_use('mixlength') ) THEN
586	CALL lbc_lnk( zmxlm(:,:,1) , 'T', 1. )
587	CALL iom_put( "mixlength" , zmxlm(:,:,1) )
588	ENDIF
589	!
590	!CDIR NOVERRCHK
591	DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2)
592	!CDIR NOVERRCHK
593	DO jj = 2, jpjm1
594	!CDIR NOVERRCHK
595	DO ji = fs_2, fs_jpim1 ! vector opt.
596	zrn2 = MAX( rn2(ji,jj,jk), rsmall )
597	zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) )
598	END DO
599	END DO
600	END DO
601	!
602	! !* Physical limits for the mixing length
603	!
604	zmxld(:,:,1 ) = zmxlm(:,:,1) ! surface set to the minimum value
605	zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value
606	!
607	SELECT CASE ( nn_mxl )
608	!
609	! where wmask = 0 set zmxlm == fse3w
610	CASE ( 0 ) ! bounded by the distance to surface and bottom
611	DO jk = 2, jpkm1
612	DO jj = 2, jpjm1
613	DO ji = fs_2, fs_jpim1 ! vector opt.
614	zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), &
615	& fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) )
616	! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
617	zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
618	zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
619	END DO
620	END DO
621	END DO
622	!
623	CASE ( 1 ) ! bounded by the vertical scale factor
624	DO jk = 2, jpkm1
625	DO jj = 2, jpjm1
626	DO ji = fs_2, fs_jpim1 ! vector opt.
627	zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
628	zmxlm(ji,jj,jk) = zemxl
629	zmxld(ji,jj,jk) = zemxl
630	END DO
631	END DO
632	END DO
633	!
634	CASE ( 2 ) ! \|dk[xml]\| bounded by e3t :
635	DO jk = 2, jpkm1 ! from the surface to the bottom :
636	DO jj = 2, jpjm1
637	DO ji = fs_2, fs_jpim1 ! vector opt.
638	zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
639	END DO
640	END DO
641	END DO
642	DO jk = jpkm1, 2, -1 ! from the bottom to the surface :
643	DO jj = 2, jpjm1
644	DO ji = fs_2, fs_jpim1 ! vector opt.
645	zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
646	zmxlm(ji,jj,jk) = zemxl
647	zmxld(ji,jj,jk) = zemxl
648	END DO
649	END DO
650	END DO
651	!
652	CASE ( 3 ) ! lup and ldown, \|dk[xml]\| bounded by e3t :
653	DO jk = 2, jpkm1 ! from the surface to the bottom : lup
654	DO jj = 2, jpjm1
655	DO ji = fs_2, fs_jpim1 ! vector opt.
656	zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
657	END DO
658	END DO
659	END DO
660	DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown
661	DO jj = 2, jpjm1
662	DO ji = fs_2, fs_jpim1 ! vector opt.
663	zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
664	END DO
665	END DO
666	END DO
667	!CDIR NOVERRCHK
668	DO jk = 2, jpkm1
669	!CDIR NOVERRCHK
670	DO jj = 2, jpjm1
671	!CDIR NOVERRCHK
672	DO ji = fs_2, fs_jpim1 ! vector opt.
673	zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) )
674	zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
675	zmxlm(ji,jj,jk) = zemlm
676	zmxld(ji,jj,jk) = zemlp
677	END DO
678	END DO
679	END DO
680	!
681	END SELECT
682	!
683	# if defined key_c1d
684	e_dis(:,:,:) = zmxld(:,:,:) ! c1d configuration : save mixing and dissipation turbulent length scales
685	e_mix(:,:,:) = zmxlm(:,:,:)
686	# endif
687
688	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
689	! ! Vertical eddy viscosity and diffusivity (avmu, avmv, avt)
690	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
691	!CDIR NOVERRCHK
692	DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points
693	!CDIR NOVERRCHK
694	DO jj = 2, jpjm1
695	!CDIR NOVERRCHK
696	DO ji = fs_2, fs_jpim1 ! vector opt.
697	zsqen = SQRT( en(ji,jj,jk) )
698	zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen
699	avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk)
700	avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
701	dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
702	END DO
703	END DO
704	END DO
705	CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged)
706	!
707	DO jk = 2, jpkm1 !* vertical eddy viscosity at wu- and wv-points
708	DO jj = 2, jpjm1
709	DO ji = fs_2, fs_jpim1 ! vector opt.
710	avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk)
711	avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk)
712	END DO
713	END DO
714	END DO
715	CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions
716	!
717	IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt
718	DO jk = 2, jpkm1
719	DO jj = 2, jpjm1
720	DO ji = fs_2, fs_jpim1 ! vector opt.
721	zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
722	! ! shear
723	zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) ) &
724	& + avmu(ji ,jj,jk) * ( un(ji ,jj,jk-1) - un(ji ,jj,jk) ) * ( ub(ji ,jj,jk-1) - ub(ji ,jj,jk) )
725	zdkv = avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) * ( vb(ji,jj-1,jk-1) - vb(ji,jj-1,jk) ) &
726	& + avmv(ji,jj ,jk) * ( vn(ji,jj ,jk-1) - vn(ji,jj ,jk) ) * ( vb(ji,jj ,jk-1) - vb(ji,jj ,jk) )
727	! ! local Richardson number
728	zri = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear )
729	zpdlr = MAX( 0.1_wp, 0.2 / MAX( 0.2 , zri ) )
730	!!gm and even better with the use of the "true" ri_crit=0.22222... (this change the results!)
731	!!gm zpdlr = MAX( 0.1_wp, ri_crit / MAX( ri_crit , zri ) )
732	avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
733	# if defined key_c1d
734	e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number
735	e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk) ! c1d config. : save Ri
736	# endif
737	END DO
738	END DO
739	END DO
740	ENDIF
741	CALL lbc_lnk( avt, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged)
742
743	IF(ln_ctl) THEN
744	CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
745	CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke - u: ', mask1=umask, &
746	& tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
747	ENDIF
748	!
749	CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )
750	!
751	IF( nn_timing == 1 ) CALL timing_stop('tke_avn')
752	!
753	END SUBROUTINE tke_avn
754
755
756	SUBROUTINE zdf_tke_init
757	!!----------------------------------------------------------------------
758	!! * ROUTINE zdf_tke_init *
759	!!
760	!! ** Purpose : Initialization of the vertical eddy diffivity and
761	!! viscosity when using a tke turbulent closure scheme
762	!!
763	!! ** Method : Read the namzdf_tke namelist and check the parameters
764	!! called at the first timestep (nit000)
765	!!
766	!! ** input : Namlist namzdf_tke
767	!!
768	!! ** Action : Increase by 1 the nstop flag is setting problem encounter
769	!!----------------------------------------------------------------------
770	INTEGER :: ji, jj, jk ! dummy loop indices
771	INTEGER :: ios, ierr
772	!!
773	NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , &
774	& rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , &
775	!CE
776	& nn_mxl0 , rn_hice, &
777	!CE
778	& rn_mxl0 , nn_pdl , ln_lc , rn_lc , &
779	& nn_etau , nn_htau , rn_efr
780	!!----------------------------------------------------------------------
781	!
782	REWIND( numnam_ref ) ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy
783	READ ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
784	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist', lwp )
785
786	REWIND( numnam_cfg ) ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
787	READ ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
788	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist', lwp )
789	IF(lwm) WRITE ( numond, namzdf_tke )
790	!
791	ri_cri = 2._wp / ( 2._wp + rn_ediss / rn_ediff ) ! resulting critical Richardson number
792	# if defined key_zdftmx_new
793	! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used
794	rn_emin = 1.e-10_wp
795	rmxl_min = 1.e-03_wp
796	IF(lwp) THEN ! Control print
797	WRITE(numout,*)
798	WRITE(numout,*) 'zdf_tke_init : New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
799	WRITE(numout,*) '~~~~~~~~~~~~'
800	ENDIF
801	# else
802	rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) ) ! resulting minimum length to recover molecular viscosity
803	# endif
804	!
805	IF(lwp) THEN !* Control print
806	WRITE(numout,*)
807	WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
808	WRITE(numout,*) '~~~~~~~~~~~~'
809	WRITE(numout,*) ' Namelist namzdf_tke : set tke mixing parameters'
810	WRITE(numout,*) ' coef. to compute avt rn_ediff = ', rn_ediff
811	WRITE(numout,*) ' Kolmogoroff dissipation coef. rn_ediss = ', rn_ediss
812	WRITE(numout,*) ' tke surface input coef. rn_ebb = ', rn_ebb
813	WRITE(numout,*) ' minimum value of tke rn_emin = ', rn_emin
814	WRITE(numout,*) ' surface minimum value of tke rn_emin0 = ', rn_emin0
815	WRITE(numout,*) ' background shear (>0) rn_bshear = ', rn_bshear
816	WRITE(numout,*) ' mixing length type nn_mxl = ', nn_mxl
817	WRITE(numout,*) ' prandl number flag nn_pdl = ', nn_pdl
818	WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0
819	!CE
820	IF( ln_mxl0 ) THEN
821	WRITE(numout,*) ' type of scaling under sea-ice nn_mxl0 = ', nn_mxl0
822	IF( nn_mxl0 == 1 ) &
823	WRITE(numout,*) ' ice thickness value when scaling under sea-ice rn_hice = ', rn_hice
824	ENDIF
825	!CE
826	WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0
827	WRITE(numout,*) ' flag to take into acc. Langmuir circ. ln_lc = ', ln_lc
828	WRITE(numout,*) ' coef to compute verticla velocity of LC rn_lc = ', rn_lc
829	WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau
830	WRITE(numout,*) ' flag for computation of exp. tke profile nn_htau = ', nn_htau
831	WRITE(numout,*) ' fraction of en which pene. the thermocline rn_efr = ', rn_efr
832	WRITE(numout,*)
833	WRITE(numout,*) ' critical Richardson nb with your parameters ri_cri = ', ri_cri
834	ENDIF
835	!
836	! ! allocate tke arrays
837	IF( zdf_tke_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
838	!
839	! !* Check of some namelist values
840	IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' )
841	IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' )
842	IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
843	IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
844
845	IF( ln_mxl0 ) THEN
846	IF(lwp) WRITE(numout,*) ' use a surface mixing length = F(stress) : set rn_mxl0 = rmxl_min'
847	rn_mxl0 = rmxl_min
848	ENDIF
849
850	IF( nn_etau == 2 ) THEN
851	ierr = zdf_mxl_alloc()
852	nmln(:,:) = nlb10 ! Initialization of nmln
853	ENDIF
854
855	! !* depth of penetration of surface tke
856	IF( nn_etau /= 0 ) THEN
857	SELECT CASE( nn_htau ) ! Choice of the depth of penetration
858	CASE( 0 ) ! constant depth penetration (here 10 meters)
859	htau(:,:) = 10._wp
860	CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees
861	htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) )
862	END SELECT
863	ENDIF
864	! !* set vertical eddy coef. to the background value
865	DO jk = 1, jpk
866	avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
867	avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
868	avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
869	avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
870	END DO
871	dissl(:,:,:) = 1.e-12_wp
872	!
873	CALL tke_rst( nit000, 'READ' ) !* read or initialize all required files
874	!
875	END SUBROUTINE zdf_tke_init
876
877
878	SUBROUTINE tke_rst( kt, cdrw )
879	!!---------------------------------------------------------------------
880	!! * ROUTINE tke_rst *
881	!!
882	!! ** Purpose : Read or write TKE file (en) in restart file
883	!!
884	!! ** Method : use of IOM library
885	!! if the restart does not contain TKE, en is either
886	!! set to rn_emin or recomputed
887	!!----------------------------------------------------------------------
888	INTEGER , INTENT(in) :: kt ! ocean time-step
889	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
890	!
891	INTEGER :: jit, jk ! dummy loop indices
892	INTEGER :: id1, id2, id3, id4, id5, id6 ! local integers
893	!!----------------------------------------------------------------------
894	!
895	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
896	! ! ---------------
897	IF( ln_rstart ) THEN !* Read the restart file
898	id1 = iom_varid( numror, 'en' , ldstop = .FALSE. )
899	id2 = iom_varid( numror, 'avt' , ldstop = .FALSE. )
900	id3 = iom_varid( numror, 'avm' , ldstop = .FALSE. )
901	id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
902	id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
903	id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
904	!
905	IF( id1 > 0 ) THEN ! 'en' exists
906	CALL iom_get( numror, jpdom_autoglo, 'en', en )
907	IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN ! all required arrays exist
908	CALL iom_get( numror, jpdom_autoglo, 'avt' , avt )
909	CALL iom_get( numror, jpdom_autoglo, 'avm' , avm )
910	CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu )
911	CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv )
912	CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
913	ELSE ! one at least array is missing
914	CALL tke_avn ! compute avt, avm, avmu, avmv and dissl (approximation)
915	ENDIF
916	ELSE ! No TKE array found: initialisation
917	IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
918	en (:,:,:) = rn_emin * tmask(:,:,:)
919	CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation)
920	!
921	avt_k (:,:,:) = avt (:,:,:)
922	avm_k (:,:,:) = avm (:,:,:)
923	avmu_k(:,:,:) = avmu(:,:,:)
924	avmv_k(:,:,:) = avmv(:,:,:)
925	!
926	DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO
927	ENDIF
928	ELSE !* Start from rest
929	en(:,:,:) = rn_emin * tmask(:,:,:)
930	DO jk = 1, jpk ! set the Kz to the background value
931	avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
932	avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
933	avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
934	avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
935	END DO
936	ENDIF
937	!
938	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
939	! ! -------------------
940	IF(lwp) WRITE(numout,*) '---- tke-rst ----'
941	CALL iom_rstput( kt, nitrst, numrow, 'en' , en )
942	CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k )
943	CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k )
944	CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
945	CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
946	CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
947	!
948	ENDIF
949	!
950	END SUBROUTINE tke_rst
951
952	#else
953	!!----------------------------------------------------------------------
954	!! Dummy module : NO TKE scheme
955	!!----------------------------------------------------------------------
956	LOGICAL, PUBLIC, PARAMETER :: lk_zdftke = .FALSE. !: TKE flag
957	CONTAINS
958	SUBROUTINE zdf_tke_init ! Dummy routine
959	END SUBROUTINE zdf_tke_init
960	SUBROUTINE zdf_tke( kt ) ! Dummy routine
961	WRITE(,) 'zdf_tke: You should not have seen this print! error?', kt
962	END SUBROUTINE zdf_tke
963	SUBROUTINE tke_rst( kt, cdrw )
964	CHARACTER(len=*) :: cdrw
965	WRITE(,) 'tke_rst: You should not have seen this print! error?', kt, cdwr
966	END SUBROUTINE tke_rst
967	#endif
968
969	!!======================================================================
970	END MODULE zdftke

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

Download in other formats: