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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

zdftke.F90 in branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

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

Last change on this file since 3229 was 3229, checked in by charris, 12 years ago
Added timing calls to most significant routines in LDF, SBC and ZDF.
Property svn:keywords set to `Id`
File size: 44.5 KB

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

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

Download in other formats: