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/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 2636

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

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

Download in other formats: