Context Navigation

zdftke.F90 @ 2690

Last change on this file since 2690 was 2690, checked in by gm, 13 years ago
dynamic mem: #785 ; homogeneization of the coding style associated with dyn allocation
Property svn:keywords set to `Id`
File size: 44.1 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	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau)
86	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation
87	#if defined key_c1d
88	! !! 1D cfg only ('key_c1d')
89	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_dis, e_mix !: dissipation and mixing turbulent lengh scales
90	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers
91	#endif
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, ONLY: zmpdl => ua , zmxlm => va , zmxld => ta ! (ua,va,ta) used as workspace
474	!!
475	INTEGER :: ji, jj, jk ! dummy loop indices
476	REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars
477	REAL(wp) :: zdku, zpdlr, zri, zsqen ! - -
478	REAL(wp) :: zdkv, zemxl, zemlm, zemlp ! - -
479	!!--------------------------------------------------------------------
480
481	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
482	! ! Mixing length
483	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
484	!
485	! !* Buoyancy length scale: l=sqrt(2e/n*2)
486	!
487	IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn2.e5taum/(rau0*g)
488	zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
489	zmxlm(:,:,1) = MAX( rn_mxl0, zraug * taum(:,:) )
490	ELSE ! surface set to the minimum value
491	zmxlm(:,:,1) = rn_mxl0
492	ENDIF
493	zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value
494	!
495	!CDIR NOVERRCHK
496	DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2)
497	!CDIR NOVERRCHK
498	DO jj = 2, jpjm1
499	!CDIR NOVERRCHK
500	DO ji = fs_2, fs_jpim1 ! vector opt.
501	zrn2 = MAX( rn2(ji,jj,jk), rsmall )
502	zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) )
503	END DO
504	END DO
505	END DO
506	!
507	! !* Physical limits for the mixing length
508	!
509	zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the zmxlm value
510	zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value
511	!
512	SELECT CASE ( nn_mxl )
513	!
514	CASE ( 0 ) ! bounded by the distance to surface and bottom
515	DO jk = 2, jpkm1
516	DO jj = 2, jpjm1
517	DO ji = fs_2, fs_jpim1 ! vector opt.
518	zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk), &
519	& fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) )
520	zmxlm(ji,jj,jk) = zemxl
521	zmxld(ji,jj,jk) = zemxl
522	END DO
523	END DO
524	END DO
525	!
526	CASE ( 1 ) ! bounded by the vertical scale factor
527	DO jk = 2, jpkm1
528	DO jj = 2, jpjm1
529	DO ji = fs_2, fs_jpim1 ! vector opt.
530	zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
531	zmxlm(ji,jj,jk) = zemxl
532	zmxld(ji,jj,jk) = zemxl
533	END DO
534	END DO
535	END DO
536	!
537	CASE ( 2 ) ! \|dk[xml]\| bounded by e3t :
538	DO jk = 2, jpkm1 ! from the surface to the bottom :
539	DO jj = 2, jpjm1
540	DO ji = fs_2, fs_jpim1 ! vector opt.
541	zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
542	END DO
543	END DO
544	END DO
545	DO jk = jpkm1, 2, -1 ! from the bottom to the surface :
546	DO jj = 2, jpjm1
547	DO ji = fs_2, fs_jpim1 ! vector opt.
548	zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
549	zmxlm(ji,jj,jk) = zemxl
550	zmxld(ji,jj,jk) = zemxl
551	END DO
552	END DO
553	END DO
554	!
555	CASE ( 3 ) ! lup and ldown, \|dk[xml]\| bounded by e3t :
556	DO jk = 2, jpkm1 ! from the surface to the bottom : lup
557	DO jj = 2, jpjm1
558	DO ji = fs_2, fs_jpim1 ! vector opt.
559	zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
560	END DO
561	END DO
562	END DO
563	DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown
564	DO jj = 2, jpjm1
565	DO ji = fs_2, fs_jpim1 ! vector opt.
566	zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
567	END DO
568	END DO
569	END DO
570	!CDIR NOVERRCHK
571	DO jk = 2, jpkm1
572	!CDIR NOVERRCHK
573	DO jj = 2, jpjm1
574	!CDIR NOVERRCHK
575	DO ji = fs_2, fs_jpim1 ! vector opt.
576	zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) )
577	zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
578	zmxlm(ji,jj,jk) = zemlm
579	zmxld(ji,jj,jk) = zemlp
580	END DO
581	END DO
582	END DO
583	!
584	END SELECT
585	!
586	# if defined key_c1d
587	e_dis(:,:,:) = zmxld(:,:,:) ! c1d configuration : save mixing and dissipation turbulent length scales
588	e_mix(:,:,:) = zmxlm(:,:,:)
589	# endif
590
591	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
592	! ! Vertical eddy viscosity and diffusivity (avmu, avmv, avt)
593	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
594	!CDIR NOVERRCHK
595	DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points
596	!CDIR NOVERRCHK
597	DO jj = 2, jpjm1
598	!CDIR NOVERRCHK
599	DO ji = fs_2, fs_jpim1 ! vector opt.
600	zsqen = SQRT( en(ji,jj,jk) )
601	zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen
602	avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)
603	avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
604	dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
605	END DO
606	END DO
607	END DO
608	CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged)
609	!
610	DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points
611	DO jj = 2, jpjm1
612	DO ji = fs_2, fs_jpim1 ! vector opt.
613	avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk)
614	avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk)
615	END DO
616	END DO
617	END DO
618	CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions
619	!
620	IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt
621	DO jk = 2, jpkm1
622	DO jj = 2, jpjm1
623	DO ji = fs_2, fs_jpim1 ! vector opt.
624	zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
625	! ! shear
626	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) ) &
627	& + avmu(ji ,jj,jk) * ( un(ji ,jj,jk-1) - un(ji ,jj,jk) ) * ( ub(ji ,jj,jk-1) - ub(ji ,jj,jk) )
628	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) ) &
629	& + avmv(ji,jj ,jk) * ( vn(ji,jj ,jk-1) - vn(ji,jj ,jk) ) * ( vb(ji,jj ,jk-1) - vb(ji,jj ,jk) )
630	! ! local Richardson number
631	zri = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear )
632	zpdlr = MAX( 0.1_wp, 0.2 / MAX( 0.2 , zri ) )
633	!!gm and even better with the use of the "true" ri_crit=0.22222... (this change the results!)
634	!!gm zpdlr = MAX( 0.1_wp, ri_crit / MAX( ri_crit , zri ) )
635	avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
636	# if defined key_c1d
637	e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number
638	e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri
639	# endif
640	END DO
641	END DO
642	END DO
643	ENDIF
644	CALL lbc_lnk( avt, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged)
645
646	IF(ln_ctl) THEN
647	CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
648	CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke - u: ', mask1=umask, &
649	& tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
650	ENDIF
651	!
652	END SUBROUTINE tke_avn
653
654
655	SUBROUTINE zdf_tke_init
656	!!----------------------------------------------------------------------
657	!! * ROUTINE zdf_tke_init *
658	!!
659	!! ** Purpose : Initialization of the vertical eddy diffivity and
660	!! viscosity when using a tke turbulent closure scheme
661	!!
662	!! ** Method : Read the namzdf_tke namelist and check the parameters
663	!! called at the first timestep (nit000)
664	!!
665	!! ** input : Namlist namzdf_tke
666	!!
667	!! ** Action : Increase by 1 the nstop flag is setting problem encounter
668	!!----------------------------------------------------------------------
669	INTEGER :: ji, jj, jk ! dummy loop indices
670	!!
671	NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , &
672	& rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , &
673	& rn_mxl0 , nn_pdl , ln_lc , rn_lc , &
674	& nn_etau , nn_htau , rn_efr
675	!!----------------------------------------------------------------------
676	!
677	REWIND ( numnam ) !* Read Namelist namzdf_tke : Turbulente Kinetic Energy
678	READ ( numnam, namzdf_tke )
679	!
680	ri_cri = 2._wp / ( 2._wp + rn_ediss / rn_ediff ) ! resulting critical Richardson number
681	rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) ) ! resulting minimum length to recover molecular viscosity
682	!
683	IF(lwp) THEN !* Control print
684	WRITE(numout,*)
685	WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
686	WRITE(numout,*) '~~~~~~~~~~~~'
687	WRITE(numout,*) ' Namelist namzdf_tke : set tke mixing parameters'
688	WRITE(numout,*) ' coef. to compute avt rn_ediff = ', rn_ediff
689	WRITE(numout,*) ' Kolmogoroff dissipation coef. rn_ediss = ', rn_ediss
690	WRITE(numout,*) ' tke surface input coef. rn_ebb = ', rn_ebb
691	WRITE(numout,*) ' minimum value of tke rn_emin = ', rn_emin
692	WRITE(numout,*) ' surface minimum value of tke rn_emin0 = ', rn_emin0
693	WRITE(numout,*) ' background shear (>0) rn_bshear = ', rn_bshear
694	WRITE(numout,*) ' mixing length type nn_mxl = ', nn_mxl
695	WRITE(numout,*) ' prandl number flag nn_pdl = ', nn_pdl
696	WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0
697	WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0
698	WRITE(numout,*) ' flag to take into acc. Langmuir circ. ln_lc = ', ln_lc
699	WRITE(numout,*) ' coef to compute verticla velocity of LC rn_lc = ', rn_lc
700	WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau
701	WRITE(numout,*) ' flag for computation of exp. tke profile nn_htau = ', nn_htau
702	WRITE(numout,*) ' fraction of en which pene. the thermocline rn_efr = ', rn_efr
703	WRITE(numout,*)
704	WRITE(numout,*) ' critical Richardson nb with your parameters ri_cri = ', ri_cri
705	ENDIF
706	!
707	! ! allocate tke arrays
708	IF( zdf_tke_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
709	!
710	! !* Check of some namelist values
711	IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' )
712	IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' )
713	IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
714	#if ! key_coupled
715	IF( nn_etau == 3 ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
716	#endif
717
718	IF( ln_mxl0 ) THEN
719	IF(lwp) WRITE(numout,*) ' use a surface mixing length = F(stress) : set rn_mxl0 = rmxl_min'
720	rn_mxl0 = rmxl_min
721	ENDIF
722
723	IF( nn_etau == 2 ) CALL zdf_mxl( nit000 ) ! Initialization of nmln
724
725	! !* depth of penetration of surface tke
726	IF( nn_etau /= 0 ) THEN
727	SELECT CASE( nn_htau ) ! Choice of the depth of penetration
728	CASE( 0 ) ! constant depth penetration (here 10 meters)
729	htau(:,:) = 10._wp
730	CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees
731	htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) )
732	END SELECT
733	ENDIF
734	! !* set vertical eddy coef. to the background value
735	DO jk = 1, jpk
736	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
737	avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
738	avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
739	avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
740	END DO
741	dissl(:,:,:) = 1.e-12_wp
742	!
743	CALL tke_rst( nit000, 'READ' ) !* read or initialize all required files
744	!
745	END SUBROUTINE zdf_tke_init
746
747
748	SUBROUTINE tke_rst( kt, cdrw )
749	!!---------------------------------------------------------------------
750	!! * ROUTINE tke_rst *
751	!!
752	!! ** Purpose : Read or write TKE file (en) in restart file
753	!!
754	!! ** Method : use of IOM library
755	!! if the restart does not contain TKE, en is either
756	!! set to rn_emin or recomputed
757	!!----------------------------------------------------------------------
758	INTEGER , INTENT(in) :: kt ! ocean time-step
759	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
760	!
761	INTEGER :: jit, jk ! dummy loop indices
762	INTEGER :: id1, id2, id3, id4, id5, id6 ! local integers
763	!!----------------------------------------------------------------------
764	!
765	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
766	! ! ---------------
767	IF( ln_rstart ) THEN !* Read the restart file
768	id1 = iom_varid( numror, 'en' , ldstop = .FALSE. )
769	id2 = iom_varid( numror, 'avt' , ldstop = .FALSE. )
770	id3 = iom_varid( numror, 'avm' , ldstop = .FALSE. )
771	id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
772	id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
773	id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
774	!
775	IF( id1 > 0 ) THEN ! 'en' exists
776	CALL iom_get( numror, jpdom_autoglo, 'en', en )
777	IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN ! all required arrays exist
778	CALL iom_get( numror, jpdom_autoglo, 'avt' , avt )
779	CALL iom_get( numror, jpdom_autoglo, 'avm' , avm )
780	CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu )
781	CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv )
782	CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
783	ELSE ! one at least array is missing
784	CALL tke_avn ! compute avt, avm, avmu, avmv and dissl (approximation)
785	ENDIF
786	ELSE ! No TKE array found: initialisation
787	IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
788	en (:,:,:) = rn_emin * tmask(:,:,:)
789	CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation)
790	DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO
791	ENDIF
792	ELSE !* Start from rest
793	en(:,:,:) = rn_emin * tmask(:,:,:)
794	DO jk = 1, jpk ! set the Kz to the background value
795	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
796	avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
797	avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
798	avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
799	END DO
800	ENDIF
801	!
802	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
803	! ! -------------------
804	IF(lwp) WRITE(numout,*) '---- tke-rst ----'
805	CALL iom_rstput( kt, nitrst, numrow, 'en' , en )
806	CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt )
807	CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm )
808	CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu )
809	CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv )
810	CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
811	!
812	ENDIF
813	!
814	END SUBROUTINE tke_rst
815
816	#else
817	!!----------------------------------------------------------------------
818	!! Dummy module : NO TKE scheme
819	!!----------------------------------------------------------------------
820	LOGICAL, PUBLIC, PARAMETER :: lk_zdftke = .FALSE. !: TKE flag
821	CONTAINS
822	SUBROUTINE zdf_tke_init ! Dummy routine
823	END SUBROUTINE zdf_tke_init
824	SUBROUTINE zdf_tke( kt ) ! Dummy routine
825	WRITE(,) 'zdf_tke: You should not have seen this print! error?', kt
826	END SUBROUTINE zdf_tke
827	SUBROUTINE tke_rst( kt, cdrw )
828	CHARACTER(len=*) :: cdrw
829	WRITE(,) 'tke_rst: You should not have seen this print! error?', kt, cdwr
830	END SUBROUTINE tke_rst
831	#endif
832
833	!!======================================================================
834	END MODULE zdftke

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

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

Download in other formats: