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 branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 3186

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

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

Download in other formats: