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

Last change on this file since 3231 was 3231, checked in by smasson, 12 years ago
dev_NEMO_MERGE_2011: supress TARGET attribute for tsa and use work arrays
Property svn:keywords set to `Id`
File size: 44.2 KB

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

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

Download in other formats: