Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

zdftke.F90

Last change on this file was 85, checked in by smasson, 10 years ago
agrif fixes
File size: 49.0 KB

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

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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

Download in other formats: