Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

zdftke.F90

Last change on this file was 15071, checked in by clem, 3 years ago
make tke and gls work with tiling in debug mode. But it needs to be checked over
Property svn:keywords set to `Id`
File size: 50.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	!! 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(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
95	REAL(wp) :: rmxl_min ! minimum mixing length value (deduced from rn_ediff and rn_emin values) [m]
96	REAL(wp) :: rhftau_add = 1.e-3_wp ! add offset applied to HF part of taum (nn_etau=3)
97	REAL(wp) :: 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(wp), 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(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: p_sh2 ! shear production term
173	REAL(wp), 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(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(in ) :: p_sh2 ! shear production term
205	REAL(wp), 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 , zind ! - -
216	REAL(wp) :: zzd_up, zzd_lw ! - -
217	REAL(wp) :: ztaui, ztauj, z1_norm
218	INTEGER , DIMENSION(A2D(nn_hls)) :: imlc
219	REAL(wp), DIMENSION(A2D(nn_hls)) :: zice_fra, zhlc, zus3, zWlc2
220	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpelc, zdiag, zd_up, zd_lw
221	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztmp ! for diags
222	!!--------------------------------------------------------------------
223	!
224	zbbrau = rn_ebb / rho0 ! Local constant initialisation
225	zbbirau = 3.75_wp / rho0
226	zfact1 = -.5_wp * rn_Dt
227	zfact2 = 1.5_wp * rn_Dt * rn_ediss
228	zfact3 = 0.5_wp * rn_ediss
229	!
230	zpelc(:,:,:) = 0._wp ! need to be initialised in case ln_lc is not used
231	!
232	! ice fraction considered for attenuation of langmuir & wave breaking
233	SELECT CASE ( nn_eice )
234	CASE( 0 ) ; zice_fra(:,:) = 0._wp
235	CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i(A2D(nn_hls)) * 10._wp )
236	CASE( 2 ) ; zice_fra(:,:) = fr_i(A2D(nn_hls))
237	CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp )
238	END SELECT
239	!
240	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
241	! ! Surface/top/bottom boundary condition on tke
242	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
243	!
244	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
245	en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) )
246	zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)
247	zd_lw(ji,jj,1) = 1._wp
248	zd_up(ji,jj,1) = 0._wp
249	END_2D
250	!
251	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
252	! ! Bottom boundary condition on tke
253	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
254	!
255	! en(bot) = (ebb0/rho0)0.5sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
256	! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75)
257	! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2
258	!
259	IF( .NOT.ln_drg_OFF ) THEN !== friction used as top/bottom boundary condition on TKE
260	!
261	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! bottom friction
262	zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
263	zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
264	! ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
265	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 &
266	& + ( zmskv( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )*2 )
267	en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj)
268	END_2D
269	IF( ln_isfcav ) THEN
270	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! top friction
271	zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
272	zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
273	! ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
274	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 &
275	& + ( zmskv( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )*2 )
276	! (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) = 1 where ice shelves are present
277	en(ji,jj,mikt(ji,jj)) = en(ji,jj,1) * tmask(ji,jj,1) &
278	& + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj)
279	END_2D
280	ENDIF
281	!
282	ENDIF
283	!
284	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
285	IF( ln_lc ) THEN ! Langmuir circulation source term added to tke (Axell JGR 2002)
286	! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
287	!
288	! !* Langmuir velocity scale
289	!
290	IF ( cpl_sdrftx ) THEN ! Surface Stokes Drift available
291	! ! Craik-Leibovich velocity scale Wlc = ( u* u_s )^1/2 with u* = (taum/rho0)^1/2
292	! ! associated kinetic energy : 1/2 (Wlc)^2 = u* u_s
293	! ! more precisely, it is the dot product that must be used :
294	! ! 1/2 (W_lc)^2 = MAX( u* u_s + v* v_s , 0 ) only the positive part
295	!!gm ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s
296	!!gm ! so we will overestimate the LC velocity.... !!gm I will do the work if !LC have an effect !
297	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
298	!!XC zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)2 +vt0sd(ji,jj)2 ) )
299	zWlc2(ji,jj) = 0.5_wp * ( ut0sd(ji,jj)2 +vt0sd(ji,jj)2 )
300	END_2D
301	!
302	! Projection of Stokes drift in the wind stress direction
303	!
304	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
305	ztaui = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) )
306	ztauj = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) )
307	z1_norm = 1._wp / MAX( SQRT(ztauiztaui+ztaujztauj), 1.e-12 ) * tmask(ji,jj,1)
308	zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)ztaui + vt0sd(ji,jj)ztauj, 0._wp ) )**2
309	END_2D
310	ELSE ! Surface Stokes drift deduced from surface stress
311	! ! Wlc = u_s with u_s = 0.016*U_10m, the surface stokes drift (Axell 2002, Eq.44)
312	! ! using \|tau\| = rho_air Cd \|U_10m\|^2 , it comes:
313	! ! Wlc = 0.016 * [\|tau\|/(rho_air Cdrag) ]^1/2 and thus:
314	! ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 \|tau\| /( rho_air Cdrag )
315	zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) ! to convert stress in 10m wind using a constant drag
316	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
317	zWlc2(ji,jj) = zcof * taum(ji,jj)
318	END_2D
319	!
320	ENDIF
321	!
322	! !* Depth of the LC circulation (Axell 2002, Eq.47)
323	! !- LHS of Eq.47
324	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
325	zpelc(ji,jj,1) = MAX( rn2b(ji,jj,1), 0._wp ) * gdepw(ji,jj,1,Kmm) * e3w(ji,jj,1,Kmm)
326	END_2D
327	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpk )
328	zpelc(ji,jj,jk) = zpelc(ji,jj,jk-1) + &
329	& MAX( rn2b(ji,jj,jk), 0._wp ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
330	END_3D
331	!
332	! !- compare LHS to RHS of Eq.47
333	imlc(:,:) = mbkt(A2D(nn_hls)) + 1 ! Initialization to the number of w ocean point (=2 over land)
334	DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )
335	IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) ) imlc(ji,jj) = jk
336	END_3D
337	! ! finite LC depth
338	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
339	zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm)
340	END_2D
341	!
342	zcof = 0.016 / SQRT( zrhoa * zcdrag )
343	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
344	zus = SQRT( 2. * zWlc2(ji,jj) ) ! Stokes drift
345	zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
346	END_2D
347	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
348	IF ( zus3(ji,jj) /= 0._wp ) THEN
349	IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN
350	! ! vertical velocity due to LC
351	zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )
352	! ! TKE Langmuir circulation source term
353	en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
354	ENDIF
355	ENDIF
356	END_3D
357	!
358	ENDIF
359	!
360	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
361	! ! Now Turbulent kinetic energy (output in en)
362	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
363	! ! Resolution of a tridiagonal linear system by a "methode de chasse"
364	! ! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ).
365	! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
366	!
367	IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri )
368	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
369	! ! local Richardson number
370	IF (rn2b(ji,jj,jk) <= 0.0_wp) then
371	zri = 0.0_wp
372	ELSE
373	zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
374	ENDIF
375	! ! inverse of Prandtl number
376	apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) )
377	END_3D
378	ENDIF
379	!
380	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* Matrix and right hand side in en
381	zcof = zfact1 * tmask(ji,jj,jk)
382	! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
383	! ! eddy coefficient (ensure numerical stability)
384	zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal
385	& / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk ,Kmm) )
386	zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal
387	& / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk ,Kmm) )
388	!
389	zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw)
390	zd_lw(ji,jj,jk) = zzd_lw
391	zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk)
392	!
393	! ! right hand side in en
394	en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * ( p_sh2(ji,jj,jk) & ! shear
395	& - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification
396	& + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation
397	& ) * wmask(ji,jj,jk)
398	END_3D
399	!
400	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
401	! ! Surface boundary condition on tke if
402	! ! coupling with waves
403	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
404	!
405	IF ( cpl_phioc .and. ln_phioc ) THEN
406	SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves
407
408	CASE ( 0 ) ! Dirichlet BC
409	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)
410	IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp
411	en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )*(2./3.) ) tmask(ji,jj,1)
412	zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) ! choose to keep coherence with former estimation of
413	END_2D
414
415	CASE ( 1 ) ! Neumann BC
416	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
417	IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp
418	en(ji,jj,2) = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm)
419	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) )
420	zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2)
421	zdiag(ji,jj,1) = 1._wp
422	zd_lw(ji,jj,2) = 0._wp
423	END_2D
424
425	END SELECT
426
427	ENDIF
428	!
429	! !* Matrix inversion from level 2 (tke prescribed at level 1)
430	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
431	zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
432	END_3D
433	!XC : commented to allow for neumann boundary condition
434	! DO_2D( 0, 0, 0, 0 )
435	! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke
436	! END_2D
437	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
438	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)
439	END_3D
440	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
441	en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
442	END_2D
443	DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 )
444	en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
445	END_3D
446	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! set the minimum value of tke
447	en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
448	END_3D
449	!
450	! Kolmogorov energy of dissipation (W/kg)
451	! ediss = Cesqrt(en)/Len
452	! dissl = sqrt(en)/L
453	IF( iom_use('ediss_k') ) THEN
454	ALLOCATE( ztmp(A2D(nn_hls),jpk) )
455	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
456	ztmp(ji,jj,jk) = zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) * wmask(ji,jj,jk)
457	END_3D
458	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
459	ztmp(ji,jj,jpk) = 0._wp
460	END_2D
461	CALL iom_put( 'ediss_k', ztmp )
462	DEALLOCATE( ztmp )
463	ENDIF
464	!
465	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
466	! ! TKE due to surface and internal wave breaking
467	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
468	!!gm BUG : in the exp remove the depth of ssh !!!
469	!!gm i.e. use gde3w in argument (gdepw(:,:,:,Kmm))
470	!
471	! penetration is partly switched off below sea-ice if nn_eice/=0
472	!
473	IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction)
474	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
475	en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) &
476	& * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
477	END_3D
478	ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction)
479	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
480	jk = nmln(ji,jj)
481	en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) &
482	& * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
483	END_2D
484	ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability)
485	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
486	ztx2 = utau(ji-1,jj ) + utau(ji,jj)
487	zty2 = vtau(ji ,jj-1) + vtau(ji,jj)
488	ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress
489	zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean
490	zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications...
491	en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) &
492	& * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
493	END_3D
494	ENDIF
495	!
496	END SUBROUTINE tke_tke
497
498
499	SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt )
500	!!----------------------------------------------------------------------
501	!! * ROUTINE tke_avn *
502	!!
503	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
504	!!
505	!! ** Method : At this stage, en, the now TKE, is known (computed in
506	!! the tke_tke routine). First, the now mixing lenth is
507	!! computed from en and the strafification (N^2), then the mixings
508	!! coefficients are computed.
509	!! - Mixing length : a first evaluation of the mixing lengh
510	!! scales is:
511	!! mxl = sqrt(2*en) / N
512	!! where N is the brunt-vaisala frequency, with a minimum value set
513	!! to rmxl_min (rn_mxl0) in the interior (surface) ocean.
514	!! The mixing and dissipative length scale are bound as follow :
515	!! nn_mxl=0 : mxl bounded by the distance to surface and bottom.
516	!! zmxld = zmxlm = mxl
517	!! nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
518	!! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
519	!! less than 1 (\|d/dz(mxl)\|<1) and zmxld = zmxlm = mxl
520	!! nn_mxl=3 : mxl is bounded from the surface to the bottom usings
521	!! \|d/dz(xml)\|<1 to obtain lup, and from the bottom to
522	!! the surface to obtain ldown. the resulting length
523	!! scales are:
524	!! zmxld = sqrt( lup * ldown )
525	!! zmxlm = min ( lup , ldown )
526	!! - Vertical eddy viscosity and diffusivity:
527	!! avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
528	!! avt = max( avmb, pdlr * avm )
529	!! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
530	!!
531	!! ** Action : - avt, avm : now vertical eddy diffusivity and viscosity (w-point)
532	!!----------------------------------------------------------------------
533	USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d ! ocean vertical physics
534	!!
535	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
536	REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points)
537	!
538	INTEGER :: ji, jj, jk ! dummy loop indices
539	REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars
540	REAL(wp) :: zdku, zdkv, zsqen ! - -
541	REAL(wp) :: zemxl, zemlm, zemlp, zmaxice ! - -
542	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zmxlm, zmxld ! 3D workspace
543	!!--------------------------------------------------------------------
544	!
545	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
546	! ! Mixing length
547	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
548	!
549	! !* Buoyancy length scale: l=sqrt(2e/n*2)
550	!
551	! initialisation of interior minimum value (avoid a 2d loop with mikt)
552	zmxlm(:,:,:) = rmxl_min
553	zmxld(:,:,:) = rmxl_min
554	!
555	IF(ln_sdw .AND. ln_mxhsw) THEN
556	zmxlm(:,:,1)= vkarmn * MAX ( 1.6 * hsw(:,:) , 0.02 ) ! surface mixing length = F(wave height)
557	! from terray et al 1999 and mellor and blumberg 2004 it should be 0.85 and not 1.6
558	zcoef = vkarmn * ( (rn_ediffrn_ediss)*0.25 ) / rn_ediff
559	zmxlm(:,:,1)= zcoef * MAX ( 1.6 * hsw(:,:) , 0.02 ) ! surface mixing length = F(wave height)
560	ELSE
561	!
562	IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn2.e5taum/(rho0*g)
563	!
564	zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
565	#if ! defined key_si3 && ! defined key_cice
566	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! No sea-ice
567	zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
568	END_2D
569	#else
570	SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice
571	!
572	CASE( 0 ) ! No scaling under sea-ice
573	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
574	zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
575	END_2D
576	!
577	CASE( 1 ) ! scaling with constant sea-ice thickness
578	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
579	zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
580	& fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1)
581	END_2D
582	!
583	CASE( 2 ) ! scaling with mean sea-ice thickness
584	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
585	#if defined key_si3
586	zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
587	& fr_i(ji,jj) * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1)
588	#elif defined key_cice
589	zmaxice = MAXVAL( h_i(ji,jj,:) )
590	zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
591	& fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1)
592	#endif
593	END_2D
594	!
595	CASE( 3 ) ! scaling with max sea-ice thickness
596	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
597	zmaxice = MAXVAL( h_i(ji,jj,:) )
598	zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
599	& fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1)
600	END_2D
601	!
602	END SELECT
603	#endif
604	!
605	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
606	zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
607	END_2D
608	!
609	ELSE
610	zmxlm(:,:,1) = rn_mxl0
611	ENDIF
612	ENDIF
613	!
614	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
615	zrn2 = MAX( rn2(ji,jj,jk), rsmall )
616	zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) )
617	END_3D
618	!
619	! !* Physical limits for the mixing length
620	!
621	zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the minimum value
622	zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value
623	!
624	SELECT CASE ( nn_mxl )
625	!
626	!!gm Not sure of that coding for ISF....
627	! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm)
628	CASE ( 0 ) ! bounded by the distance to surface and bottom
629	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
630	zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk), &
631	& gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) )
632	! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
633	zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) &
634	& + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
635	zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) &
636	& + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
637	END_3D
638	!
639	CASE ( 1 ) ! bounded by the vertical scale factor
640	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
641	zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) )
642	zmxlm(ji,jj,jk) = zemxl
643	zmxld(ji,jj,jk) = zemxl
644	END_3D
645	!
646	CASE ( 2 ) ! \|dk[xml]\| bounded by e3t :
647	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! from the surface to the bottom :
648	zmxlm(ji,jj,jk) = &
649	& MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
650	END_3D
651	DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! from the bottom to the surface :
652	zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
653	zmxlm(ji,jj,jk) = zemxl
654	zmxld(ji,jj,jk) = zemxl
655	END_3D
656	!
657	CASE ( 3 ) ! lup and ldown, \|dk[xml]\| bounded by e3t :
658	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! from the surface to the bottom : lup
659	zmxld(ji,jj,jk) = &
660	& MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
661	END_3D
662	DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! from the bottom to the surface : ldown
663	zmxlm(ji,jj,jk) = &
664	& MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
665	END_3D
666	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
667	zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) )
668	zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
669	zmxlm(ji,jj,jk) = zemlm
670	zmxld(ji,jj,jk) = zemlp
671	END_3D
672	!
673	END SELECT
674	!
675	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
676	! ! Vertical eddy viscosity and diffusivity (avm and avt)
677	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
678	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* vertical eddy viscosity & diffivity at w-points
679	zsqen = SQRT( en(ji,jj,jk) )
680	zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen
681	p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk)
682	p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
683	dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
684	END_3D
685	!
686	!
687	IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt
688	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
689	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)
690	END_3D
691	ENDIF
692	!
693	IF(sn_cfctl%l_prtctl) THEN
694	CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk)
695	CALL prt_ctl( tab3d_1=p_avm, clinfo1=' tke - m: ', kdim=jpk )
696	ENDIF
697	!
698	END SUBROUTINE tke_avn
699
700
701	SUBROUTINE zdf_tke_init( Kmm )
702	!!----------------------------------------------------------------------
703	!! * ROUTINE zdf_tke_init *
704	!!
705	!! ** Purpose : Initialization of the vertical eddy diffivity and
706	!! viscosity when using a tke turbulent closure scheme
707	!!
708	!! ** Method : Read the namzdf_tke namelist and check the parameters
709	!! called at the first timestep (nit000)
710	!!
711	!! ** input : Namlist namzdf_tke
712	!!
713	!! ** Action : Increase by 1 the nstop flag is setting problem encounter
714	!!----------------------------------------------------------------------
715	USE zdf_oce , ONLY : ln_zdfiwm ! Internal Wave Mixing flag
716	!!
717	INTEGER, INTENT(in) :: Kmm ! time level index
718	INTEGER :: ji, jj, jk ! dummy loop indices
719	INTEGER :: ios
720	!!
721	NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , &
722	& rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , &
723	& rn_mxl0 , nn_mxlice, rn_mxlice, &
724	& nn_pdl , ln_lc , rn_lc , &
725	& nn_etau , nn_htau , rn_efr , nn_eice , &
726	& nn_bc_surf, nn_bc_bot, ln_mxhsw
727	!!----------------------------------------------------------------------
728	!
729	READ ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
730	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist' )
731
732	READ ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
733	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist' )
734	IF(lwm) WRITE ( numond, namzdf_tke )
735	!
736	ri_cri = 2._wp / ( 2._wp + rn_ediss / rn_ediff ) ! resulting critical Richardson number
737	!
738	IF(lwp) THEN !* Control print
739	WRITE(numout,*)
740	WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
741	WRITE(numout,*) '~~~~~~~~~~~~'
742	WRITE(numout,*) ' Namelist namzdf_tke : set tke mixing parameters'
743	WRITE(numout,*) ' coef. to compute avt rn_ediff = ', rn_ediff
744	WRITE(numout,*) ' Kolmogoroff dissipation coef. rn_ediss = ', rn_ediss
745	WRITE(numout,*) ' tke surface input coef. rn_ebb = ', rn_ebb
746	WRITE(numout,*) ' minimum value of tke rn_emin = ', rn_emin
747	WRITE(numout,*) ' surface minimum value of tke rn_emin0 = ', rn_emin0
748	WRITE(numout,*) ' prandl number flag nn_pdl = ', nn_pdl
749	WRITE(numout,*) ' background shear (>0) rn_bshear = ', rn_bshear
750	WRITE(numout,*) ' mixing length type nn_mxl = ', nn_mxl
751	WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0
752	WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0
753	IF( ln_mxl0 ) THEN
754	WRITE(numout,*) ' type of scaling under sea-ice nn_mxlice = ', nn_mxlice
755	IF( nn_mxlice == 1 ) &
756	WRITE(numout,*) ' ice thickness when scaling under sea-ice rn_mxlice = ', rn_mxlice
757	SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice
758	CASE( 0 ) ; WRITE(numout,*) ' ==>>> No scaling under sea-ice'
759	CASE( 1 ) ; WRITE(numout,*) ' ==>>> scaling with constant sea-ice thickness'
760	CASE( 2 ) ; WRITE(numout,*) ' ==>>> scaling with mean sea-ice thickness'
761	CASE( 3 ) ; WRITE(numout,*) ' ==>>> scaling with max sea-ice thickness'
762	CASE DEFAULT
763	CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4')
764	END SELECT
765	ENDIF
766	WRITE(numout,*) ' Langmuir cells parametrization ln_lc = ', ln_lc
767	WRITE(numout,*) ' coef to compute vertical velocity of LC rn_lc = ', rn_lc
768	IF ( cpl_phioc .and. ln_phioc ) THEN
769	SELECT CASE( nn_bc_surf) ! Type of scaling under sea-ice
770	CASE( 0 ) ; WRITE(numout,*) ' nn_bc_surf=0 ==>>> DIRICHLET SBC using surface TKE flux from waves'
771	CASE( 1 ) ; WRITE(numout,*) ' nn_bc_surf=1 ==>>> NEUMANN SBC using surface TKE flux from waves'
772	END SELECT
773	ENDIF
774	WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau
775	WRITE(numout,*) ' type of tke penetration profile nn_htau = ', nn_htau
776	WRITE(numout,*) ' fraction of TKE that penetrates rn_efr = ', rn_efr
777	WRITE(numout,*) ' langmuir & surface wave breaking under ice nn_eice = ', nn_eice
778	SELECT CASE( nn_eice )
779	CASE( 0 ) ; WRITE(numout,*) ' ==>>> no impact of ice cover on langmuir & surface wave breaking'
780	CASE( 1 ) ; WRITE(numout,) ' ==>>> weigthed by 1-TANH( fr_i(:,:) 10 )'
781	CASE( 2 ) ; WRITE(numout,*) ' ==>>> weighted by 1-fr_i(:,:)'
782	CASE( 3 ) ; WRITE(numout,) ' ==>>> weighted by 1-MIN( 1, 4 fr_i(:,:) )'
783	CASE DEFAULT
784	CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3')
785	END SELECT
786	WRITE(numout,*)
787	WRITE(numout,*) ' ==>>> critical Richardson nb with your parameters ri_cri = ', ri_cri
788	WRITE(numout,*)
789	ENDIF
790	!
791	IF( ln_zdfiwm ) THEN ! Internal wave-driven mixing
792	rn_emin = 1.e-10_wp ! specific values of rn_emin & rmxl_min are used
793	rmxl_min = 1.e-03_wp ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s)
794	IF(lwp) WRITE(numout,*) ' ==>>> Internal wave-driven mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3'
795	ELSE ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s)
796	rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) ) ! resulting minimum length to recover molecular viscosity
797	IF(lwp) WRITE(numout,*) ' ==>>> minimum mixing length with your parameters rmxl_min = ', rmxl_min
798	ENDIF
799	!
800	! ! allocate tke arrays
801	IF( zdf_tke_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
802	!
803	! !* Check of some namelist values
804	IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1, 2 or 3' )
805	IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1' )
806	IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0 or 1' )
807	IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
808	!
809	IF( ln_mxl0 ) THEN
810	IF(lwp) WRITE(numout,*)
811	IF(lwp) WRITE(numout,*) ' ==>>> use a surface mixing length = F(stress) : set rn_mxl0 = rmxl_min'
812	rn_mxl0 = rmxl_min
813	ENDIF
814	! !* depth of penetration of surface tke
815	IF( nn_etau /= 0 ) THEN
816	SELECT CASE( nn_htau ) ! Choice of the depth of penetration
817	CASE( 0 ) ! constant depth penetration (here 10 meters)
818	htau(:,:) = 10._wp
819	CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees
820	htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) )
821	END SELECT
822	ENDIF
823	! !* read or initialize all required files
824	CALL tke_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, dissl)
825	!
826	END SUBROUTINE zdf_tke_init
827
828
829	SUBROUTINE tke_rst( kt, cdrw )
830	!!---------------------------------------------------------------------
831	!! * ROUTINE tke_rst *
832	!!
833	!! ** Purpose : Read or write TKE file (en) in restart file
834	!!
835	!! ** Method : use of IOM library
836	!! if the restart does not contain TKE, en is either
837	!! set to rn_emin or recomputed
838	!!----------------------------------------------------------------------
839	USE zdf_oce , ONLY : en, avt_k, avm_k ! ocean vertical physics
840	!!
841	INTEGER , INTENT(in) :: kt ! ocean time-step
842	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
843	!
844	INTEGER :: jit, jk ! dummy loop indices
845	INTEGER :: id1, id2, id3, id4 ! local integers
846	!!----------------------------------------------------------------------
847	!
848	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
849	! ! ---------------
850	IF( ln_rstart ) THEN !* Read the restart file
851	id1 = iom_varid( numror, 'en' , ldstop = .FALSE. )
852	id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
853	id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
854	id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
855	!
856	IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! fields exist
857	CALL iom_get( numror, jpdom_auto, 'en' , en )
858	CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k )
859	CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k )
860	CALL iom_get( numror, jpdom_auto, 'dissl', dissl )
861	ELSE ! start TKE from rest
862	IF(lwp) WRITE(numout,*)
863	IF(lwp) WRITE(numout,*) ' ==>>> previous run without TKE scheme, set en to background values'
864	en (:,:,:) = rn_emin * wmask(:,:,:)
865	dissl(:,:,:) = 1.e-12_wp
866	! avt_k, avm_k already set to the background value in zdf_phy_init
867	ENDIF
868	ELSE !* Start from rest
869	IF(lwp) WRITE(numout,*)
870	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set en to the background value'
871	en (:,:,:) = rn_emin * wmask(:,:,:)
872	dissl(:,:,:) = 1.e-12_wp
873	! avt_k, avm_k already set to the background value in zdf_phy_init
874	ENDIF
875	!
876	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
877	! ! -------------------
878	IF(lwp) WRITE(numout,*) '---- tke_rst ----'
879	CALL iom_rstput( kt, nitrst, numrow, 'en' , en )
880	CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
881	CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k )
882	CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
883	!
884	ENDIF
885	!
886	END SUBROUTINE tke_rst
887
888	!!======================================================================
889	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/trunk/src/OCE/ZDF/zdftke.F90

Download in other formats: