Context Navigation

zdftke.F90 @ 8762

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

Download in other formats: