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

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