Context Navigation

zdfgls.F90 @ 8279

Last change on this file since 8279 was 8279, checked in by mocavero, 7 years ago
Implementation of OMP coarse-grained parallelization on ZDF new package
Property svn:keywords set to `Id`
File size: 54.4 KB

Line
1	MODULE zdfgls
2	!!======================================================================
3	!! * MODULE zdfgls *
4	!! Ocean physics: vertical mixing coefficient computed from the gls
5	!! turbulent closure parameterization
6	!!======================================================================
7	!! History : 3.0 ! 2009-09 (G. Reffray) Original code
8	!! 3.3 ! 2010-10 (C. Bricaud) Add in the reference
9	!! 4.0 ! 2017-04 (G. Madec) remove CPP keys & avm at t-point only
10	!!----------------------------------------------------------------------
11
12	!!----------------------------------------------------------------------
13	!! zdf_gls : update momentum and tracer Kz from a gls scheme
14	!! zdf_gls_init : initialization, namelist read, and parameters control
15	!! gls_rst : read/write gls restart in ocean restart file
16	!!----------------------------------------------------------------------
17	USE oce ! ocean dynamics and active tracers
18	USE dom_oce ! ocean space and time domain
19	USE domvvl ! ocean space and time domain : variable volume layer
20	USE zdf_oce ! ocean vertical physics
21	USE zdfbfr ! bottom friction (only for rn_bfrz0)
22	USE sbc_oce ! surface boundary condition: ocean
23	USE phycst ! physical constants
24	USE zdfmxl ! mixed layer
25	USE sbcwave , ONLY: hsw ! significant wave height
26	!
27	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
28	USE lib_mpp ! MPP manager
29	USE wrk_nemo ! work arrays
30	USE prtctl ! Print control
31	USE in_out_manager ! I/O manager
32	USE iom ! I/O manager library
33	USE timing ! Timing
34	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
35
36	IMPLICIT NONE
37	PRIVATE
38
39	PUBLIC zdf_gls ! called in zdfphy
40	PUBLIC zdf_gls_init ! called in zdfphy
41	PUBLIC gls_rst ! called in zdfphy
42
43	!
44	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hmxl_n !: now mixing length
45	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zwall !: wall function
46	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points
47	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustarb2 !: Squared bottom velocity scale at T-points
48
49	! !! Namelist namzdf_gls
50	LOGICAL :: ln_length_lim ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988)
51	LOGICAL :: ln_sigpsi ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing
52	INTEGER :: nn_bc_surf ! surface boundary condition (=0/1)
53	INTEGER :: nn_bc_bot ! bottom boundary condition (=0/1)
54	INTEGER :: nn_z0_met ! Method for surface roughness computation
55	INTEGER :: nn_stab_func ! stability functions G88, KC or Canuto (=0/1/2)
56	INTEGER :: nn_clos ! closure 0/1/2/3 MY82/k-eps/k-w/gen
57	REAL(wp) :: rn_clim_galp ! Holt 2008 value for k-eps: 0.267
58	REAL(wp) :: rn_epsmin ! minimum value of dissipation (m2/s3)
59	REAL(wp) :: rn_emin ! minimum value of TKE (m2/s2)
60	REAL(wp) :: rn_charn ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value)
61	REAL(wp) :: rn_crban ! Craig and Banner constant for surface breaking waves mixing
62	REAL(wp) :: rn_hsro ! Minimum surface roughness
63	REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met > 1)
64
65	REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters
66	REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1
67	REAL(wp) :: rl_sf = 0.2_wp ! 0 <rl_sf<vkarmn
68	REAL(wp) :: rghmin = -0.28_wp
69	REAL(wp) :: rgh0 = 0.0329_wp
70	REAL(wp) :: rghcri = 0.03_wp
71	REAL(wp) :: ra1 = 0.92_wp
72	REAL(wp) :: ra2 = 0.74_wp
73	REAL(wp) :: rb1 = 16.60_wp
74	REAL(wp) :: rb2 = 10.10_wp
75	REAL(wp) :: re2 = 1.33_wp
76	REAL(wp) :: rl1 = 0.107_wp
77	REAL(wp) :: rl2 = 0.0032_wp
78	REAL(wp) :: rl3 = 0.0864_wp
79	REAL(wp) :: rl4 = 0.12_wp
80	REAL(wp) :: rl5 = 11.9_wp
81	REAL(wp) :: rl6 = 0.4_wp
82	REAL(wp) :: rl7 = 0.0_wp
83	REAL(wp) :: rl8 = 0.48_wp
84	REAL(wp) :: rm1 = 0.127_wp
85	REAL(wp) :: rm2 = 0.00336_wp
86	REAL(wp) :: rm3 = 0.0906_wp
87	REAL(wp) :: rm4 = 0.101_wp
88	REAL(wp) :: rm5 = 11.2_wp
89	REAL(wp) :: rm6 = 0.4_wp
90	REAL(wp) :: rm7 = 0.0_wp
91	REAL(wp) :: rm8 = 0.318_wp
92	REAL(wp) :: rtrans = 0.1_wp
93	REAL(wp) :: rc02, rc02r, rc03, rc04 ! coefficients deduced from above parameters
94	REAL(wp) :: rsbc_tke1, rsbc_tke2, rfact_tke ! - - - -
95	REAL(wp) :: rsbc_psi1, rsbc_psi2, rfact_psi ! - - - -
96	REAL(wp) :: rsbc_zs1, rsbc_zs2 ! - - - -
97	REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff ! - - - -
98	REAL(wp) :: rs0, rs1, rs2, rs4, rs5, rs6 ! - - - -
99	REAL(wp) :: rd0, rd1, rd2, rd3, rd4, rd5 ! - - - -
100	REAL(wp) :: rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0 ! - - - -
101	REAL(wp) :: rpsi3m, rpsi3p, rpp, rmm, rnn ! - - - -
102
103	!! * Substitutions
104	# include "domain_substitute.h90"
105	!!----------------------------------------------------------------------
106	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
107	!! $Id$
108	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
109	!!----------------------------------------------------------------------
110	CONTAINS
111
112	INTEGER FUNCTION zdf_gls_alloc()
113	!!----------------------------------------------------------------------
114	!! * FUNCTION zdf_gls_alloc *
115	!!----------------------------------------------------------------------
116	ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustars2(jpi,jpj) , &
117	& zwall (jpi,jpj,jpk) , ustarb2(jpi,jpj) , STAT= zdf_gls_alloc )
118	!
119	IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc )
120	IF( zdf_gls_alloc /= 0 ) CALL ctl_warn('zdf_gls_alloc: failed to allocate arrays')
121	END FUNCTION zdf_gls_alloc
122
123
124	SUBROUTINE zdf_gls( ARG_2D, kt, psh2, p_avm, p_avt )
125	!!----------------------------------------------------------------------
126	!! * ROUTINE zdf_gls *
127	!!
128	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
129	!! coefficients using the GLS turbulent closure scheme.
130	!!----------------------------------------------------------------------
131	INTEGER , INTENT(in ) :: ARG_2D ! inner domain start-end i-indices
132	INTEGER , INTENT(in ) :: kt ! ocean time step
133	REAL(wp), DIMENSION(WRK_3D), INTENT(in ) :: psh2 ! shear production term
134	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)
135	!
136	INTEGER :: ji, jj, jk, ibot, ibotm1, dir ! dummy loop arguments
137	REAL(wp) :: zesh2, zsigpsi, zcoef, zex1, zex2 ! local scalars
138	REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - -
139	REAL(wp) :: zratio, zrn2, zflxb, sh ! - -
140	REAL(wp) :: prod, buoy, diss, zdiss, sm ! - -
141	REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - -
142	REAL(wp), DIMENSION(WRK_2D) :: zdep
143	REAL(wp), DIMENSION(WRK_2D) :: zkar
144	REAL(wp), DIMENSION(WRK_2D) :: zflxs ! Turbulence fluxed induced by internal waves
145	REAL(wp), DIMENSION(WRK_2D) :: zhsro ! Surface roughness (surface waves)
146	REAL(wp), DIMENSION(WRK_3D) :: eb ! tke at time before
147	REAL(wp), DIMENSION(WRK_3D) :: hmxl_b ! mixing length at time before
148	REAL(wp), DIMENSION(WRK_3D) :: eps ! dissipation rate
149	REAL(wp), DIMENSION(WRK_3D) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi)
150	REAL(wp), DIMENSION(WRK_3D) :: psi ! psi at time now
151	REAL(wp), DIMENSION(WRK_3D) :: z_elem_a ! element of the first matrix diagonal
152	REAL(wp), DIMENSION(WRK_3D) :: z_elem_b ! element of the second matrix diagonal
153	REAL(wp), DIMENSION(WRK_3D) :: z_elem_c ! element of the third matrix diagonal
154	REAL(wp), DIMENSION(WRK_3D) :: zstt, zstm ! stability function on tracer and momentum
155	!!--------------------------------------------------------------------
156	!
157	IF( nn_timing == 1 ) CALL timing_start('zdf_gls')
158	!
159	! Preliminary computing
160
161	ustars2(WRK_2D) = 0._wp ; ustarb2(WRK_2D) = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp
162
163
164	! Compute surface and bottom friction at T-points
165	DO jj = tnldj, tnlej
166	DO ji = tnldi, tnlei
167	!
168	! surface friction
169	ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
170	!
171	! bottom friction (explicit before friction)
172	! Note that we chose here not to bound the friction as in dynbfr)
173	ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) &
174	& * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) )
175	zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) &
176	& * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) )
177	ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
178	END DO
179	END DO
180
181	SELECT CASE ( nn_z0_met ) !== Set surface roughness length ==!
182	CASE ( 0 ) ! Constant roughness
183	zhsro(WRK_2D) = rn_hsro
184	CASE ( 1 ) ! Standard Charnock formula
185	zhsro(WRK_2D) = MAX(rsbc_zs1 * ustars2(WRK_2D), rn_hsro)
186	CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008)
187	!!gm zcof = 2._wp * 0.6_wp / 28._wp
188	!!gm zdep(WRK_2D) = 30._wp * TANH( zcof/ SQRT( MAX(ustars2(WRK_2D),rsmall) ) ) ! Wave age (eq. 10)
189	zdep(WRK_2D) = 30.TANH(2.0.3/(28.*SQRT(MAX(ustars2(WRK_2D),rsmall)))) ! Wave age (eq. 10)
190	zhsro(WRK_2D) = MAX(rsbc_zs2 * ustars2(WRK_2D) * zdep(WRK_2D)*1.5, rn_hsro) ! zhsro = rn_frac_hs Hsw (eq. 11)
191	CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file)
192	!!gm BUG missing a multiplicative coefficient....
193	zhsro(WRK_2D) = hsw(WRK_2D)
194	END SELECT
195	!
196	DO jk = 2, jpkm1 !== Compute dissipation rate ==!
197	DO jj = tnldj, tnlej
198	DO ji = tnldi, tnlei
199	eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk)
200	END DO
201	END DO
202	END DO
203
204	! Save tke at before time step
205	eb (WRK_2D,:) = en (WRK_2D,:)
206	hmxl_b(WRK_2D,:) = hmxl_n(WRK_2D,:)
207
208	IF( nn_clos == 0 ) THEN ! Mellor-Yamada
209	DO jk = 2, jpkm1
210	DO jj = tnldj, tnlej
211	DO ji = tnldi, tnlei
212	zup = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)
213	zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) )
214	zcoef = ( zup / MAX( zdown, rsmall ) )
215	zwall (ji,jj,jk) = ( 1._wp + re2 * zcoefzcoef ) tmask(ji,jj,jk)
216	END DO
217	END DO
218	END DO
219	ENDIF
220
221	!!---------------------------------!!
222	!! Equation to prognostic k !!
223	!!---------------------------------!!
224	!
225	! Now Turbulent kinetic energy (output in en)
226	! -------------------------------
227	! Resolution of a tridiagonal linear system by a "methode de chasse"
228	! computation from level 2 to jpkm1 (e(1) computed after and e(jpk)=0 ).
229	! The surface boundary condition are set after
230	! The bottom boundary condition are also set after. In standard e(bottom)=0.
231	! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal
232	! Warning : after this step, en : right hand side of the matrix
233
234	DO jk = 2, jpkm1
235	DO jj = tnldj, tnlej
236	DO ji = tnldi, tnlei
237	!
238	buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction
239	!
240	diss = eps(ji,jj,jk) ! dissipation
241	!
242	dir = 0.5_wp + SIGN( 0.5_wp, psh2(ji,jj,jk) + buoy ) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
243	!
244	zesh2 = dir(psh2(ji,jj,jk)+buoy)+(1._wp-dir)psh2(ji,jj,jk) ! production term
245	zdiss = dir(diss/en(ji,jj,jk)) +(1._wp-dir)(diss-buoy)/en(ji,jj,jk) ! dissipation term
246	!
247	! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0
248	! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style)
249	! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries.
250	! Otherwise, this should be rsc_psi/rsc_psi0
251	IF( ln_sigpsi ) THEN
252	zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1.
253	zwall_psi(ji,jj,jk) = rsc_psi / &
254	& ( zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) )
255	ELSE
256	zwall_psi(ji,jj,jk) = 1._wp
257	ENDIF
258	!
259	! building the matrix
260	zcof = rfact_tke * tmask(ji,jj,jk)
261	! ! lower diagonal
262	z_elem_a(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
263	! ! upper diagonal
264	z_elem_c(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) )
265	! ! diagonal
266	z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk)
267	! ! right hand side in en
268	en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk)
269	END DO
270	END DO
271	END DO
272	!
273	z_elem_b(WRK_2D,jpk) = 1._wp
274	!
275	! Set surface condition on zwall_psi (1 at the bottom)
276	zwall_psi(WRK_2D, 1 ) = zwall_psi(WRK_2D,2)
277	zwall_psi(WRK_2D,jpk) = 1._wp
278	!
279	! Surface boundary condition on tke
280	! ---------------------------------
281	!
282	SELECT CASE ( nn_bc_surf )
283	!
284	CASE ( 0 ) ! Dirichlet case
285	! First level
286	en(WRK_2D,1) = rc02r * ustars2(WRK_2D) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
287	en(WRK_2D,1) = MAX(en(WRK_2D,1), rn_emin)
288	z_elem_a(WRK_2D,1) = en(WRK_2D,1)
289	z_elem_c(WRK_2D,1) = 0._wp
290	z_elem_b(WRK_2D,1) = 1._wp
291	!
292	! One level below
293	en(WRK_2D,2) = rc02r * ustars2(WRK_2D) * (1._wp + rsbc_tke1 * ((zhsro(WRK_2D)+gdepw_n(WRK_2D,2)) &
294	& / zhsro(WRK_2D) )*(1.5_wpra_sf))**(2._wp/3._wp)
295	en(WRK_2D,2) = MAX(en(WRK_2D,2), rn_emin )
296	z_elem_a(WRK_2D,2) = 0._wp
297	z_elem_c(WRK_2D,2) = 0._wp
298	z_elem_b(WRK_2D,2) = 1._wp
299	!
300	!
301	CASE ( 1 ) ! Neumann boundary condition on d(e)/dz
302	!
303	! Dirichlet conditions at k=1
304	en(WRK_2D,1) = rc02r * ustars2(WRK_2D) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
305	en(WRK_2D,1) = MAX(en(WRK_2D,1), rn_emin)
306	z_elem_a(WRK_2D,1) = en(WRK_2D,1)
307	z_elem_c(WRK_2D,1) = 0._wp
308	z_elem_b(WRK_2D,1) = 1._wp
309	!
310	! at k=2, set de/dz=Fw
311	!cbr
312	z_elem_b(WRK_2D,2) = z_elem_b(WRK_2D,2) + z_elem_a(WRK_2D,2) ! Remove z_elem_a from z_elem_b
313	z_elem_a(WRK_2D,2) = 0._wp
314	zkar(WRK_2D) = (rl_sf + (vkarmn-rl_sf)(1.-exp(-rtransgdept_n(WRK_2D,1)/zhsro(WRK_2D)) ))
315	zflxs(WRK_2D) = rsbc_tke2 * ustars2(WRK_2D)*1.5_wp zkar(WRK_2D) &
316	& * ((zhsro(WRK_2D)+gdept_n(WRK_2D,1)) / zhsro(WRK_2D) )*(1.5_wpra_sf)
317
318	en(WRK_2D,2) = en(WRK_2D,2) + zflxs(WRK_2D)/e3w_n(WRK_2D,2)
319	!
320	!
321	END SELECT
322
323	! Bottom boundary condition on tke
324	! --------------------------------
325	!
326	SELECT CASE ( nn_bc_bot )
327	!
328	CASE ( 0 ) ! Dirichlet
329	! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
330	! ! Balance between the production and the dissipation terms
331	DO jj = tnldj, tnlej
332	DO ji = tnldi, tnlei
333	ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point
334	ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1
335	!
336	! Bottom level Dirichlet condition:
337	z_elem_a(ji,jj,ibot ) = 0._wp
338	z_elem_c(ji,jj,ibot ) = 0._wp
339	z_elem_b(ji,jj,ibot ) = 1._wp
340	en(ji,jj,ibot ) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
341	!
342	! Just above last level, Dirichlet condition again
343	z_elem_a(ji,jj,ibotm1) = 0._wp
344	z_elem_c(ji,jj,ibotm1) = 0._wp
345	z_elem_b(ji,jj,ibotm1) = 1._wp
346	en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
347	END DO
348	END DO
349	!
350	CASE ( 1 ) ! Neumman boundary condition
351	!
352	DO jj = tnldj, tnlej
353	DO ji = tnldi, tnlei
354	ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point
355	ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1
356	!
357	! Bottom level Dirichlet condition:
358	z_elem_a(ji,jj,ibot) = 0._wp
359	z_elem_c(ji,jj,ibot) = 0._wp
360	z_elem_b(ji,jj,ibot) = 1._wp
361	en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
362	!
363	! Just above last level: Neumann condition
364	z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b
365	z_elem_c(ji,jj,ibotm1) = 0._wp
366	END DO
367	END DO
368	!
369	END SELECT
370
371	! Matrix inversion (en prescribed at surface and the bottom)
372	! ----------------------------------------------------------
373	!
374	DO jk = 2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
375	DO jj = tnldj, tnlej
376	DO ji = tnldi, tnlei
377	z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)
378	END DO
379	END DO
380	END DO
381	DO jk = 2, jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
382	DO jj = tnldj, tnlej
383	DO ji = tnldi, tnlei
384	z_elem_a(ji,jj,jk) = en(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)
385	END DO
386	END DO
387	END DO
388	DO jk = jpk-1, 2, -1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
389	DO jj = tnldj, tnlej
390	DO ji = tnldi, tnlei
391	en(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * en(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)
392	END DO
393	END DO
394	END DO
395	! ! set the minimum value of tke
396	en(WRK_3D) = MAX( en(WRK_3D), rn_emin )
397
398	!!----------------------------------------!!
399	!! Solve prognostic equation for psi !!
400	!!----------------------------------------!!
401
402	! Set psi to previous time step value
403	!
404	SELECT CASE ( nn_clos )
405	!
406	CASE( 0 ) ! k-kl (Mellor-Yamada)
407	DO jk = 2, jpkm1
408	DO jj = tnldj, tnlej
409	DO ji = tnldi, tnlei
410	psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
411	END DO
412	END DO
413	END DO
414	!
415	CASE( 1 ) ! k-eps
416	DO jk = 2, jpkm1
417	DO jj = tnldj, tnlej
418	DO ji = tnldi, tnlei
419	psi(ji,jj,jk) = eps(ji,jj,jk)
420	END DO
421	END DO
422	END DO
423	!
424	CASE( 2 ) ! k-w
425	DO jk = 2, jpkm1
426	DO jj = tnldj, tnlej
427	DO ji = tnldi, tnlei
428	psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
429	END DO
430	END DO
431	END DO
432	!
433	CASE( 3 ) ! generic
434	DO jk = 2, jpkm1
435	DO jj = tnldj, tnlej
436	DO ji = tnldi, tnlei
437	psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn
438	END DO
439	END DO
440	END DO
441	!
442	END SELECT
443	!
444	! Now gls (output in psi)
445	! -------------------------------
446	! Resolution of a tridiagonal linear system by a "methode de chasse"
447	! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ).
448	! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal
449	! Warning : after this step, en : right hand side of the matrix
450
451	DO jk = 2, jpkm1
452	DO jj = tnldj, tnlej
453	DO ji = tnldi, tnlei
454	!
455	! psi / k
456	zratio = psi(ji,jj,jk) / eb(ji,jj,jk)
457	!
458	! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable)
459	dir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )
460	!
461	rpsi3 = dir * rpsi3m + ( 1._wp - dir ) * rpsi3p
462	!
463	! shear prod. - stratif. destruction
464	prod = rpsi1 * zratio * psh2(ji,jj,jk)
465	!
466	! stratif. destruction
467	buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) )
468	!
469	! shear prod. - stratif. destruction
470	diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
471	!
472	dir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
473	!
474	zesh2 = dir * ( prod + buoy ) + (1._wp - dir ) * prod ! production term
475	zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
476	!
477	! building the matrix
478	zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
479	! ! lower diagonal
480	z_elem_a(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
481	! ! upper diagonal
482	z_elem_c(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) )
483	! ! diagonal
484	z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk)
485	! ! right hand side in psi
486	psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk)
487	END DO
488	END DO
489	END DO
490	!
491	z_elem_b(WRK_2D,jpk) = 1._wp
492
493	! Surface boundary condition on psi
494	! ---------------------------------
495	!
496	SELECT CASE ( nn_bc_surf )
497	!
498	CASE ( 0 ) ! Dirichlet boundary conditions
499	!
500	! Surface value
501	zdep(WRK_2D) = zhsro(WRK_2D) * rl_sf ! Cosmetic
502	psi (WRK_2D,1) = rc0*rpp en(WRK_2D,1)*rmm zdep(WRK_2D)*rnn tmask(WRK_2D,1)
503	z_elem_a(WRK_2D,1) = psi(WRK_2D,1)
504	z_elem_c(WRK_2D,1) = 0._wp
505	z_elem_b(WRK_2D,1) = 1._wp
506	!
507	! One level below
508	zkar(WRK_2D) = (rl_sf + (vkarmn-rl_sf)(1._wp-exp(-rtransgdepw_n(WRK_2D,2)/zhsro(WRK_2D) )))
509	zdep(WRK_2D) = (zhsro(WRK_2D) + gdepw_n(WRK_2D,2)) * zkar(WRK_2D)
510	psi (WRK_2D,2) = rc0*rpp en(WRK_2D,2)*rmm zdep(WRK_2D)*rnn tmask(WRK_2D,1)
511	z_elem_a(WRK_2D,2) = 0._wp
512	z_elem_c(WRK_2D,2) = 0._wp
513	z_elem_b(WRK_2D,2) = 1._wp
514	!
515	!
516	CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz
517	!
518	! Surface value: Dirichlet
519	zdep(WRK_2D) = zhsro(WRK_2D) * rl_sf
520	psi (WRK_2D,1) = rc0*rpp en(WRK_2D,1)*rmm zdep(WRK_2D)*rnn tmask(WRK_2D,1)
521	z_elem_a(WRK_2D,1) = psi(WRK_2D,1)
522	z_elem_c(WRK_2D,1) = 0._wp
523	z_elem_b(WRK_2D,1) = 1._wp
524	!
525	! Neumann condition at k=2
526	z_elem_b(WRK_2D,2) = z_elem_b(WRK_2D,2) + z_elem_a(WRK_2D,2) ! Remove z_elem_a from z_elem_b
527	z_elem_a(WRK_2D,2) = 0._wp
528	!
529	! Set psi vertical flux at the surface:
530	zkar(WRK_2D) = rl_sf + (vkarmn-rl_sf)(1._wp-exp(-rtransgdept_n(WRK_2D,1)/zhsro(WRK_2D) )) ! Lengh scale slope
531	zdep(WRK_2D) = ((zhsro(WRK_2D) + gdept_n(WRK_2D,1)) / zhsro(WRK_2D))*(rmmra_sf)
532	zflxs(WRK_2D) = (rnn + rsbc_tke1 * (rnn + rmmra_sf) zdep(WRK_2D))(1._wp + rsbc_tke1zdep(WRK_2D))*(2._wprmm/3._wp-1_wp)
533	zdep(WRK_2D) = rsbc_psi1 * (zwall_psi(WRK_2D,1)p_avm(WRK_2D,1)+zwall_psi(WRK_2D,2)p_avm(WRK_2D,2)) * &
534	& ustars2(WRK_2D)*rmm zkar(WRK_2D)*rnn (zhsro(WRK_2D) + gdept_n(WRK_2D,1))**(rnn-1.)
535	zflxs(WRK_2D) = zdep(WRK_2D) * zflxs(WRK_2D)
536	psi(WRK_2D,2) = psi(WRK_2D,2) + zflxs(WRK_2D) / e3w_n(WRK_2D,2)
537	!
538	!
539	END SELECT
540
541	! Bottom boundary condition on psi
542	! --------------------------------
543	!
544	SELECT CASE ( nn_bc_bot )
545	!
546	!
547	CASE ( 0 ) ! Dirichlet
548	! ! en(ibot) = u^2 / Co2 and hmxl_n(ibot) = vkarmn rn_bfrz0
549	! ! Balance between the production and the dissipation terms
550	DO jj = tnldj, tnlej
551	DO ji = tnldi, tnlei
552	ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point
553	ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1
554	zdep(ji,jj) = vkarmn * rn_bfrz0
555	psi (ji,jj,ibot) = rc0*rpp en(ji,jj,ibot)*rmm zdep(ji,jj)**rnn
556	z_elem_a(ji,jj,ibot) = 0._wp
557	z_elem_c(ji,jj,ibot) = 0._wp
558	z_elem_b(ji,jj,ibot) = 1._wp
559	!
560	! Just above last level, Dirichlet condition again (GOTM like)
561	zdep(ji,jj) = vkarmn * ( rn_bfrz0 + e3t_n(ji,jj,ibotm1) )
562	psi (ji,jj,ibotm1) = rc0*rpp en(ji,jj,ibot )*rmm zdep(ji,jj)**rnn
563	z_elem_a(ji,jj,ibotm1) = 0._wp
564	z_elem_c(ji,jj,ibotm1) = 0._wp
565	z_elem_b(ji,jj,ibotm1) = 1._wp
566	END DO
567	END DO
568	!
569	CASE ( 1 ) ! Neumman boundary condition
570	!
571	DO jj = tnldj, tnlej
572	DO ji = tnldi, tnlei
573	ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point
574	ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1
575	!
576	! Bottom level Dirichlet condition:
577	zdep(ji,jj) = vkarmn * rn_bfrz0
578	psi (ji,jj,ibot) = rc0*rpp en(ji,jj,ibot)*rmm zdep(ji,jj)**rnn
579	!
580	z_elem_a(ji,jj,ibot) = 0._wp
581	z_elem_c(ji,jj,ibot) = 0._wp
582	z_elem_b(ji,jj,ibot) = 1._wp
583	!
584	! Just above last level: Neumann condition with flux injection
585	z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b
586	z_elem_c(ji,jj,ibotm1) = 0.
587	!
588	! Set psi vertical flux at the bottom:
589	zdep(ji,jj) = rn_bfrz0 + 0.5_wp*e3t_n(ji,jj,ibotm1)
590	zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) ) &
591	& * (0.5_wp(en(ji,jj,ibot)+en(ji,jj,ibotm1)))rmm zdep(ji,jj)**(rnn-1._wp)
592	psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1)
593	END DO
594	END DO
595	!
596	END SELECT
597
598	! Matrix inversion
599	! ----------------
600	!
601	DO jk = 2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
602	DO jj = tnldj, tnlej
603	DO ji = tnldi, tnlei
604	z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)
605	END DO
606	END DO
607	END DO
608	DO jk = 2, jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
609	DO jj = tnldj, tnlej
610	DO ji = tnldi, tnlei
611	z_elem_a(ji,jj,jk) = psi(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)
612	END DO
613	END DO
614	END DO
615	DO jk = jpk-1, 2, -1 ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
616	DO jj = tnldj, tnlej
617	DO ji = tnldi, tnlei
618	psi(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * psi(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)
619	END DO
620	END DO
621	END DO
622
623	! Set dissipation
624	!----------------
625
626	SELECT CASE ( nn_clos )
627	!
628	CASE( 0 ) ! k-kl (Mellor-Yamada)
629	DO jk = 1, jpkm1
630	DO jj = tnldj, tnlej
631	DO ji = tnldi, tnlei
632	eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin)
633	END DO
634	END DO
635	END DO
636	!
637	CASE( 1 ) ! k-eps
638	DO jk = 1, jpkm1
639	DO jj = tnldj, tnlej
640	DO ji = tnldi, tnlei
641	eps(ji,jj,jk) = psi(ji,jj,jk)
642	END DO
643	END DO
644	END DO
645	!
646	CASE( 2 ) ! k-w
647	DO jk = 1, jpkm1
648	DO jj = tnldj, tnlej
649	DO ji = tnldi, tnlei
650	eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)
651	END DO
652	END DO
653	END DO
654	!
655	CASE( 3 ) ! generic
656	zcoef = rc0**( 3._wp + rpp/rnn )
657	zex1 = ( 1.5_wp + rmm/rnn )
658	zex2 = -1._wp / rnn
659	DO jk = 1, jpkm1
660	DO jj = tnldj, tnlej
661	DO ji = tnldi, tnlei
662	eps(ji,jj,jk) = zcoef * en(ji,jj,jk)*zex1 psi(ji,jj,jk)**zex2
663	END DO
664	END DO
665	END DO
666	!
667	END SELECT
668
669	! Limit dissipation rate under stable stratification
670	! --------------------------------------------------
671	DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time
672	DO jj = tnldj, tnlej
673	DO ji = tnldi, tnlei
674	! limitation
675	eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin )
676	hmxl_n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
677	! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
678	zrn2 = MAX( rn2(ji,jj,jk), rsmall )
679	IF( ln_length_lim ) hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) )
680	END DO
681	END DO
682	END DO
683
684	!
685	! Stability function and vertical viscosity and diffusivity
686	! ---------------------------------------------------------
687	!
688	SELECT CASE ( nn_stab_func )
689	!
690	CASE ( 0 , 1 ) ! Galperin or Kantha-Clayson stability functions
691	DO jk = 2, jpkm1
692	DO jj = tnldj, tnlej
693	DO ji = tnldi, tnlei
694	! zcof = l²/q²
695	zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
696	! Gh = -N²l²/q²
697	gh = - rn2(ji,jj,jk) * zcof
698	gh = MIN( gh, rgh0 )
699	gh = MAX( gh, rghmin )
700	! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
701	sh = ra2( 1._wp-6._wpra1/rb1 ) / ( 1.-3.ra2gh(6.ra1+rb2*( 1._wp-rc3 ) ) )
702	sm = ( rb1*(-1._wp/3._wp) + ( 18._wpra1ra1 + 9._wpra1ra2(1._wp-rc2) )shgh ) / (1._wp-9._wpra1ra2*gh)
703	!
704	! Store stability function in zstt and zstm
705	zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
706	zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
707	END DO
708	END DO
709	END DO
710	!
711	CASE ( 2, 3 ) ! Canuto stability functions
712	DO jk = 2, jpkm1
713	DO jj = tnldj, tnlej
714	DO ji = tnldi, tnlei
715	! zcof = l²/q²
716	zcof = hmxl_b(ji,jj,jk)hmxl_b(ji,jj,jk) / ( 2._wp eb(ji,jj,jk) )
717	! Gh = -N²l²/q²
718	gh = - rn2(ji,jj,jk) * zcof
719	gh = MIN( gh, rgh0 )
720	gh = MAX( gh, rghmin )
721	gh = gh * rf6
722	! Gm = M²l²/q² Shear number
723	shr = psh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall )
724	gm = MAX( shr * zcof , 1.e-10 )
725	gm = gm * rf6
726	gm = MIN ( (rd0 - rd1gh + rd3ghgh) / (rd2-rd4gh) , gm )
727	! Stability functions from Canuto
728	rcff = rd0 - rd1gh +rd2gm + rd3ghgh - rd4ghgm + rd5gmgm
729	sm = (rs0 - rs1gh + rs2gm) / rcff
730	sh = (rs4 - rs5gh + rs6gm) / rcff
731	!
732	! Store stability function in zstt and zstm
733	zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
734	zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
735	END DO
736	END DO
737	END DO
738	!
739	END SELECT
740
741	! Boundary conditions on stability functions for momentum (Neumann):
742	! Lines below are useless if GOTM style Dirichlet conditions are used
743
744	zstm(WRK_2D,1) = zstm(WRK_2D,2)
745
746	!!gm should be done for ISF (top boundary cond.)
747	DO jj = tnldj, tnlej
748	DO ji = tnldi, tnlei
749	zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
750	END DO
751	END DO
752
753	! Compute diffusivities/viscosities
754	! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
755	DO jk = 1, jpk
756	DO jj = tnldj, tnlej
757	DO ji = tnldi, tnlei
758	zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
759	zav = zsqen * zstt(ji,jj,jk)
760	p_avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
761	zav = zsqen * zstm(ji,jj,jk)
762	p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom
763	END DO
764	END DO
765	END DO
766	p_avt(WRK_2D,1) = 0._wp
767	!
768	IF(ln_ctl) THEN
769	CALL prt_ctl( tab3d_1=en , clinfo1=' gls - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
770	CALL prt_ctl( tab3d_1=avm, clinfo1=' gls - m: ', ovlap=1, kdim=jpk )
771	ENDIF
772	!
773	IF( nn_timing == 1 ) CALL timing_stop('zdf_gls')
774	!
775	END SUBROUTINE zdf_gls
776
777
778	SUBROUTINE zdf_gls_init
779	!!----------------------------------------------------------------------
780	!! * ROUTINE zdf_gls_init *
781	!!
782	!! ** Purpose : Initialization of the vertical eddy diffivity and
783	!! viscosity computed using a GLS turbulent closure scheme
784	!!
785	!! ** Method : Read the namzdf_gls namelist and check the parameters
786	!!
787	!! ** input : Namlist namzdf_gls
788	!!
789	!! ** Action : Increase by 1 the nstop flag is setting problem encounter
790	!!
791	!!----------------------------------------------------------------------
792	USE dynzdf_exp
793	USE trazdf_exp
794	!
795	INTEGER :: jk ! dummy loop indices
796	INTEGER :: ios ! Local integer output status for namelist read
797	REAL(wp):: zcr ! local scalar
798	!!
799	NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
800	& rn_clim_galp, ln_sigpsi, rn_hsro, &
801	& rn_crban, rn_charn, rn_frac_hs, &
802	& nn_bc_surf, nn_bc_bot, nn_z0_met, &
803	& nn_stab_func, nn_clos
804	!!----------------------------------------------------------
805	!
806	IF( nn_timing == 1 ) CALL timing_start('zdf_gls_init')
807	!
808	REWIND( numnam_ref ) ! Namelist namzdf_gls in reference namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme
809	READ ( numnam_ref, namzdf_gls, IOSTAT = ios, ERR = 901)
810	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in reference namelist', lwp )
811
812	REWIND( numnam_cfg ) ! Namelist namzdf_gls in configuration namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme
813	READ ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 )
814	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist', lwp )
815	IF(lwm) WRITE ( numond, namzdf_gls )
816
817	IF(lwp) THEN !* Control print
818	WRITE(numout,*)
819	WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme'
820	WRITE(numout,*) '~~~~~~~~~~~~'
821	WRITE(numout,*) ' Namelist namzdf_gls : set gls mixing parameters'
822	WRITE(numout,*) ' minimum value of en rn_emin = ', rn_emin
823	WRITE(numout,*) ' minimum value of eps rn_epsmin = ', rn_epsmin
824	WRITE(numout,*) ' Limit dissipation rate under stable stratif. ln_length_lim = ', ln_length_lim
825	WRITE(numout,*) ' Galperin limit (Standard: 0.53, Holt: 0.26) rn_clim_galp = ', rn_clim_galp
826	WRITE(numout,*) ' TKE Surface boundary condition nn_bc_surf = ', nn_bc_surf
827	WRITE(numout,*) ' TKE Bottom boundary condition nn_bc_bot = ', nn_bc_bot
828	WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi
829	WRITE(numout,*) ' Craig and Banner coefficient rn_crban = ', rn_crban
830	WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn
831	WRITE(numout,*) ' Surface roughness formula nn_z0_met = ', nn_z0_met
832	WRITE(numout,*) ' Wave height frac. (used if nn_z0_met=2) rn_frac_hs = ', rn_frac_hs
833	WRITE(numout,*) ' Stability functions nn_stab_func = ', nn_stab_func
834	WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos
835	WRITE(numout,*) ' Surface roughness (m) rn_hsro = ', rn_hsro
836	WRITE(numout,*) ' Bottom roughness (m) (nambfr namelist) rn_bfrz0 = ', rn_bfrz0
837	WRITE(numout,*)
838	ENDIF
839
840	! !* allocate GLS arrays
841	IF( zdf_gls_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
842
843	! !* Check of some namelist values
844	IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )
845	IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )
846	IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' )
847	IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )
848	IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' )
849	IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
850
851	SELECT CASE ( nn_clos ) !* set the parameters for the chosen closure
852	!
853	CASE( 0 ) ! k-kl (Mellor-Yamada)
854	!
855	IF(lwp) WRITE(numout,*) ' ==>> k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)'
856	IF(lwp) WRITE(numout,*)
857	rpp = 0._wp
858	rmm = 1._wp
859	rnn = 1._wp
860	rsc_tke = 1.96_wp
861	rsc_psi = 1.96_wp
862	rpsi1 = 0.9_wp
863	rpsi3p = 1._wp
864	rpsi2 = 0.5_wp
865	!
866	SELECT CASE ( nn_stab_func )
867	CASE( 0, 1 ) ; rpsi3m = 2.53_wp ! G88 or KC stability functions
868	CASE( 2 ) ; rpsi3m = 2.62_wp ! Canuto A stability functions
869	CASE( 3 ) ; rpsi3m = 2.38 ! Canuto B stability functions (caution : constant not identified)
870	END SELECT
871	!
872	CASE( 1 ) ! k-eps
873	!
874	IF(lwp) WRITE(numout,*) ' ==>> k-eps closure chosen'
875	IF(lwp) WRITE(numout,*)
876	rpp = 3._wp
877	rmm = 1.5_wp
878	rnn = -1._wp
879	rsc_tke = 1._wp
880	rsc_psi = 1.2_wp ! Schmidt number for psi
881	rpsi1 = 1.44_wp
882	rpsi3p = 1._wp
883	rpsi2 = 1.92_wp
884	!
885	SELECT CASE ( nn_stab_func )
886	CASE( 0, 1 ) ; rpsi3m = -0.52_wp ! G88 or KC stability functions
887	CASE( 2 ) ; rpsi3m = -0.629_wp ! Canuto A stability functions
888	CASE( 3 ) ; rpsi3m = -0.566 ! Canuto B stability functions
889	END SELECT
890	!
891	CASE( 2 ) ! k-omega
892	!
893	IF(lwp) WRITE(numout,*) ' ==>> k-omega closure chosen'
894	IF(lwp) WRITE(numout,*)
895	rpp = -1._wp
896	rmm = 0.5_wp
897	rnn = -1._wp
898	rsc_tke = 2._wp
899	rsc_psi = 2._wp
900	rpsi1 = 0.555_wp
901	rpsi3p = 1._wp
902	rpsi2 = 0.833_wp
903	!
904	SELECT CASE ( nn_stab_func )
905	CASE( 0, 1 ) ; rpsi3m = -0.58_wp ! G88 or KC stability functions
906	CASE( 2 ) ; rpsi3m = -0.64_wp ! Canuto A stability functions
907	CASE( 3 ) ; rpsi3m = -0.64_wp ! Canuto B stability functions caution : constant not identified)
908	END SELECT
909	!
910	CASE( 3 ) ! generic
911	!
912	IF(lwp) WRITE(numout,*) ' ==>> generic closure chosen'
913	IF(lwp) WRITE(numout,*)
914	rpp = 2._wp
915	rmm = 1._wp
916	rnn = -0.67_wp
917	rsc_tke = 0.8_wp
918	rsc_psi = 1.07_wp
919	rpsi1 = 1._wp
920	rpsi3p = 1._wp
921	rpsi2 = 1.22_wp
922	!
923	SELECT CASE ( nn_stab_func )
924	CASE( 0, 1 ) ; rpsi3m = 0.1_wp ! G88 or KC stability functions
925	CASE( 2 ) ; rpsi3m = 0.05_wp ! Canuto A stability functions
926	CASE( 3 ) ; rpsi3m = 0.05_wp ! Canuto B stability functions caution : constant not identified)
927	END SELECT
928	!
929	END SELECT
930
931	!
932	SELECT CASE ( nn_stab_func ) !* set the parameters of the stability functions
933	!
934	CASE ( 0 ) ! Galperin stability functions
935	!
936	IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Galperin'
937	rc2 = 0._wp
938	rc3 = 0._wp
939	rc_diff = 1._wp
940	rc0 = 0.5544_wp
941	rcm_sf = 0.9884_wp
942	rghmin = -0.28_wp
943	rgh0 = 0.0233_wp
944	rghcri = 0.02_wp
945	!
946	CASE ( 1 ) ! Kantha-Clayson stability functions
947	!
948	IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Kantha-Clayson'
949	rc2 = 0.7_wp
950	rc3 = 0.2_wp
951	rc_diff = 1._wp
952	rc0 = 0.5544_wp
953	rcm_sf = 0.9884_wp
954	rghmin = -0.28_wp
955	rgh0 = 0.0233_wp
956	rghcri = 0.02_wp
957	!
958	CASE ( 2 ) ! Canuto A stability functions
959	!
960	IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Canuto A'
961	rs0 = 1.5_wp * rl1 * rl5*rl5
962	rs1 = -rl4(rl6+rl7) + 2._wprl4rl5(rl1-(1._wp/3._wp)rl2-rl3) + 1.5_wprl1rl5rl8
963	rs2 = -(3._wp/8._wp) * rl1(rl6rl6-rl7*rl7)
964	rs4 = 2._wp * rl5
965	rs5 = 2._wp * rl4
966	rs6 = (2._wp/3._wp) * rl5 * ( 3._wprl3rl3 - rl2rl2 ) - 0.5_wp rl5rl1 (3._wp*rl3-rl2) &
967	& + 0.75_wp * rl1 * ( rl6 - rl7 )
968	rd0 = 3._wp * rl5*rl5
969	rd1 = rl5 * ( 7._wprl4 + 3._wprl8 )
970	rd2 = rl5rl5 ( 3._wprl3rl3 - rl2rl2 ) - 0.75_wp(rl6rl6 - rl7rl7 )
971	rd3 = rl4 * ( 4._wprl4 + 3._wprl8)
972	rd4 = rl4 * ( rl2 * rl6 - 3._wprl3rl7 - rl5(rl2rl2 - rl3rl3 ) ) + rl5rl8 * ( 3._wprl3rl3 - rl2*rl2 )
973	rd5 = 0.25_wp * ( rl2rl2 - 3._wp rl3rl3 ) ( rl6rl6 - rl7rl7 )
974	rc0 = 0.5268_wp
975	rf6 = 8._wp / (rc0**6._wp)
976	rc_diff = SQRT(2._wp) / (rc0**3._wp)
977	rcm_sf = 0.7310_wp
978	rghmin = -0.28_wp
979	rgh0 = 0.0329_wp
980	rghcri = 0.03_wp
981	!
982	CASE ( 3 ) ! Canuto B stability functions
983	!
984	IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Canuto B'
985	rs0 = 1.5_wp * rm1 * rm5*rm5
986	rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4rm5(rm1-(1._wp/3._wp)rm2-rm3) + 1.5_wp rm1rm5rm8
987	rs2 = -(3._wp/8._wp) * rm1 * (rm6rm6-rm7rm7 )
988	rs4 = 2._wp * rm5
989	rs5 = 2._wp * rm4
990	rs6 = (2._wp/3._wp) * rm5 * (3._wprm3rm3-rm2rm2) - 0.5_wp rm5rm1(3._wprm3-rm2) + 0.75_wp rm1*(rm6-rm7)
991	rd0 = 3._wp * rm5*rm5
992	rd1 = rm5 * (7._wprm4 + 3._wprm8)
993	rd2 = rm5rm5 (3._wprm3rm3 - rm2rm2) - 0.75_wp (rm6rm6 - rm7rm7)
994	rd3 = rm4 * ( 4._wprm4 + 3._wprm8 )
995	rd4 = rm4 * ( rm2rm6 -3._wprm3rm7 - rm5(rm2rm2 - rm3rm3) ) + rm5 * rm8 * ( 3._wprm3rm3 - rm2*rm2 )
996	rd5 = 0.25_wp * ( rm2rm2 - 3._wprm3rm3 ) ( rm6rm6 - rm7rm7 )
997	rc0 = 0.5268_wp !! rc0 = 0.5540_wp (Warner ...) to verify !
998	rf6 = 8._wp / ( rc0**6._wp )
999	rc_diff = SQRT(2._wp)/(rc0**3.)
1000	rcm_sf = 0.7470_wp
1001	rghmin = -0.28_wp
1002	rgh0 = 0.0444_wp
1003	rghcri = 0.0414_wp
1004	!
1005	END SELECT
1006
1007	! !* Set Schmidt number for psi diffusion in the wave breaking case
1008	! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
1009	! ! or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
1010	IF( ln_sigpsi ) THEN
1011	ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf
1012	! Verification: retrieve Burchard (2001) results by uncomenting the line below:
1013	! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work
1014	! ra_sf = -SQRT(2./3.rc03./rn_cm_sfrn_sc_tke)/vkarmn
1015	rsc_psi0 = rsc_tke/(24.rpsi2)(-1.+(4.rnn + ra_sf(1.+4.rmm))2./(ra_sf*2.))
1016	ELSE
1017	rsc_psi0 = rsc_psi
1018	ENDIF
1019
1020	! !* Shear free turbulence parameters
1021	!
1022	ra_sf = -4._wprnnSQRT(rsc_tke) / ( (1._wp+4._wprmm)SQRT(rsc_tke) &
1023	& - SQRT(rsc_tke + 24._wprsc_psi0rpsi2 ) )
1024
1025	IF ( rn_crban==0._wp ) THEN
1026	rl_sf = vkarmn
1027	ELSE
1028	rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wprmm + 8._wprmm*2_wp)rsc_tke &
1029	& + 12._wp * rsc_psi0rpsi2 - (1._wp + 4._wprmm) &
1030	& SQRT(rsc_tke(rsc_tke &
1031	& + 24._wprsc_psi0rpsi2)) ) &
1032	& /(12._wprnn*2.) &
1033	& )
1034	ENDIF
1035
1036	!
1037	IF(lwp) THEN !* Control print
1038	WRITE(numout,*)
1039	WRITE(numout,*) ' Limit values :'
1040	WRITE(numout,*) ' Parameter m = ',rmm
1041	WRITE(numout,*) ' Parameter n = ',rnn
1042	WRITE(numout,*) ' Parameter p = ',rpp
1043	WRITE(numout,*) ' rpsi1 = ',rpsi1
1044	WRITE(numout,*) ' rpsi2 = ',rpsi2
1045	WRITE(numout,*) ' rpsi3m = ',rpsi3m
1046	WRITE(numout,*) ' rpsi3p = ',rpsi3p
1047	WRITE(numout,*) ' rsc_tke = ',rsc_tke
1048	WRITE(numout,*) ' rsc_psi = ',rsc_psi
1049	WRITE(numout,*) ' rsc_psi0 = ',rsc_psi0
1050	WRITE(numout,*) ' rc0 = ',rc0
1051	WRITE(numout,*)
1052	WRITE(numout,*) ' Shear free turbulence parameters:'
1053	WRITE(numout,*) ' rcm_sf = ',rcm_sf
1054	WRITE(numout,*) ' ra_sf = ',ra_sf
1055	WRITE(numout,*) ' rl_sf = ',rl_sf
1056	ENDIF
1057
1058	! !* Constants initialization
1059	rc02 = rc0 * rc0 ; rc02r = 1. / rc02
1060	rc03 = rc02 * rc0
1061	rc04 = rc03 * rc0
1062	rsbc_tke1 = -3._wp/2._wprn_crbanra_sf*rl_sf ! Dirichlet + Wave breaking
1063	rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking
1064	zcr = MAX(rsmall, rsbc_tke1*(1./(-ra_sf3._wp/2._wp))-1._wp )
1065	rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer
1066	rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness
1067	rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness
1068	rsbc_psi1 = -0.5_wp * rdt * rc0*(rpp-2._wprmm) / rsc_psi
1069	rsbc_psi2 = -0.5_wp * rdt * rc0*rpp rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking
1070	!
1071	rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke
1072	rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke
1073	!
1074	! !* Wall proximity function
1075	zwall(:,:,:) = 1._wp * tmask(:,:,:)
1076
1077	! !* read or initialize all required files
1078	CALL gls_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, hmxl_n)
1079	!
1080	IF( nn_timing == 1 ) CALL timing_stop('zdf_gls_init')
1081	!
1082	END SUBROUTINE zdf_gls_init
1083
1084
1085	SUBROUTINE gls_rst( kt, cdrw )
1086	!!---------------------------------------------------------------------
1087	!! * ROUTINE ts_rst *
1088	!!
1089	!! ** Purpose : Read or write TKE file (en) in restart file
1090	!!
1091	!! ** Method : use of IOM library
1092	!! if the restart does not contain TKE, en is either
1093	!! set to rn_emin or recomputed (nn_igls/=0)
1094	!!----------------------------------------------------------------------
1095	INTEGER , INTENT(in) :: kt ! ocean time-step
1096	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1097	!
1098	INTEGER :: jit, jk ! dummy loop indices
1099	INTEGER :: id1, id2, id3, id4
1100	INTEGER :: ji, jj, ikbu, ikbv
1101	REAL(wp):: cbx, cby
1102	!!----------------------------------------------------------------------
1103	!
1104	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
1105	! ! ---------------
1106	IF( ln_rstart ) THEN !* Read the restart file
1107	id1 = iom_varid( numror, 'en' , ldstop = .FALSE. )
1108	id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. )
1109	id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. )
1110	id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. )
1111	!
1112	IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! all required arrays exist
1113	CALL iom_get( numror, jpdom_autoglo, 'en' , en )
1114	CALL iom_get( numror, jpdom_autoglo, 'avt_k' , avt_k )
1115	CALL iom_get( numror, jpdom_autoglo, 'avm_k' , avm_k )
1116	CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n )
1117	ELSE
1118	IF(lwp) WRITE(numout,*)
1119	IF(lwp) WRITE(numout,*) ' ==>> previous run without GLS scheme, set en and hmxl_n to background values'
1120	en (:,:,:) = rn_emin
1121	hmxl_n(:,:,:) = 0.05_wp
1122	! avt_k, avm_k already set to the background value in zdf_phy_init
1123	ENDIF
1124	ELSE !* Start from rest
1125	IF(lwp) WRITE(numout,*)
1126	IF(lwp) WRITE(numout,*) ' ==>> start from rest, set en and hmxl_n by background values'
1127	en (:,:,:) = rn_emin
1128	hmxl_n(:,:,:) = 0.05_wp
1129	! avt_k, avm_k already set to the background value in zdf_phy_init
1130	ENDIF
1131	!
1132	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
1133	! ! -------------------
1134	IF(lwp) WRITE(numout,*) '---- gls-rst ----'
1135	CALL iom_rstput( kt, nitrst, numrow, 'en' , en )
1136	CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k )
1137	CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k )
1138	CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n )
1139	!
1140	ENDIF
1141	!
1142	END SUBROUTINE gls_rst
1143
1144	!!======================================================================
1145	END MODULE zdfgls
1146

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 8279

Download in other formats: