Context Navigation

zdftke.F90 @ 15540

Last change on this file since 15540 was 15540, checked in by sparonuz, 3 years ago
Mixed precision version, tested up to 30 years on ORCA2.
Property svn:keywords set to `Id`
File size: 50.1 KB

Line
1	MODULE zdftke
2	!!======================================================================
3	!! * MODULE zdftke *
4	!! Ocean physics: vertical mixing coefficient computed from the tke
5	!! turbulent closure parameterization
6	!!=====================================================================
7	!! History : OPA ! 1991-03 (b. blanke) Original code
8	!! 7.0 ! 1991-11 (G. Madec) bug fix
9	!! 7.1 ! 1992-10 (G. Madec) new mixing length and eav
10	!! 7.2 ! 1993-03 (M. Guyon) symetrical conditions
11	!! 7.3 ! 1994-08 (G. Madec, M. Imbard) nn_pdl flag
12	!! 7.5 ! 1996-01 (G. Madec) s-coordinates
13	!! 8.0 ! 1997-07 (G. Madec) lbc
14	!! 8.1 ! 1999-01 (E. Stretta) new option for the mixing length
15	!! NEMO 1.0 ! 2002-06 (G. Madec) add tke_init routine
16	!! - ! 2004-10 (C. Ethe ) 1D configuration
17	!! 2.0 ! 2006-07 (S. Masson) distributed restart using iom
18	!! 3.0 ! 2008-05 (C. Ethe, G.Madec) : update TKE physics:
19	!! ! - tke penetration (wind steering)
20	!! ! - suface condition for tke & mixing length
21	!! ! - Langmuir cells
22	!! - ! 2008-05 (J.-M. Molines, G. Madec) 2D form of avtb
23	!! - ! 2008-06 (G. Madec) style + DOCTOR name for namelist parameters
24	!! - ! 2008-12 (G. Reffray) stable discretization of the production term
25	!! 3.2 ! 2009-06 (G. Madec, S. Masson) TKE restart compatible with key_cpl
26	!! ! + cleaning of the parameters + bugs correction
27	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
28	!! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability
29	!! 4.0 ! 2017-04 (G. Madec) remove CPP ddm key & avm at t-point only
30	!! - ! 2017-05 (G. Madec) add top/bottom friction as boundary condition
31	!! 4.2 ! 2020-12 (G. Madec, E. Clementi) add wave coupling
32	! ! following Couvelard et al., 2019
33	!!----------------------------------------------------------------------
34
35	!!----------------------------------------------------------------------
36	!! zdf_tke : update momentum and tracer Kz from a tke scheme
37	!! tke_tke : tke time stepping: update tke at now time step (en)
38	!! tke_avn : compute mixing length scale and deduce avm and avt
39	!! zdf_tke_init : initialization, namelist read, and parameters control
40	!! tke_rst : read/write tke restart in ocean restart file
41	!!----------------------------------------------------------------------
42	USE oce ! ocean: dynamics and active tracers variables
43	USE phycst ! physical constants
44	USE dom_oce ! domain: ocean
45	USE domvvl ! domain: variable volume layer
46	USE sbc_oce ! surface boundary condition: ocean
47	USE zdfdrg ! vertical physics: top/bottom drag coef.
48	USE zdfmxl ! vertical physics: mixed layer
49	#if defined key_si3
50	USE ice, ONLY: hm_i, h_i
51	#endif
52	#if defined key_cice
53	USE sbc_ice, ONLY: h_i
54	#endif
55	!
56	USE in_out_manager ! I/O manager
57	USE iom ! I/O manager library
58	USE lib_mpp ! MPP library
59	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
60	USE prtctl ! Print control
61	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
62	USE sbcwave ! Surface boundary waves
63
64	IMPLICIT NONE
65	PRIVATE
66
67	PUBLIC zdf_tke ! routine called in step module
68	PUBLIC zdf_tke_init ! routine called in opa module
69	PUBLIC tke_rst ! routine called in step module
70
71	! !! Namelist namzdf_tke
72	LOGICAL :: ln_mxl0 ! mixing length scale surface value as function of wind stress or not
73	LOGICAL :: ln_mxhsw ! mixing length scale surface value as a fonction of wave height
74	INTEGER :: nn_mxlice ! type of scaling under sea-ice (=0/1/2/3)
75	REAL(wp) :: rn_mxlice ! ice thickness value when scaling under sea-ice
76	INTEGER :: nn_mxl ! type of mixing length (=0/1/2/3)
77	REAL(wp) :: rn_mxl0 ! surface min value of mixing length (kappaz_o=0.40.1 m) [m]
78	INTEGER :: nn_pdl ! Prandtl number or not (ratio avt/avm) (=0/1)
79	REAL(wp) :: rn_ediff ! coefficient for avt: avt=rn_ediffmxlsqrt(e)
80	REAL(wp) :: rn_ediss ! coefficient of the Kolmogoroff dissipation
81	REAL(wp) :: rn_ebb ! coefficient of the surface input of tke
82	REAL(wp) :: rn_emin ! minimum value of tke [m2/s2]
83	REAL(wp) :: rn_emin0 ! surface minimum value of tke [m2/s2]
84	REAL(wp) :: rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
85	INTEGER :: nn_etau ! type of depth penetration of surface tke (=0/1/2/3)
86	INTEGER :: nn_htau ! type of tke profile of penetration (=0/1)
87	INTEGER :: nn_bc_surf! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling
88	INTEGER :: nn_bc_bot ! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling
89	REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean
90	LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not
91	REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells
92	INTEGER :: nn_eice ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3)
93
94	REAL(dp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
95	REAL(dp) :: rmxl_min ! minimum mixing length value (deduced from rn_ediff and rn_emin values) [m]
96	REAL(dp) :: rhftau_add = 1.e-3_wp ! add offset applied to HF part of taum (nn_etau=3)
97	REAL(dp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3)
98
99	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau)
100	REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation
101	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation
102
103	!! * Substitutions
104	# include "do_loop_substitute.h90"
105	# include "domzgr_substitute.h90"
106	!!----------------------------------------------------------------------
107	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
108	!! $Id$
109	!! Software governed by the CeCILL license (see ./LICENSE)
110	!!----------------------------------------------------------------------
111	CONTAINS
112
113	INTEGER FUNCTION zdf_tke_alloc()
114	!!----------------------------------------------------------------------
115	!! * FUNCTION zdf_tke_alloc *
116	!!----------------------------------------------------------------------
117	ALLOCATE( htau(jpi,jpj) , dissl(jpi,jpj,jpk) , apdlr(jpi,jpj,jpk) , STAT= zdf_tke_alloc )
118	!
119	CALL mpp_sum ( 'zdftke', zdf_tke_alloc )
120	IF( zdf_tke_alloc /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tke_alloc: failed to allocate arrays' )
121	!
122	END FUNCTION zdf_tke_alloc
123
124
125	SUBROUTINE zdf_tke( kt, Kbb, Kmm, p_sh2, p_avm, p_avt )
126	!!----------------------------------------------------------------------
127	!! * ROUTINE zdf_tke *
128	!!
129	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
130	!! coefficients using a turbulent closure scheme (TKE).
131	!!
132	!! ** Method : The time evolution of the turbulent kinetic energy (tke)
133	!! is computed from a prognostic equation :
134	!! d(en)/dt = avm (d(u)/dz)**2 ! shear production
135	!! + d( avm d(en)/dz )/dz ! diffusion of tke
136	!! + avt N^2 ! stratif. destruc.
137	!! - rn_ediss / emxl en**(2/3) ! Kolmogoroff dissipation
138	!! with the boundary conditions:
139	!! surface: en = max( rn_emin0, rn_ebb * taum )
140	!! bottom : en = rn_emin
141	!! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
142	!!
143	!! The now Turbulent kinetic energy is computed using the following
144	!! time stepping: implicit for vertical diffusion term, linearized semi
145	!! implicit for kolmogoroff dissipation term, and explicit forward for
146	!! both buoyancy and shear production terms. Therefore a tridiagonal
147	!! linear system is solved. Note that buoyancy and shear terms are
148	!! discretized in a energy conserving form (Bruchard 2002).
149	!!
150	!! The dissipative and mixing length scale are computed from en and
151	!! the stratification (see tke_avn)
152	!!
153	!! The now vertical eddy vicosity and diffusivity coefficients are
154	!! given by:
155	!! avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
156	!! avt = max( avmb, pdl * avm )
157	!! eav = max( avmb, avm )
158	!! where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
159	!! given by an empirical funtion of the localRichardson number if nn_pdl=1
160	!!
161	!! ** Action : compute en (now turbulent kinetic energy)
162	!! update avt, avm (before vertical eddy coef.)
163	!!
164	!! References : Gaspar et al., JGR, 1990,
165	!! Blanke and Delecluse, JPO, 1991
166	!! Mellor and Blumberg, JPO 2004
167	!! Axell, JGR, 2002
168	!! Bruchard OM 2002
169	!!----------------------------------------------------------------------
170	INTEGER , INTENT(in ) :: kt ! ocean time step
171	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
172	REAL(dp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: p_sh2 ! shear production term
173	REAL(dp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)
174	!!----------------------------------------------------------------------
175	!
176	CALL tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt ) ! now tke (en)
177	!
178	CALL tke_avn( Kbb, Kmm, p_avm, p_avt ) ! now avt, avm, dissl
179	!
180	END SUBROUTINE zdf_tke
181
182
183	SUBROUTINE tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt )
184	!!----------------------------------------------------------------------
185	!! * ROUTINE tke_tke *
186	!!
187	!! ** Purpose : Compute the now Turbulente Kinetic Energy (TKE)
188	!!
189	!! ** Method : - TKE surface boundary condition
190	!! - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
191	!! - source term due to shear (= Kz dz[Ub] * dz[Un] )
192	!! - Now TKE : resolution of the TKE equation by inverting
193	!! a tridiagonal linear system by a "methode de chasse"
194	!! - increase TKE due to surface and internal wave breaking
195	!! NB: when sea-ice is present, both LC parameterization
196	!! and TKE penetration are turned off when the ice fraction
197	!! is smaller than 0.25
198	!!
199	!! ** Action : - en : now turbulent kinetic energy)
200	!! ---------------------------------------------------------------------
201	USE zdf_oce , ONLY : en ! ocean vertical physics
202	!!
203	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
204	REAL(dp), DIMENSION(A2D(nn_hls),jpk) , INTENT(in ) :: p_sh2 ! shear production term
205	REAL(dp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points)
206	!
207	INTEGER :: ji, jj, jk ! dummy loop arguments
208	REAL(wp) :: zetop, zebot, zmsku, zmskv ! local scalars
209	REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3
210	REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient
211	REAL(wp) :: zbbrau, zbbirau, zri ! local scalars
212	REAL(wp) :: zfact1, zfact2, zfact3 ! - -
213	REAL(wp) :: ztx2 , zty2 , zcof ! - -
214	REAL(wp) :: ztau , zdif ! - -
215	REAL(wp) :: zus, zwlc ! - -
216	REAL(dp) :: zind ! - -
217	REAL(wp) :: zzd_up, zzd_lw ! - -
218	REAL(wp) :: ztaui, ztauj, z1_norm
219	INTEGER , DIMENSION(A2D(nn_hls)) :: imlc
220	REAL(wp), DIMENSION(A2D(nn_hls)) :: zice_fra, zhlc, zus3, zWlc2
221	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpelc, zdiag, zd_up, zd_lw
222	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: ztmp ! for diags
223	!!--------------------------------------------------------------------
224	!
225	zbbrau = rn_ebb / rho0 ! Local constant initialisation
226	zbbirau = 3.75_wp / rho0
227	zfact1 = -.5_wp * rn_Dt
228	zfact2 = 1.5_wp * rn_Dt * rn_ediss
229	zfact3 = 0.5_wp * rn_ediss
230	!
231	zpelc(:,:,:) = 0._wp ! need to be initialised in case ln_lc is not used
232	!
233	! ice fraction considered for attenuation of langmuir & wave breaking
234	SELECT CASE ( nn_eice )
235	CASE( 0 ) ; zice_fra(:,:) = 0._wp
236	CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i(A2D(nn_hls)) * 10._wp )
237	CASE( 2 ) ; zice_fra(:,:) = fr_i(A2D(nn_hls))
238	CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp )
239	END SELECT
240	!
241	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
242	! ! Surface/top/bottom boundary condition on tke
243	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
244	!
245	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
246	en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) )
247	zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)
248	zd_lw(ji,jj,1) = 1._wp
249	zd_up(ji,jj,1) = 0._wp
250	END_2D
251	!
252	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
253	! ! Bottom boundary condition on tke
254	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
255	!
256	! en(bot) = (ebb0/rho0)0.5sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
257	! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75)
258	! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2
259	!
260	IF( .NOT.ln_drg_OFF ) THEN !== friction used as top/bottom boundary condition on TKE
261	!
262	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! bottom friction
263	zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
264	zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
265	! ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
266	zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT( ( zmsku( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )*2 &
267	& + ( zmskv( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )*2 )
268	en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj)
269	END_2D
270	IF( ln_isfcav ) THEN
271	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! top friction
272	zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
273	zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
274	! ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
275	zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )*2 &
276	& + ( zmskv( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )*2 )
277	! (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) = 1 where ice shelves are present
278	en(ji,jj,mikt(ji,jj)) = en(ji,jj,1) * tmask(ji,jj,1) &
279	& + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj)
280	END_2D
281	ENDIF
282	!
283	ENDIF
284	!
285	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
286	IF( ln_lc ) THEN ! Langmuir circulation source term added to tke (Axell JGR 2002)
287	! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
288	!
289	! !* Langmuir velocity scale
290	!
291	IF ( cpl_sdrftx ) THEN ! Surface Stokes Drift available
292	! ! Craik-Leibovich velocity scale Wlc = ( u* u_s )^1/2 with u* = (taum/rho0)^1/2
293	! ! associated kinetic energy : 1/2 (Wlc)^2 = u* u_s
294	! ! more precisely, it is the dot product that must be used :
295	! ! 1/2 (W_lc)^2 = MAX( u* u_s + v* v_s , 0 ) only the positive part
296	!!gm ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s
297	!!gm ! so we will overestimate the LC velocity.... !!gm I will do the work if !LC have an effect !
298	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
299	!!XC zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)2 +vt0sd(ji,jj)2 ) )
300	zWlc2(ji,jj) = 0.5_wp * ( ut0sd(ji,jj)2 +vt0sd(ji,jj)2 )
301	END_2D
302	!
303	! Projection of Stokes drift in the wind stress direction
304	!
305	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
306	ztaui = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) )
307	ztauj = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) )
308	z1_norm = 1._wp / MAX( SQRT(ztauiztaui+ztaujztauj), 1.e-12 ) * tmask(ji,jj,1)
309	zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)ztaui + vt0sd(ji,jj)ztauj, 0._wp ) )**2
310	END_2D
311	ELSE ! Surface Stokes drift deduced from surface stress
312	! ! Wlc = u_s with u_s = 0.016*U_10m, the surface stokes drift (Axell 2002, Eq.44)
313	! ! using \|tau\| = rho_air Cd \|U_10m\|^2 , it comes:
314	! ! Wlc = 0.016 * [\|tau\|/(rho_air Cdrag) ]^1/2 and thus:
315	! ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 \|tau\| /( rho_air Cdrag )
316	zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) ! to convert stress in 10m wind using a constant drag
317	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
318	zWlc2(ji,jj) = zcof * taum(ji,jj)
319	END_2D
320	!
321	ENDIF
322	!
323	! !* Depth of the LC circulation (Axell 2002, Eq.47)
324	! !- LHS of Eq.47
325	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
326	zpelc(ji,jj,1) = MAX( rn2b(ji,jj,1), 0._wp ) * gdepw(ji,jj,1,Kmm) * e3w(ji,jj,1,Kmm)
327	END_2D
328	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpk )
329	zpelc(ji,jj,jk) = zpelc(ji,jj,jk-1) + &
330	& MAX( rn2b(ji,jj,jk), 0._wp ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
331	END_3D
332	!
333	! !- compare LHS to RHS of Eq.47
334	imlc(:,:) = mbkt(A2D(nn_hls)) + 1 ! Initialization to the number of w ocean point (=2 over land)
335	DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )
336	IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) ) imlc(ji,jj) = jk
337	END_3D
338	! ! finite LC depth
339	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
340	zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm)
341	END_2D
342	!
343	zcof = 0.016 / SQRT( zrhoa * zcdrag )
344	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
345	zus = SQRT( 2. * zWlc2(ji,jj) ) ! Stokes drift
346	zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
347	END_2D
348	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* TKE Langmuir circulation source term added to en
349	IF ( zus3(ji,jj) /= 0._wp ) THEN
350	IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN
351	! ! vertical velocity due to LC
352	zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )
353	! ! TKE Langmuir circulation source term
354	en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
355	ENDIF
356	ENDIF
357	END_3D
358	!
359	ENDIF
360	!
361	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
362	! ! Now Turbulent kinetic energy (output in en)
363	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
364	! ! Resolution of a tridiagonal linear system by a "methode de chasse"
365	! ! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ).
366	! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
367	!
368	IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri )
369	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
370	! ! local Richardson number
371	IF (rn2b(ji,jj,jk) <= 0.0_wp) then
372	zri = 0.0_wp
373	ELSE
374	zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
375	ENDIF
376	! ! inverse of Prandtl number
377	apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) )
378	END_3D
379	ENDIF
380	!
381	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* Matrix and right hand side in en
382	zcof = zfact1 * tmask(ji,jj,jk)
383	! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
384	! ! eddy coefficient (ensure numerical stability)
385	zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal
386	& / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk ,Kmm) )
387	zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal
388	& / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk ,Kmm) )
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) * wmask(ji,jj,jk)
393	!
394	! ! right hand side in en
395	en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * ( p_sh2(ji,jj,jk) & ! shear
396	& - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification
397	& + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation
398	& ) * wmask(ji,jj,jk)
399	END_3D
400	!
401	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
402	! ! Surface boundary condition on tke if
403	! ! coupling with waves
404	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
405	!
406	IF ( cpl_phioc .and. ln_phioc ) THEN
407	SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves
408
409	CASE ( 0 ) ! Dirichlet BC
410	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! en(1) = rn_ebb taum / rho0 (min value rn_emin0)
411	IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp
412	en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )*(2./3.) ) tmask(ji,jj,1)
413	zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) ! choose to keep coherence with former estimation of
414	END_2D
415
416	CASE ( 1 ) ! Neumann BC
417	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
418	IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp
419	en(ji,jj,2) = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm)
420	en(ji,jj,1) = en(ji,jj,2) + (2 * e3t(ji,jj,1,Kmm) * phioc(ji,jj)/rho0) / ( p_avm(ji,jj,1) + p_avm(ji,jj,2) )
421	zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2)
422	zdiag(ji,jj,1) = 1._wp
423	zd_lw(ji,jj,2) = 0._wp
424	END_2D
425
426	END SELECT
427
428	ENDIF
429	!
430	! !* Matrix inversion from level 2 (tke prescribed at level 1)
431	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
432	zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
433	END_3D
434	!XC : commented to allow for neumann boundary condition
435	! DO_2D( 0, 0, 0, 0 )
436	! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke
437	! END_2D
438	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
439	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)
440	END_3D
441	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
442	en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
443	END_2D
444	DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 )
445	en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
446	END_3D
447	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! set the minimum value of tke
448	en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
449	END_3D
450	!
451	! Kolmogorov energy of dissipation (W/kg)
452	! ediss = Cesqrt(en)/Len
453	! dissl = sqrt(en)/L
454	IF( iom_use('ediss_k') ) THEN
455	ALLOCATE( ztmp(A2D(nn_hls),jpk) )
456	ztmp(:,:,:) = zfact3 * dissl * en * wmask
457	CALL lbc_lnk( 'zdftke', ztmp, 'W', 1._wp )
458	CALL iom_put( 'ediss_k', ztmp )
459	DEALLOCATE( ztmp )
460	ENDIF
461	!
462	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
463	! ! TKE due to surface and internal wave breaking
464	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
465	!!gm BUG : in the exp remove the depth of ssh !!!
466	!!gm i.e. use gde3w in argument (gdepw(:,:,:,Kmm))
467	!
468	! penetration is partly switched off below sea-ice if nn_eice/=0
469	!
470	IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction)
471	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
472	en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) &
473	& * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
474	END_3D
475	ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction)
476	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
477	jk = nmln(ji,jj)
478	en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) &
479	& * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
480	END_2D
481	ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability)
482	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
483	ztx2 = utau(ji-1,jj ) + utau(ji,jj)
484	zty2 = vtau(ji ,jj-1) + vtau(ji,jj)
485	ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress
486	zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean
487	zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications...
488	en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) &
489	& * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
490	END_3D
491	ENDIF
492	!
493	END SUBROUTINE tke_tke
494
495
496	SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt )
497	!!----------------------------------------------------------------------
498	!! * ROUTINE tke_avn *
499	!!
500	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
501	!!
502	!! ** Method : At this stage, en, the now TKE, is known (computed in
503	!! the tke_tke routine). First, the now mixing lenth is
504	!! computed from en and the strafification (N^2), then the mixings
505	!! coefficients are computed.
506	!! - Mixing length : a first evaluation of the mixing lengh
507	!! scales is:
508	!! mxl = sqrt(2*en) / N
509	!! where N is the brunt-vaisala frequency, with a minimum value set
510	!! to rmxl_min (rn_mxl0) in the interior (surface) ocean.
511	!! The mixing and dissipative length scale are bound as follow :
512	!! nn_mxl=0 : mxl bounded by the distance to surface and bottom.
513	!! zmxld = zmxlm = mxl
514	!! nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
515	!! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
516	!! less than 1 (\|d/dz(mxl)\|<1) and zmxld = zmxlm = mxl
517	!! nn_mxl=3 : mxl is bounded from the surface to the bottom usings
518	!! \|d/dz(xml)\|<1 to obtain lup, and from the bottom to
519	!! the surface to obtain ldown. the resulting length
520	!! scales are:
521	!! zmxld = sqrt( lup * ldown )
522	!! zmxlm = min ( lup , ldown )
523	!! - Vertical eddy viscosity and diffusivity:
524	!! avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
525	!! avt = max( avmb, pdlr * avm )
526	!! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
527	!!
528	!! ** Action : - avt, avm : now vertical eddy diffusivity and viscosity (w-point)
529	!!----------------------------------------------------------------------
530	USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d ! ocean vertical physics
531	!!
532	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
533	REAL(dp), DIMENSION(:,:,:), INTENT( out) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points)
534	!
535	INTEGER :: ji, jj, jk ! dummy loop indices
536	REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars
537	REAL(wp) :: zsqen ! - -
538	REAL(dp) :: zdku, zdkv ! - -
539	REAL(wp) :: zemxl, zemlm, zemlp ! - -
540	REAL(dp) :: zmaxice ! - -
541	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zmxlm, zmxld ! 3D workspace
542	!!--------------------------------------------------------------------
543	!
544	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
545	! ! Mixing length
546	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
547	!
548	! !* Buoyancy length scale: l=sqrt(2e/n*2)
549	!
550	! initialisation of interior minimum value (avoid a 2d loop with mikt)
551	zmxlm(:,:,:) = rmxl_min
552	zmxld(:,:,:) = rmxl_min
553	!
554	IF(ln_sdw .AND. ln_mxhsw) THEN
555	zmxlm(:,:,1)= vkarmn * MAX ( 1.6 * hsw(:,:) , 0.02 ) ! surface mixing length = F(wave height)
556	! from terray et al 1999 and mellor and blumberg 2004 it should be 0.85 and not 1.6
557	zcoef = vkarmn * ( (rn_ediffrn_ediss)*0.25 ) / rn_ediff
558	zmxlm(:,:,1)= zcoef * MAX ( 1.6 * hsw(:,:) , 0.02 ) ! surface mixing length = F(wave height)
559	ELSE
560	!
561	IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn2.e5taum/(rho0*g)
562	!
563	zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
564	#if ! defined key_si3 && ! defined key_cice
565	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! No sea-ice
566	zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
567	END_2D
568	#else
569	SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice
570	!
571	CASE( 0 ) ! No scaling under sea-ice
572	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
573	zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
574	END_2D
575	!
576	CASE( 1 ) ! scaling with constant sea-ice thickness
577	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
578	zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
579	& fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1)
580	END_2D
581	!
582	CASE( 2 ) ! scaling with mean sea-ice thickness
583	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
584	#if defined key_si3
585	zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
586	& fr_i(ji,jj) * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1)
587	#elif defined key_cice
588	zmaxice = MAXVAL( h_i(ji,jj,:) )
589	zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
590	& fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1)
591	#endif
592	END_2D
593	!
594	CASE( 3 ) ! scaling with max sea-ice thickness
595	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
596	zmaxice = MAXVAL( h_i(ji,jj,:) )
597	zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
598	& fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1)
599	END_2D
600	!
601	END SELECT
602	#endif
603	!
604	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
605	zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
606	END_2D
607	!
608	ELSE
609	zmxlm(:,:,1) = rn_mxl0
610	ENDIF
611	ENDIF
612	!
613	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
614	zrn2 = MAX( rn2(ji,jj,jk), rsmall )
615	zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) )
616	END_3D
617	!
618	! !* Physical limits for the mixing length
619	!
620	zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the minimum value
621	zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value
622	!
623	SELECT CASE ( nn_mxl )
624	!
625	!!gm Not sure of that coding for ISF....
626	! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm)
627	CASE ( 0 ) ! bounded by the distance to surface and bottom
628	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
629	zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk), &
630	& gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) )
631	! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
632	zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) &
633	& + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
634	zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) &
635	& + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
636	END_3D
637	!
638	CASE ( 1 ) ! bounded by the vertical scale factor
639	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
640	zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) )
641	zmxlm(ji,jj,jk) = zemxl
642	zmxld(ji,jj,jk) = zemxl
643	END_3D
644	!
645	CASE ( 2 ) ! \|dk[xml]\| bounded by e3t :
646	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! from the surface to the bottom :
647	zmxlm(ji,jj,jk) = &
648	& MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
649	END_3D
650	DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! from the bottom to the surface :
651	zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
652	zmxlm(ji,jj,jk) = zemxl
653	zmxld(ji,jj,jk) = zemxl
654	END_3D
655	!
656	CASE ( 3 ) ! lup and ldown, \|dk[xml]\| bounded by e3t :
657	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! from the surface to the bottom : lup
658	zmxld(ji,jj,jk) = &
659	& MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
660	END_3D
661	DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! from the bottom to the surface : ldown
662	zmxlm(ji,jj,jk) = &
663	& MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
664	END_3D
665	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
666	zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) )
667	zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
668	zmxlm(ji,jj,jk) = zemlm
669	zmxld(ji,jj,jk) = zemlp
670	END_3D
671	!
672	END SELECT
673	!
674	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
675	! ! Vertical eddy viscosity and diffusivity (avm and avt)
676	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
677	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* vertical eddy viscosity & diffivity at w-points
678	zsqen = SQRT( en(ji,jj,jk) )
679	zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen
680	p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk)
681	p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
682	dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
683	END_3D
684	!
685	!
686	IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt
687	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
688	p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
689	END_3D
690	ENDIF
691	!
692	IF(sn_cfctl%l_prtctl) THEN
693	CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk)
694	CALL prt_ctl( tab3d_1=p_avm, clinfo1=' tke - m: ', kdim=jpk )
695	ENDIF
696	!
697	END SUBROUTINE tke_avn
698
699
700	SUBROUTINE zdf_tke_init( Kmm )
701	!!----------------------------------------------------------------------
702	!! * ROUTINE zdf_tke_init *
703	!!
704	!! ** Purpose : Initialization of the vertical eddy diffivity and
705	!! viscosity when using a tke turbulent closure scheme
706	!!
707	!! ** Method : Read the namzdf_tke namelist and check the parameters
708	!! called at the first timestep (nit000)
709	!!
710	!! ** input : Namlist namzdf_tke
711	!!
712	!! ** Action : Increase by 1 the nstop flag is setting problem encounter
713	!!----------------------------------------------------------------------
714	USE zdf_oce , ONLY : ln_zdfiwm ! Internal Wave Mixing flag
715	!!
716	INTEGER, INTENT(in) :: Kmm ! time level index
717	INTEGER :: ji, jj, jk ! dummy loop indices
718	INTEGER :: ios
719	!!
720	NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , &
721	& rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , &
722	& rn_mxl0 , nn_mxlice, rn_mxlice, &
723	& nn_pdl , ln_lc , rn_lc , &
724	& nn_etau , nn_htau , rn_efr , nn_eice , &
725	& nn_bc_surf, nn_bc_bot, ln_mxhsw
726	!!----------------------------------------------------------------------
727	!
728	READ ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
729	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist' )
730
731	READ ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
732	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist' )
733	IF(lwm) WRITE ( numond, namzdf_tke )
734	!
735	ri_cri = 2._wp / ( 2._wp + rn_ediss / rn_ediff ) ! resulting critical Richardson number
736	!
737	IF(lwp) THEN !* Control print
738	WRITE(numout,*)
739	WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
740	WRITE(numout,*) '~~~~~~~~~~~~'
741	WRITE(numout,*) ' Namelist namzdf_tke : set tke mixing parameters'
742	WRITE(numout,*) ' coef. to compute avt rn_ediff = ', rn_ediff
743	WRITE(numout,*) ' Kolmogoroff dissipation coef. rn_ediss = ', rn_ediss
744	WRITE(numout,*) ' tke surface input coef. rn_ebb = ', rn_ebb
745	WRITE(numout,*) ' minimum value of tke rn_emin = ', rn_emin
746	WRITE(numout,*) ' surface minimum value of tke rn_emin0 = ', rn_emin0
747	WRITE(numout,*) ' prandl number flag nn_pdl = ', nn_pdl
748	WRITE(numout,*) ' background shear (>0) rn_bshear = ', rn_bshear
749	WRITE(numout,*) ' mixing length type nn_mxl = ', nn_mxl
750	WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0
751	WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0
752	IF( ln_mxl0 ) THEN
753	WRITE(numout,*) ' type of scaling under sea-ice nn_mxlice = ', nn_mxlice
754	IF( nn_mxlice == 1 ) &
755	WRITE(numout,*) ' ice thickness when scaling under sea-ice rn_mxlice = ', rn_mxlice
756	SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice
757	CASE( 0 ) ; WRITE(numout,*) ' ==>>> No scaling under sea-ice'
758	CASE( 1 ) ; WRITE(numout,*) ' ==>>> scaling with constant sea-ice thickness'
759	CASE( 2 ) ; WRITE(numout,*) ' ==>>> scaling with mean sea-ice thickness'
760	CASE( 3 ) ; WRITE(numout,*) ' ==>>> scaling with max sea-ice thickness'
761	CASE DEFAULT
762	CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4')
763	END SELECT
764	ENDIF
765	WRITE(numout,*) ' Langmuir cells parametrization ln_lc = ', ln_lc
766	WRITE(numout,*) ' coef to compute vertical velocity of LC rn_lc = ', rn_lc
767	IF ( cpl_phioc .and. ln_phioc ) THEN
768	SELECT CASE( nn_bc_surf) ! Type of scaling under sea-ice
769	CASE( 0 ) ; WRITE(numout,*) ' nn_bc_surf=0 ==>>> DIRICHLET SBC using surface TKE flux from waves'
770	CASE( 1 ) ; WRITE(numout,*) ' nn_bc_surf=1 ==>>> NEUMANN SBC using surface TKE flux from waves'
771	END SELECT
772	ENDIF
773	WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau
774	WRITE(numout,*) ' type of tke penetration profile nn_htau = ', nn_htau
775	WRITE(numout,*) ' fraction of TKE that penetrates rn_efr = ', rn_efr
776	WRITE(numout,*) ' langmuir & surface wave breaking under ice nn_eice = ', nn_eice
777	SELECT CASE( nn_eice )
778	CASE( 0 ) ; WRITE(numout,*) ' ==>>> no impact of ice cover on langmuir & surface wave breaking'
779	CASE( 1 ) ; WRITE(numout,) ' ==>>> weigthed by 1-TANH( fr_i(:,:) 10 )'
780	CASE( 2 ) ; WRITE(numout,*) ' ==>>> weighted by 1-fr_i(:,:)'
781	CASE( 3 ) ; WRITE(numout,) ' ==>>> weighted by 1-MIN( 1, 4 fr_i(:,:) )'
782	CASE DEFAULT
783	CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3')
784	END SELECT
785	WRITE(numout,*)
786	WRITE(numout,*) ' ==>>> critical Richardson nb with your parameters ri_cri = ', ri_cri
787	WRITE(numout,*)
788	ENDIF
789	!
790	IF( ln_zdfiwm ) THEN ! Internal wave-driven mixing
791	rn_emin = 1.e-10_wp ! specific values of rn_emin & rmxl_min are used
792	rmxl_min = 1.e-03_wp ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s)
793	IF(lwp) WRITE(numout,*) ' ==>>> Internal wave-driven mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3'
794	ELSE ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s)
795	rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) ) ! resulting minimum length to recover molecular viscosity
796	IF(lwp) WRITE(numout,*) ' ==>>> minimum mixing length with your parameters rmxl_min = ', rmxl_min
797	ENDIF
798	!
799	! ! allocate tke arrays
800	IF( zdf_tke_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
801	!
802	! !* Check of some namelist values
803	IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1, 2 or 3' )
804	IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1' )
805	IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0 or 1' )
806	IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
807	!
808	IF( ln_mxl0 ) THEN
809	IF(lwp) WRITE(numout,*)
810	IF(lwp) WRITE(numout,*) ' ==>>> use a surface mixing length = F(stress) : set rn_mxl0 = rmxl_min'
811	rn_mxl0 = rmxl_min
812	ENDIF
813	! !* depth of penetration of surface tke
814	IF( nn_etau /= 0 ) THEN
815	SELECT CASE( nn_htau ) ! Choice of the depth of penetration
816	CASE( 0 ) ! constant depth penetration (here 10 meters)
817	htau(:,:) = 10._wp
818	CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees
819	htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) )
820	END SELECT
821	ENDIF
822	! !* read or initialize all required files
823	CALL tke_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, dissl)
824	!
825	END SUBROUTINE zdf_tke_init
826
827
828	SUBROUTINE tke_rst( kt, cdrw )
829	!!---------------------------------------------------------------------
830	!! * ROUTINE tke_rst *
831	!!
832	!! ** Purpose : Read or write TKE file (en) in restart file
833	!!
834	!! ** Method : use of IOM library
835	!! if the restart does not contain TKE, en is either
836	!! set to rn_emin or recomputed
837	!!----------------------------------------------------------------------
838	USE zdf_oce , ONLY : en, avt_k, avm_k ! ocean vertical physics
839	!!
840	INTEGER , INTENT(in) :: kt ! ocean time-step
841	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
842	!
843	INTEGER :: jit, jk ! dummy loop indices
844	INTEGER :: id1, id2, id3, id4 ! local integers
845	!!----------------------------------------------------------------------
846	!
847	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
848	! ! ---------------
849	IF( ln_rstart ) THEN !* Read the restart file
850	id1 = iom_varid( numror, 'en' , ldstop = .FALSE. )
851	id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
852	id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
853	id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
854	!
855	IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! fields exist
856	CALL iom_get( numror, jpdom_auto, 'en' , en )
857	CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k )
858	CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k )
859	CALL iom_get( numror, jpdom_auto, 'dissl', dissl )
860	ELSE ! start TKE from rest
861	IF(lwp) WRITE(numout,*)
862	IF(lwp) WRITE(numout,*) ' ==>>> previous run without TKE scheme, set en to background values'
863	en (:,:,:) = rn_emin * wmask(:,:,:)
864	dissl(:,:,:) = 1.e-12_wp
865	! avt_k, avm_k already set to the background value in zdf_phy_init
866	ENDIF
867	ELSE !* Start from rest
868	IF(lwp) WRITE(numout,*)
869	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set en to the background value'
870	en (:,:,:) = rn_emin * wmask(:,:,:)
871	dissl(:,:,:) = 1.e-12_wp
872	! avt_k, avm_k already set to the background value in zdf_phy_init
873	ENDIF
874	!
875	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
876	! ! -------------------
877	IF(lwp) WRITE(numout,*) '---- tke_rst ----'
878	CALL iom_rstput( kt, nitrst, numrow, 'en' , en )
879	CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
880	CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k )
881	CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
882	!
883	ENDIF
884	!
885	END SUBROUTINE tke_rst
886
887	!!======================================================================
888	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: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/ZDF/zdftke.F90 @ 15540

Download in other formats: