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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

zdfgls.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 4443

Last change on this file since 4443 was 4443, checked in by trackstand2, 10 years ago
Use of mbkmax in zdfgls
Property svn:keywords set to `Id`
File size: 66.9 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	!!----------------------------------------------------------------------
10	#if defined key_zdfgls \|\| defined key_esopa
11	!!----------------------------------------------------------------------
12	!! 'key_zdfgls' Generic Length Scale vertical physics
13	!!----------------------------------------------------------------------
14	!! zdf_gls : update momentum and tracer Kz from a gls scheme
15	!! zdf_gls_init : initialization, namelist read, and parameters control
16	!! gls_rst : read/write gls restart in ocean restart file
17	!!----------------------------------------------------------------------
18	USE oce ! ocean dynamics and active tracers
19	USE dom_oce ! ocean space and time domain
20	USE domvvl ! ocean space and time domain : variable volume layer
21	USE zdf_oce ! ocean vertical physics
22	USE sbc_oce ! surface boundary condition: ocean
23	USE phycst ! physical constants
24	USE zdfmxl ! mixed layer
25	USE restart ! only for lrst_oce
26	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
27	USE lib_mpp ! MPP manager
28	USE prtctl ! Print control
29	USE in_out_manager ! I/O manager
30	USE iom ! I/O manager library
31
32	IMPLICIT NONE
33	PRIVATE
34
35	PUBLIC zdf_gls ! routine called in step module
36	PUBLIC zdf_gls_init ! routine called in opa module
37	PUBLIC gls_rst ! routine called in step module
38
39	LOGICAL , PUBLIC, PARAMETER :: lk_zdfgls = .TRUE. !: TKE vertical mixing flag
40	!
41	!! DCSE_NEMO: does not need to be public
42	! REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy
43	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy
44
45	!! DCSE_NEMO: does not need to be public
46	! REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mxln !: now mixing length
47	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mxln !: now mixing length
48
49	!! DCSE_NEMO: does not need to be public
50	! REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zwall !: wall function
51	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zwall !: wall function
52
53	!! DCSE_NEMO: does not need to be public
54	! REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points
55	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points
56
57	!! DCSE_NEMO: does not need to be public
58	! REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustarb2 !: Squared bottom velocity scale at T-points
59	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustarb2 !: Squared bottom velocity scale at T-points
60
61	! !!! Namelist namzdf_gls
62	LOGICAL :: ln_crban = .FALSE. ! =T use Craig and Banner scheme
63	LOGICAL :: ln_length_lim = .FALSE. ! use limit on the dissipation rate under
64	! stable stratification (Galperin et al. 1988)
65	LOGICAL :: ln_sigpsi = .FALSE. ! Activate Burchard (2003) modification for
66	! k-eps closure & wave breaking mixing
67	INTEGER :: nn_tkebc_surf = 0 ! TKE surface boundary condition (=0/1)
68	INTEGER :: nn_tkebc_bot = 0 ! TKE bottom boundary condition (=0/1)
69	INTEGER :: nn_psibc_surf = 0 ! PSI surface boundary condition (=0/1)
70	INTEGER :: nn_psibc_bot = 0 ! PSI bottom boundary condition (=0/1)
71	INTEGER :: nn_stab_func = 0 ! stability functions G88, KC or Canuto (=0/1/2)
72	INTEGER :: nn_clos = 0 ! closure 0/1/2/3 MY82/k-eps/k-w/gen
73	REAL(wp) :: rn_clim_galp = 0.53_wp ! Holt 2008 value for k-eps: 0.267
74	REAL(wp) :: rn_epsmin = 1.e-12_wp ! minimum value of dissipation (m2/s3)
75	REAL(wp) :: rn_emin = 1.e-6_wp ! minimum value of TKE (m2/s2)
76	REAL(wp) :: rn_charn = 2.e+5_wp ! Charnock constant for surface breaking waves
77	! mixing : 1400. (standard) or 2.e5 (Stacey value)
78	REAL(wp) :: rn_crban = 100._wp ! Craig and Banner constant for surface breaking waves mixing
79
80	REAL(wp) :: hsro = 0.003_wp ! Minimum surface roughness
81	REAL(wp) :: hbro = 0.003_wp ! Bottom roughness (m)
82	REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters
83	REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1
84	REAL(wp) :: rl_sf = 0.2_wp ! 0 <rl_sf<vkarmn
85	REAL(wp) :: rghmin = -0.28_wp
86	REAL(wp) :: rgh0 = 0.0329_wp
87	REAL(wp) :: rghcri = 0.03_wp
88	REAL(wp) :: ra1 = 0.92_wp
89	REAL(wp) :: ra2 = 0.74_wp
90	REAL(wp) :: rb1 = 16.60_wp
91	REAL(wp) :: rb2 = 10.10_wp
92	REAL(wp) :: re2 = 1.33_wp
93	REAL(wp) :: rl1 = 0.107_wp
94	REAL(wp) :: rl2 = 0.0032_wp
95	REAL(wp) :: rl3 = 0.0864_wp
96	REAL(wp) :: rl4 = 0.12_wp
97	REAL(wp) :: rl5 = 11.9_wp
98	REAL(wp) :: rl6 = 0.4_wp
99	REAL(wp) :: rl7 = 0.0_wp
100	REAL(wp) :: rl8 = 0.48_wp
101	REAL(wp) :: rm1 = 0.127_wp
102	REAL(wp) :: rm2 = 0.00336_wp
103	REAL(wp) :: rm3 = 0.0906_wp
104	REAL(wp) :: rm4 = 0.101_wp
105	REAL(wp) :: rm5 = 11.2_wp
106	REAL(wp) :: rm6 = 0.4_wp
107	REAL(wp) :: rm7 = 0.0_wp
108	REAL(wp) :: rm8 = 0.318_wp
109
110	REAL(wp) :: rc02, rc02r, rc03, rc04 ! coefficients deduced from above parameters
111	REAL(wp) :: rc03_sqrt2_galp ! - - - -
112	REAL(wp) :: rsbc_tke1, rsbc_tke2, rsbc_tke3, rfact_tke ! - - - -
113	REAL(wp) :: rsbc_psi1, rsbc_psi2, rsbc_psi3, rfact_psi ! - - - -
114	REAL(wp) :: rsbc_mb , rsbc_std , rsbc_zs ! - - - -
115	REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff ! - - - -
116	REAL(wp) :: rs0, rs1, rs2, rs4, rs5, rs6 ! - - - -
117	REAL(wp) :: rd0, rd1, rd2, rd3, rd4, rd5 ! - - - -
118	REAL(wp) :: rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0 ! - - - -
119	REAL(wp) :: rpsi3m, rpsi3p, rpp, rmm, rnn ! - - - -
120
121	!! * Control permutation of array indices
122	# include "oce_ftrans.h90"
123	# include "dom_oce_ftrans.h90"
124	# include "domvvl_ftrans.h90"
125	# include "zdf_oce_ftrans.h90"
126	# include "sbc_oce_ftrans.h90"
127	!! DCSE_NEMO: private module variables do not need their own directives file
128	!FTRANS en mxln zwall :I :I :z
129
130	!! * Substitutions
131	# include "domzgr_substitute.h90"
132	# include "vectopt_loop_substitute.h90"
133	!!----------------------------------------------------------------------
134	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
135	!! $Id$
136	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
137	!!----------------------------------------------------------------------
138	CONTAINS
139
140	INTEGER FUNCTION zdf_gls_alloc()
141	!!----------------------------------------------------------------------
142	!! * FUNCTION zdf_gls_alloc *
143	!!----------------------------------------------------------------------
144	ALLOCATE( en(jpi,jpj,jpk), mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , &
145	& ustars2(jpi,jpj), ustarb2(jpi,jpj) , STAT= zdf_gls_alloc )
146	!
147	IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc )
148	IF( zdf_gls_alloc /= 0 ) CALL ctl_warn('zdf_gls_alloc: failed to allocate arrays')
149	END FUNCTION zdf_gls_alloc
150
151
152	SUBROUTINE zdf_gls( kt )
153	!!----------------------------------------------------------------------
154	!! * ROUTINE zdf_gls *
155	!!
156	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
157	!! coefficients using the GLS turbulent closure scheme.
158	!!----------------------------------------------------------------------
159	USE oce, z_elem_a => ua ! use ua as workspace
160	USE oce, z_elem_b => va ! use va as workspace
161	USE oce, z_elem_c => ta ! use ta as workspace
162	USE oce, psi => sa ! use sa as workspace
163	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
164	USE wrk_nemo, ONLY: zdep => wrk_2d_1
165	USE wrk_nemo, ONLY: zflxs => wrk_2d_2 ! Turbulence fluxed induced by internal waves
166	USE wrk_nemo, ONLY: zhsro => wrk_2d_3 ! Surface roughness (surface waves)
167	USE wrk_nemo, ONLY: eb => wrk_3d_1 ! tke at time before
168	USE wrk_nemo, ONLY: mxlb => wrk_3d_2 ! mixing length at time before
169	USE wrk_nemo, ONLY: shear => wrk_3d_3 ! vertical shear
170	USE wrk_nemo, ONLY: eps => wrk_3d_4 ! dissipation rate
171	USE wrk_nemo, ONLY: zwall_psi => wrk_3d_5 ! Wall function use in the wb case (ln_sigpsi.AND.ln_crban=T)
172	USE timing, ONLY: timing_start, timing_stop
173
174	!! DCSE_NEMO: need additional directives for renamed module variables
175	!FTRANS z_elem_a z_elem_b z_elem_c psi :I :I :z
176	!FTRANS eb mxlb shear eps zwall_psi :I :I :z
177	!
178	INTEGER, INTENT(in) :: kt ! ocean time step
179	INTEGER :: ji, jj, jk, ibot, ibotm1, dir ! dummy loop arguments
180	REAL(wp) :: zesh2, zsigpsi, zcoef, zex1, zex2 ! local scalars
181	REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - -
182	REAL(wp) :: zratio, zrn2, zflxb, sh ! - -
183	REAL(wp) :: prod, buoy, diss, zdiss, sm ! - -
184	REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - -
185	!!--------------------------------------------------------------------
186
187	CALL timing_start('zdf_gls')
188
189	IF( wrk_in_use(2, 1,2,3) .OR. wrk_in_use(3, 1,2,3,4,5) ) THEN
190	CALL ctl_stop('zdf_gls: requested workspace arrays unavailable.') ; RETURN
191	END IF
192
193	! Preliminary computing
194
195	ustars2 = 0._wp ; ustarb2 = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp
196
197	! Compute surface and bottom friction at T-points
198	!CDIR NOVERRCHK
199	DO jj = 2, jpjm1
200	!CDIR NOVERRCHK
201	DO ji = fs_2, fs_jpim1 ! vector opt.
202	!
203	! surface friction
204	#if defined key_z_first
205	ustars2(ji,jj) = rau0r * taum(ji,jj) * tmask_1(ji,jj)
206	#else
207	ustars2(ji,jj) = rau0r * taum(ji,jj) * tmask(ji,jj,1)
208	#endif
209	!
210	! bottom friction (explicit before friction)
211	! Note that we chose here not to bound the friction as in dynbfr)
212	#if defined key_z_first
213	ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) &
214	& * ( 1._wp - 0.5_wp * umask_1(ji,jj) * umask_1(ji-1,jj) )
215	zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) &
216	& * ( 1._wp - 0.5_wp * vmask_1(ji,jj) * vmask_1(ji,jj-1) )
217	ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask_1(ji,jj)
218	#else
219	ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) &
220	& * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) )
221	zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) &
222	& * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) )
223	ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
224	#endif
225	END DO
226	END DO
227
228	! In case of breaking surface waves mixing,
229	! Compute surface roughness length according to Charnock formula:
230	IF( ln_crban ) THEN ; zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro)
231	ELSE ; zhsro(:,:) = hsro
232	ENDIF
233
234	! Compute shear and dissipation rate
235	#if defined key_z_first
236	DO jj = 2, jpjm1
237	DO ji = 2, jpim1
238	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
239	#else
240	DO jk = 2, jpkm1
241	DO jj = 2, jpjm1
242	DO ji = fs_2, fs_jpim1 ! vector opt.
243	#endif
244	avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) &
245	& * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) &
246	& / ( fse3uw_n(ji,jj,jk) &
247	& * fse3uw_b(ji,jj,jk) )
248	avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) &
249	& * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) &
250	& / ( fse3vw_n(ji,jj,jk) &
251	& * fse3vw_b(ji,jj,jk) )
252	eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk)
253	END DO
254	END DO
255	END DO
256	!
257	! Lateral boundary conditions (avmu,avmv) (sign unchanged)
258	CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. )
259
260	! Save tke at before time step
261	eb (:,:,:) = en (:,:,:)
262	mxlb(:,:,:) = mxln(:,:,:)
263
264	IF( nn_clos == 0 ) THEN ! Mellor-Yamada
265	#if defined key_z_first
266	DO jj = 2, jpjm1
267	DO ji = 2, jpim1
268	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
269	#else
270	DO jk = 2, jpkm1
271	DO jj = 2, jpjm1
272	DO ji = fs_2, fs_jpim1 ! vector opt.
273	#endif
274	zup = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1)
275	zdown = vkarmn * fsdepw(ji,jj,jk) * ( -fsdepw(ji,jj,jk) + fsdepw(ji,jj,mbkt(ji,jj)+1) )
276	zcoef = ( zup / MAX( zdown, rsmall ) )
277	zwall (ji,jj,jk) = ( 1._wp + re2 * zcoefzcoef ) tmask(ji,jj,jk)
278	END DO
279	END DO
280	END DO
281	ENDIF
282
283	!!---------------------------------!!
284	!! Equation to prognostic k !!
285	!!---------------------------------!!
286	!
287	! Now Turbulent kinetic energy (output in en)
288	! -------------------------------
289	! Resolution of a tridiagonal linear system by a "methode de chasse"
290	! computation from level 2 to jpkm1 (e(1) computed after and e(jpk)=0 ).
291	! The surface boundary condition are set after
292	! The bottom boundary condition are also set after. In standard e(bottom)=0.
293	! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal
294	! Warning : after this step, en : right hand side of the matrix
295
296	#if defined key_z_first
297	DO jj = 2, jpjm1
298	DO ji = 2, jpim1
299	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
300	#else
301	DO jk = 2, jpkm1
302	DO jj = 2, jpjm1
303	DO ji = fs_2, fs_jpim1 ! vector opt.
304	#endif
305	!
306	! shear prod. at w-point weightened by mask
307	shear(ji,jj,jk) = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &
308	& + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
309	!
310	! stratif. destruction
311	buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk)
312	!
313	! shear prod. - stratif. destruction
314	diss = eps(ji,jj,jk)
315	!
316	dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy ) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
317	!
318	zesh2 = dir(shear(ji,jj,jk)+buoy)+(1._wp-dir)shear(ji,jj,jk) ! production term
319	zdiss = dir(diss/en(ji,jj,jk)) +(1._wp-dir)(diss-buoy)/en(ji,jj,jk) ! dissipation term
320	!
321	! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0
322	! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style)
323	! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries.
324	! Otherwise, this should be rsc_psi/rsc_psi0
325	IF( ln_sigpsi ) THEN
326	zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1.
327	zwall_psi(ji,jj,jk) = rsc_psi / ( zsigpsi * rsc_psi + &
328	(1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) )
329	ELSE
330	zwall_psi(ji,jj,jk) = 1._wp
331	ENDIF
332	!
333	! building the matrix
334	zcof = rfact_tke * tmask(ji,jj,jk)
335	!
336	! lower diagonal
337	z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) &
338	& / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )
339	!
340	! upper diagonal
341	z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) &
342	& / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) )
343	!
344	! diagonal
345	z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) &
346	& + rdt * zdiss * tmask(ji,jj,jk)
347	!
348	! right hand side in en
349	en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk)
350	END DO
351	END DO
352	END DO
353	!
354	#if defined key_z_first
355	DO jj = 1, jpj, 1
356	DO ji = 1, jpi, 1
357	z_elem_b(ji,jj,mbkmax(ji,jj)) = 1._wp
358	END DO
359	END DO
360	#else
361	z_elem_b(:,:,jpk) = 1._wp
362	#endif
363	!
364	! Set surface condition on zwall_psi (1 at the bottom)
365	IF( ln_sigpsi ) THEN
366	zcoef = rsc_psi / rsc_psi0
367	DO jj = 2, jpjm1
368	DO ji = fs_2, fs_jpim1 ! vector opt.
369	zwall_psi(ji,jj,1) = zcoef
370	END DO
371	END DO
372	ENDIF
373
374	! Surface boundary condition on tke
375	! ---------------------------------
376	!
377	SELECT CASE ( nn_tkebc_surf )
378	!
379	CASE ( 0 ) ! Dirichlet case
380	!
381	IF (ln_crban) THEN ! Wave induced mixing case
382	! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2
383	! ! balance between the production and the dissipation terms including the wave effect
384	en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin )
385	z_elem_a(:,:,1) = en(:,:,1)
386	z_elem_c(:,:,1) = 0._wp
387	z_elem_b(:,:,1) = 1._wp
388	!
389	! one level below
390	en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin )
391	z_elem_a(:,:,2) = 0._wp
392	z_elem_c(:,:,2) = 0._wp
393	z_elem_b(:,:,2) = 1._wp
394	!
395	ELSE ! No wave induced mixing case
396	! ! en(1) = u^2/C0^2 & l(1) = Kzs
397	! ! balance between the production and the dissipation terms
398	en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin )
399	z_elem_a(:,:,1) = en(:,:,1)
400	z_elem_c(:,:,1) = 0._wp
401	z_elem_b(:,:,1) = 1._wp
402	!
403	! one level below
404	en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin )
405	z_elem_a(:,:,2) = 0._wp
406	z_elem_c(:,:,2) = 0._wp
407	z_elem_b(:,:,2) = 1._wp
408	!
409	ENDIF
410	!
411	CASE ( 1 ) ! Neumann boundary condition on d(e)/dz
412	!
413	IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw
414	!
415	! Dirichlet conditions at k=1 (Cosmetic)
416	en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin )
417	z_elem_a(:,:,1) = en(:,:,1)
418	z_elem_c(:,:,1) = 0._wp
419	z_elem_b(:,:,1) = 1._wp
420	! at k=2, set de/dz=Fw
421	z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
422	z_elem_a(:,:,2) = 0._wp
423	zflxs(:,:) = rsbc_tke3 * ustars2(:,:)*1.5_wp ( (zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )*(1.5ra_sf)
424	en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
425	!
426	ELSE ! No wave induced mixing case: d(e)/dz=0.
427	!
428	! Dirichlet conditions at k=1 (Cosmetic)
429	en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin )
430	z_elem_a(:,:,1) = en(:,:,1)
431	z_elem_c(:,:,1) = 0._wp
432	z_elem_b(:,:,1) = 1._wp
433	! at k=2 set de/dz=0.:
434	z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
435	z_elem_a(:,:,2) = 0._wp
436	!
437	ENDIF
438	!
439	END SELECT
440
441	! Bottom boundary condition on tke
442	! --------------------------------
443	!
444	SELECT CASE ( nn_tkebc_bot )
445	!
446	CASE ( 0 ) ! Dirichlet
447	! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin
448	! ! Balance between the production and the dissipation terms
449	!CDIR NOVERRCHK
450	DO jj = 2, jpjm1
451	!CDIR NOVERRCHK
452	DO ji = fs_2, fs_jpim1 ! vector opt.
453	ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point
454	ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1
455	!
456	! Bottom level Dirichlet condition:
457	z_elem_a(ji,jj,ibot ) = 0._wp
458	z_elem_c(ji,jj,ibot ) = 0._wp
459	z_elem_b(ji,jj,ibot ) = 1._wp
460	en(ji,jj,ibot ) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
461	!
462	! Just above last level, Dirichlet condition again
463	z_elem_a(ji,jj,ibotm1) = 0._wp
464	z_elem_c(ji,jj,ibotm1) = 0._wp
465	z_elem_b(ji,jj,ibotm1) = 1._wp
466	en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
467	END DO
468	END DO
469	!
470	CASE ( 1 ) ! Neumman boundary condition
471	!
472	!CDIR NOVERRCHK
473	DO jj = 2, jpjm1
474	!CDIR NOVERRCHK
475	DO ji = fs_2, fs_jpim1 ! vector opt.
476	ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point
477	ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1
478	!
479	! Bottom level Dirichlet condition:
480	z_elem_a(ji,jj,ibot) = 0._wp
481	z_elem_c(ji,jj,ibot) = 0._wp
482	z_elem_b(ji,jj,ibot) = 1._wp
483	en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
484	!
485	! Just above last level: Neumann condition
486	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
487	z_elem_c(ji,jj,ibotm1) = 0._wp
488	END DO
489	END DO
490	!
491	END SELECT
492
493	! Matrix inversion (en prescribed at surface and the bottom)
494	! ----------------------------------------------------------
495	!
496	#if defined key_z_first
497	DO jj = 2, jpjm1
498	DO ji = 2, jpim1
499	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
500	#else
501	DO jk = 2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
502	DO jj = 2, jpjm1
503	DO ji = fs_2, fs_jpim1 ! vector opt.
504	#endif
505	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)
506	END DO
507	END DO
508	END DO
509	#if defined key_z_first
510	DO jj = 2, jpjm1
511	DO ji = 2, jpim1
512	DO jk = 2, mbkmax(ji,jj) ! jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
513	#else
514	DO jk = 2, jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
515	DO jj = 2, jpjm1
516	DO ji = fs_2, fs_jpim1 ! vector opt.
517	#endif
518	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)
519	END DO
520	END DO
521	END DO
522	#if defined key_z_first
523	DO jj = 2, jpjm1
524	DO ji = 2, jpim1
525	DO jk = mbkmax(ji,jj)-1, 2, -1 ! jpk-1, 2, -1 ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
526	#else
527	DO jk = jpk-1, 2, -1 ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
528	DO jj = 2, jpjm1
529	DO ji = fs_2, fs_jpim1 ! vector opt.
530	#endif
531	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)
532	END DO
533	END DO
534	END DO
535	! ! set the minimum value of tke
536	en(:,:,:) = MAX( en(:,:,:), rn_emin )
537
538	!!----------------------------------------!!
539	!! Solve prognostic equation for psi !!
540	!!----------------------------------------!!
541
542	! Set psi to previous time step value
543	!
544	SELECT CASE ( nn_clos )
545	!
546	CASE( 0 ) ! k-kl (Mellor-Yamada)
547	#if defined key_z_first
548	DO jj = 2, jpjm1
549	DO ji = 2, jpim1
550	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
551	#else
552	DO jk = 2, jpkm1
553	DO jj = 2, jpjm1
554	DO ji = fs_2, fs_jpim1 ! vector opt.
555	#endif
556	psi(ji,jj,jk) = en(ji,jj,jk) * mxln(ji,jj,jk)
557	END DO
558	END DO
559	END DO
560	!
561	CASE( 1 ) ! k-eps
562	#if defined key_z_first
563	DO jj = 2, jpjm1
564	DO ji = 2, jpim1
565	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
566	#else
567	DO jk = 2, jpkm1
568	DO jj = 2, jpjm1
569	DO ji = fs_2, fs_jpim1 ! vector opt.
570	#endif
571	psi(ji,jj,jk) = eps(ji,jj,jk)
572	END DO
573	END DO
574	END DO
575	!
576	CASE( 2 ) ! k-w
577	#if defined key_z_first
578	DO jj = 2, jpjm1
579	DO ji = 2, jpim1
580	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
581	#else
582	DO jk = 2, jpkm1
583	DO jj = 2, jpjm1
584	DO ji = fs_2, fs_jpim1 ! vector opt.
585	#endif
586	psi(ji,jj,jk) = SQRT( en(ji,jj,jk) ) / ( rc0 * mxln(ji,jj,jk) )
587	END DO
588	END DO
589	END DO
590	!
591	CASE( 3 ) ! generic
592	#if defined key_z_first
593	DO jj = 2, jpjm1
594	DO ji = 2, jpim1
595	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
596	#else
597	DO jk = 2, jpkm1
598	DO jj = 2, jpjm1
599	DO ji = fs_2, fs_jpim1 ! vector opt.
600	#endif
601	psi(ji,jj,jk) = rc02 * en(ji,jj,jk) * mxln(ji,jj,jk)**rnn
602	END DO
603	END DO
604	END DO
605	!
606	END SELECT
607	!
608	! Now gls (output in psi)
609	! -------------------------------
610	! Resolution of a tridiagonal linear system by a "methode de chasse"
611	! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ).
612	! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal
613	! Warning : after this step, en : right hand side of the matrix
614
615	#if defined key_z_first
616	DO jj = 2, jpjm1
617	DO ji = 2, jpim1
618	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
619	#else
620	DO jk = 2, jpkm1
621	DO jj = 2, jpjm1
622	DO ji = fs_2, fs_jpim1 ! vector opt.
623	#endif
624	!
625	! psi / k
626	zratio = psi(ji,jj,jk) / eb(ji,jj,jk)
627	!
628	! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable)
629	dir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )
630	!
631	rpsi3 = dir * rpsi3m + ( 1._wp - dir ) * rpsi3p
632	!
633	! shear prod. - stratif. destruction
634	prod = rpsi1 * zratio * shear(ji,jj,jk)
635	!
636	! stratif. destruction
637	buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) )
638	!
639	! shear prod. - stratif. destruction
640	diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
641	!
642	dir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
643	!
644	zesh2 = dir * ( prod + buoy ) + (1._wp - dir ) * prod ! production term
645	zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
646	!
647	! building the matrix
648	zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
649	! lower diagonal
650	z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) &
651	& / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )
652	! upper diagonal
653	z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) &
654	& / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) )
655	! diagonal
656	z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) &
657	& + rdt * zdiss * tmask(ji,jj,jk)
658	!
659	! right hand side in psi
660	psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk)
661	END DO
662	END DO
663	END DO
664	!
665	#if defined key_z_first
666	DO jj = 1, jpj, 1
667	DO ji = 1, jpi, 1
668	z_elem_b(ji,jj,mbkmax(ji,jj)) = 1.0_wp
669	END DO
670	END DO
671	#else
672	z_elem_b(:,:,jpk) = 1._wp
673	#endif
674
675	! Surface boundary condition on psi
676	! ---------------------------------
677	!
678	SELECT CASE ( nn_psibc_surf )
679	!
680	CASE ( 0 ) ! Dirichlet boundary conditions
681	!
682	IF( ln_crban ) THEN ! Wave induced mixing case
683	! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2
684	! ! balance between the production and the dissipation terms including the wave effect
685	zdep(:,:) = rl_sf * zhsro(:,:)
686	#if defined key_z_first
687	psi (:,:,1) = rc0*rpp en(:,:,1)*rmm zdep(:,:)*rnn tmask_1(:,:)
688	#else
689	psi (:,:,1) = rc0*rpp en(:,:,1)*rmm zdep(:,:)*rnn tmask(:,:,1)
690	#endif
691	z_elem_a(:,:,1) = psi(:,:,1)
692	z_elem_c(:,:,1) = 0._wp
693	z_elem_b(:,:,1) = 1._wp
694	!
695	! one level below
696	zex1 = (rmm*ra_sf+rnn)
697	zex2 = (rmm*ra_sf)
698	zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))zex1 ) / zhsro(:,:)zex2
699	#if defined key_z_first
700	psi (:,:,2) = rsbc_psi1 * ustars2(:,:)*rmm zdep(:,:) * tmask_1(:,:)
701	#else
702	psi (:,:,2) = rsbc_psi1 * ustars2(:,:)*rmm zdep(:,:) * tmask(:,:,1)
703	#endif
704	z_elem_a(:,:,2) = 0._wp
705	z_elem_c(:,:,2) = 0._wp
706	z_elem_b(:,:,2) = 1._wp
707	!
708	ELSE ! No wave induced mixing case
709	! ! en(1) = u^2/C0^2 & l(1) = Kzs
710	! ! balance between the production and the dissipation terms
711	!
712	zdep(:,:) = vkarmn * zhsro(:,:)
713	#if defined key_z_first
714	psi (:,:,1) = rc0*rpp en(:,:,1)*rmm zdep(:,:)*rnn tmask_1(:,:)
715	#else
716	psi (:,:,1) = rc0*rpp en(:,:,1)*rmm zdep(:,:)*rnn tmask(:,:,1)
717	#endif
718	z_elem_a(:,:,1) = psi(:,:,1)
719	z_elem_c(:,:,1) = 0._wp
720	z_elem_b(:,:,1) = 1._wp
721	!
722	! one level below
723	zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) )
724	#if defined key_z_first
725	psi (:,:,2) = rc0*rpp en(:,:,1)*rmm zdep(:,:)*rnn tmask_1(:,:)
726	#else
727	psi (:,:,2) = rc0*rpp en(:,:,1)*rmm zdep(:,:)*rnn tmask(:,:,1)
728	#endif
729	z_elem_a(:,:,2) = 0._wp
730	z_elem_c(:,:,2) = 0._wp
731	z_elem_b(:,:,2) = 1.
732	!
733	ENDIF
734	!
735	CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz
736	!
737	IF( ln_crban ) THEN ! Wave induced mixing case
738	!
739	zdep(:,:) = rl_sf * zhsro(:,:)
740	#if defined key_z_first
741	psi (:,:,1) = rc0*rpp en(:,:,1)*rmm zdep(:,:)*rnn tmask_1(:,:)
742	#else
743	psi (:,:,1) = rc0*rpp en(:,:,1)*rmm zdep(:,:)*rnn tmask(:,:,1)
744	#endif
745	z_elem_a(:,:,1) = psi(:,:,1)
746	z_elem_c(:,:,1) = 0._wp
747	z_elem_b(:,:,1) = 1._wp
748	!
749	! Neumann condition at k=2
750	z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
751	z_elem_a(:,:,2) = 0._wp
752	!
753	! Set psi vertical flux at the surface:
754	zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))*(rmmra_sf+rnn-1._wp) / zhsro(:,:)*(rmmra_sf)
755	zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)avm(:,:,1) + zwall_psi(:,:,2)avm(:,:,2) ) &
756	& * en(:,:,1)*rmm zdep
757	psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
758	!
759	ELSE ! No wave induced mixing
760	!
761	zdep(:,:) = vkarmn * zhsro(:,:)
762	#if defined key_z_first
763	psi (:,:,1) = rc0*rpp en(:,:,1)*rmm zdep(:,:)*rnn tmask_1(:,:)
764	#else
765	psi (:,:,1) = rc0*rpp en(:,:,1)*rmm zdep(:,:)*rnn tmask(:,:,1)
766	#endif
767	z_elem_a(:,:,1) = psi(:,:,1)
768	z_elem_c(:,:,1) = 0._wp
769	z_elem_b(:,:,1) = 1._wp
770	!
771	! Neumann condition at k=2
772	z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
773	z_elem_a(ji,jj,2) = 0._wp
774	!
775	! Set psi vertical flux at the surface:
776	zdep(:,:) = zhsro(:,:) + fsdept(:,:,1)
777	zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)*rmm zdep**(rnn-1._wp)
778	psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
779	!
780	ENDIF
781	!
782	END SELECT
783
784	! Bottom boundary condition on psi
785	! --------------------------------
786	!
787	SELECT CASE ( nn_psibc_bot )
788	!
789	CASE ( 0 ) ! Dirichlet
790	! ! en(ibot) = u^2 / Co2 and mxln(ibot) = vkarmn hbro
791	! ! Balance between the production and the dissipation terms
792	!CDIR NOVERRCHK
793	DO jj = 2, jpjm1
794	!CDIR NOVERRCHK
795	DO ji = fs_2, fs_jpim1 ! vector opt.
796	ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point
797	ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1
798	zdep(ji,jj) = vkarmn * hbro
799	psi (ji,jj,ibot) = rc0*rpp en(ji,jj,ibot)*rmm zdep(ji,jj)**rnn
800	z_elem_a(ji,jj,ibot) = 0._wp
801	z_elem_c(ji,jj,ibot) = 0._wp
802	z_elem_b(ji,jj,ibot) = 1._wp
803	!
804	! Just above last level, Dirichlet condition again (GOTM like)
805	zdep(ji,jj) = vkarmn * ( hbro + fse3t(ji,jj,ibotm1) )
806	psi (ji,jj,ibotm1) = rc0*rpp en(ji,jj,ibot )*rmm zdep(ji,jj)**rnn
807	z_elem_a(ji,jj,ibotm1) = 0._wp
808	z_elem_c(ji,jj,ibotm1) = 0._wp
809	z_elem_b(ji,jj,ibotm1) = 1._wp
810	END DO
811	END DO
812	!
813	CASE ( 1 ) ! Neumman boundary condition
814	!
815	!CDIR NOVERRCHK
816	DO jj = 2, jpjm1
817	!CDIR NOVERRCHK
818	DO ji = fs_2, fs_jpim1 ! vector opt.
819	ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point
820	ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1
821	!
822	! Bottom level Dirichlet condition:
823	zdep(ji,jj) = vkarmn * hbro
824	psi (ji,jj,ibot) = rc0*rpp en(ji,jj,ibot)*rmm zdep(ji,jj)**rnn
825	!
826	z_elem_a(ji,jj,ibot) = 0._wp
827	z_elem_c(ji,jj,ibot) = 0._wp
828	z_elem_b(ji,jj,ibot) = 1._wp
829	!
830	! Just above last level: Neumann condition with flux injection
831	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
832	z_elem_c(ji,jj,ibotm1) = 0.
833	!
834	! Set psi vertical flux at the bottom:
835	zdep(ji,jj) = hbro + 0.5_wp*fse3t(ji,jj,ibotm1)
836	zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) &
837	& * (0.5_wp(en(ji,jj,ibot)+en(ji,jj,ibotm1)))rmm zdep(ji,jj)**(rnn-1._wp)
838	psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1)
839	END DO
840	END DO
841	!
842	END SELECT
843
844	! Matrix inversion
845	! ----------------
846	!
847	#if defined key_z_first
848	DO jj = 2, jpjm1
849	DO ji = 2, jpim1
850	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
851	#else
852	DO jk = 2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
853	DO jj = 2, jpjm1
854	DO ji = fs_2, fs_jpim1 ! vector opt.
855	#endif
856	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)
857	END DO
858	END DO
859	END DO
860	#if defined key_z_first
861	DO jj = 2, jpjm1
862	DO ji = 2, jpim1
863	DO jk = 2, mbkmax(ji,jj) ! jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
864	#else
865	DO jk = 2, jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
866	DO jj = 2, jpjm1
867	DO ji = fs_2, fs_jpim1 ! vector opt.
868	#endif
869	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)
870	END DO
871	END DO
872	END DO
873	#if defined key_z_first
874	DO jj = 2, jpjm1
875	DO ji = 2, jpim1
876	DO jk = mbkmax(ji,jj)-1, 2, -1 ! jpk-1, 2, -1 ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
877	#else
878	DO jk = jpk-1, 2, -1 ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
879	DO jj = 2, jpjm1
880	DO ji = fs_2, fs_jpim1 ! vector opt.
881	#endif
882	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)
883	END DO
884	END DO
885	END DO
886
887	! Set dissipation
888	!----------------
889
890	SELECT CASE ( nn_clos )
891	!
892	CASE( 0 ) ! k-kl (Mellor-Yamada)
893	#if defined key_z_first
894	DO jj = 2, jpjm1
895	DO ji = 2, jpim1
896	DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
897	#else
898	DO jk = 1, jpkm1
899	DO jj = 2, jpjm1
900	DO ji = fs_2, fs_jpim1 ! vector opt.
901	#endif
902	eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / psi(ji,jj,jk)
903	END DO
904	END DO
905	END DO
906	!
907	CASE( 1 ) ! k-eps
908	#if defined key_z_first
909	DO jj = 2, jpjm1
910	DO ji = 2, jpim1
911	DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
912	#else
913	DO jk = 1, jpkm1
914	DO jj = 2, jpjm1
915	DO ji = fs_2, fs_jpim1 ! vector opt.
916	#endif
917	eps(ji,jj,jk) = psi(ji,jj,jk)
918	END DO
919	END DO
920	END DO
921	!
922	CASE( 2 ) ! k-w
923	#if defined key_z_first
924	DO jj = 2, jpjm1
925	DO ji = 2, jpim1
926	DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
927	#else
928	DO jk = 1, jpkm1
929	DO jj = 2, jpjm1
930	DO ji = fs_2, fs_jpim1 ! vector opt.
931	#endif
932	eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)
933	END DO
934	END DO
935	END DO
936	!
937	CASE( 3 ) ! generic
938	zcoef = rc0**( 3._wp + rpp/rnn )
939	zex1 = ( 1.5_wp + rmm/rnn )
940	zex2 = -1._wp / rnn
941	#if defined key_z_first
942	DO jj = 2, jpjm1
943	DO ji = 2, jpim1
944	DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
945	#else
946	DO jk = 1, jpkm1
947	DO jj = 2, jpjm1
948	DO ji = fs_2, fs_jpim1 ! vector opt.
949	#endif
950	eps(ji,jj,jk) = zcoef * en(ji,jj,jk)*zex1 psi(ji,jj,jk)**zex2
951	END DO
952	END DO
953	END DO
954	!
955	END SELECT
956
957	! Limit dissipation rate under stable stratification
958	! --------------------------------------------------
959	#if defined key_z_first
960	DO jj = 2, jpjm1
961	DO ji = 2, jpim1
962	DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1 ! Note that this set boundary conditions on mxln at the same time
963	#else
964	DO jk = 1, jpkm1 ! Note that this set boundary conditions on mxln at the same time
965	DO jj = 2, jpjm1
966	DO ji = fs_2, fs_jpim1 ! vector opt.
967	#endif
968	! limitation
969	eps(ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin )
970	mxln(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
971	! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
972	zrn2 = MAX( rn2(ji,jj,jk), rsmall )
973	mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) )
974	END DO
975	END DO
976	END DO
977
978	!
979	! Stability function and vertical viscosity and diffusivity
980	! ---------------------------------------------------------
981	!
982	SELECT CASE ( nn_stab_func )
983	!
984	CASE ( 0 , 1 ) ! Galperin or Kantha-Clayson stability functions
985	#if defined key_z_first
986	DO jj = 2, jpjm1
987	DO ji = 2, jpim1
988	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
989	#else
990	DO jk = 2, jpkm1
991	DO jj = 2, jpjm1
992	DO ji = fs_2, fs_jpim1 ! vector opt.
993	#endif
994	! zcof = l²/q²
995	zcof = mxlb(ji,jj,jk) * mxlb(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
996	! Gh = -N²l²/q²
997	gh = - rn2(ji,jj,jk) * zcof
998	gh = MIN( gh, rgh0 )
999	gh = MAX( gh, rghmin )
1000	! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
1001	sh = ra2( 1._wp-6._wpra1/rb1 ) / ( 1.-3.ra2gh(6.ra1+rb2*( 1._wp-rc3 ) ) )
1002	sm = ( rb1*(-1._wp/3._wp) + ( 18._wpra1ra1 + 9._wpra1ra2(1._wp-rc2) )shgh ) / (1._wp-9._wpra1ra2*gh)
1003	!
1004	! Store stability function in avmu and avmv
1005	avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
1006	avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
1007	END DO
1008	END DO
1009	END DO
1010	!
1011	CASE ( 2, 3 ) ! Canuto stability functions
1012	#if defined key_z_first
1013	DO jj = 2, jpjm1
1014	DO ji = 2, jpim1
1015	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
1016	#else
1017	DO jk = 2, jpkm1
1018	DO jj = 2, jpjm1
1019	DO ji = fs_2, fs_jpim1 ! vector opt.
1020	#endif
1021	! zcof = l²/q²
1022	zcof = mxlb(ji,jj,jk)mxlb(ji,jj,jk) / ( 2._wp eb(ji,jj,jk) )
1023	! Gh = -N²l²/q²
1024	gh = - rn2(ji,jj,jk) * zcof
1025	gh = MIN( gh, rgh0 )
1026	gh = MAX( gh, rghmin )
1027	gh = gh * rf6
1028	! Gm = M²l²/q² Shear number
1029	shr = shear(ji,jj,jk) / MAX( avm(ji,jj,jk), rsmall )
1030	gm = MAX( shr * zcof , 1.e-10 )
1031	gm = gm * rf6
1032	gm = MIN ( (rd0 - rd1gh + rd3ghgh) / (rd2-rd4gh) , gm )
1033	! Stability functions from Canuto
1034	rcff = rd0 - rd1gh +rd2gm + rd3ghgh - rd4ghgm + rd5gmgm
1035	sm = (rs0 - rs1gh + rs2gm) / rcff
1036	sh = (rs4 - rs5gh + rs6gm) / rcff
1037	!
1038	! Store stability function in avmu and avmv
1039	avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
1040	avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
1041	END DO
1042	END DO
1043	END DO
1044	!
1045	END SELECT
1046
1047	! Boundary conditions on stability functions for momentum (Neumann):
1048	! Lines below are useless if GOTM style Dirichlet conditions are used
1049	zcoef = rcm_sf / SQRT( 2._wp )
1050	DO jj = 2, jpjm1
1051	DO ji = fs_2, fs_jpim1 ! vector opt.
1052	avmv(ji,jj,1) = zcoef
1053	END DO
1054	END DO
1055	zcoef = rc0 / SQRT( 2._wp )
1056	DO jj = 2, jpjm1
1057	DO ji = fs_2, fs_jpim1 ! vector opt.
1058	avmv(ji,jj,mbkt(ji,jj)+1) = zcoef
1059	END DO
1060	END DO
1061
1062	! Compute diffusivities/viscosities
1063	! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
1064	#if defined key_z_first
1065	DO jj = 2, jpjm1
1066	DO ji = 2, jpim1
1067	DO jk = 1, mbkmax(ji,jj) ! jpk
1068	#else
1069	DO jk = 1, jpk
1070	DO jj = 2, jpjm1
1071	DO ji = fs_2, fs_jpim1 ! vector opt.
1072	#endif
1073	zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * mxln(ji,jj,jk)
1074	zav = zsqen * avmu(ji,jj,jk)
1075	avt(ji,jj,jk) = MAX( zav, avtb(jk) )*tmask(ji,jj,jk) ! apply mask for zdfmxl routine
1076	zav = zsqen * avmv(ji,jj,jk)
1077	avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom
1078	END DO
1079	END DO
1080	END DO
1081	!
1082	! Lateral boundary conditions (sign unchanged)
1083	avt(:,:,1) = 0._wp
1084	CALL lbc_lnk( avm, 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. )
1085
1086	#if defined key_z_first
1087	DO jj = 2, jpjm1
1088	DO ji = 2, jpim1
1089	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1 !* vertical eddy viscosity at u- and v-points
1090	#else
1091	DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points
1092	DO jj = 2, jpjm1
1093	DO ji = fs_2, fs_jpim1 ! vector opt.
1094	#endif
1095	avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk)
1096	avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk)
1097	END DO
1098	END DO
1099	END DO
1100	avmu(:,:,1) = 0._wp ; avmv(:,:,1) = 0._wp ! set surface to zero
1101	CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions
1102
1103	IF(ln_ctl) THEN
1104	CALL prt_ctl( tab3d_1=en , clinfo1=' gls - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
1105	CALL prt_ctl( tab3d_1=avmu, clinfo1=' gls - u: ', mask1=umask, &
1106	& tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
1107	ENDIF
1108	!
1109	IF( wrk_not_released(2, 1,2,3) .OR. &
1110	wrk_not_released(3, 1,2,3,4,5) ) CALL ctl_stop('zdf_gls: failed to release workspace arrays')
1111	!
1112	CALL timing_stop('zdf_gls','section')
1113	!
1114	END SUBROUTINE zdf_gls
1115
1116	!! * Reset control of array index permutation
1117	!FTRANS CLEAR
1118	# include "oce_ftrans.h90"
1119	# include "dom_oce_ftrans.h90"
1120	# include "domvvl_ftrans.h90"
1121	# include "zdf_oce_ftrans.h90"
1122	# include "sbc_oce_ftrans.h90"
1123	!! DCSE_NEMO: private module variables do not need their own directives file
1124	!FTRANS en mxln zwall :I :I :z
1125
1126	SUBROUTINE zdf_gls_init
1127	!!----------------------------------------------------------------------
1128	!! * ROUTINE zdf_gls_init *
1129	!!
1130	!! ** Purpose : Initialization of the vertical eddy diffivity and
1131	!! viscosity when using a gls turbulent closure scheme
1132	!!
1133	!! ** Method : Read the namzdf_gls namelist and check the parameters
1134	!! called at the first timestep (nit000)
1135	!!
1136	!! ** input : Namlist namzdf_gls
1137	!!
1138	!! ** Action : Increase by 1 the nstop flag is setting problem encounter
1139	!!
1140	!!----------------------------------------------------------------------
1141	USE dynzdf_exp
1142	USE trazdf_exp
1143	!
1144	INTEGER :: ji, jj, jk ! dummy loop indices
1145	REAL(wp):: zcr ! local scalar
1146	!!
1147	NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
1148	& rn_clim_galp, ln_crban, ln_sigpsi, &
1149	& rn_crban, rn_charn, &
1150	& nn_tkebc_surf, nn_tkebc_bot, &
1151	& nn_psibc_surf, nn_psibc_bot, &
1152	& nn_stab_func, nn_clos
1153	!!----------------------------------------------------------
1154
1155	REWIND( numnam ) !* Read Namelist namzdf_gls
1156	READ ( numnam, namzdf_gls )
1157
1158	IF(lwp) THEN !* Control print
1159	WRITE(numout,*)
1160	WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme'
1161	WRITE(numout,*) '~~~~~~~~~~~~'
1162	WRITE(numout,*) ' Namelist namzdf_gls : set gls mixing parameters'
1163	WRITE(numout,*) ' minimum value of en rn_emin = ', rn_emin
1164	WRITE(numout,*) ' minimum value of eps rn_epsmin = ', rn_epsmin
1165	WRITE(numout,*) ' Limit dissipation rate under stable stratif. ln_length_lim = ', ln_length_lim
1166	WRITE(numout,*) ' Galperin limit (Standard: 0.53, Holt: 0.26) rn_clim_galp = ', rn_clim_galp
1167	WRITE(numout,*) ' TKE Surface boundary condition nn_tkebc_surf = ', nn_tkebc_surf
1168	WRITE(numout,*) ' TKE Bottom boundary condition nn_tkebc_bot = ', nn_tkebc_bot
1169	WRITE(numout,*) ' PSI Surface boundary condition nn_psibc_surf = ', nn_psibc_surf
1170	WRITE(numout,*) ' PSI Bottom boundary condition nn_psibc_bot = ', nn_psibc_bot
1171	WRITE(numout,*) ' Craig and Banner scheme ln_crban = ', ln_crban
1172	WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi
1173	WRITE(numout,*) ' Craig and Banner coefficient rn_crban = ', rn_crban
1174	WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn
1175	WRITE(numout,*) ' Stability functions nn_stab_func = ', nn_stab_func
1176	WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos
1177	WRITE(numout,*) ' Hard coded parameters'
1178	WRITE(numout,*) ' Surface roughness (m) hsro = ', hsro
1179	WRITE(numout,*) ' Bottom roughness (m) hbro = ', hbro
1180	ENDIF
1181
1182	! !* allocate gls arrays
1183	IF( zdf_gls_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
1184
1185	! !* Check of some namelist values
1186	IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' )
1187	IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' )
1188	IF( nn_tkebc_bot < 0 .OR. nn_tkebc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' )
1189	IF( nn_psibc_bot < 0 .OR. nn_psibc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' )
1190	IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' )
1191	IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' )
1192
1193	SELECT CASE ( nn_clos ) !* set the parameters for the chosen closure
1194	!
1195	CASE( 0 ) ! k-kl (Mellor-Yamada)
1196	!
1197	IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada'
1198	rpp = 0._wp
1199	rmm = 1._wp
1200	rnn = 1._wp
1201	rsc_tke = 1.96_wp
1202	rsc_psi = 1.96_wp
1203	rpsi1 = 0.9_wp
1204	rpsi3p = 1._wp
1205	rpsi2 = 0.5_wp
1206	!
1207	SELECT CASE ( nn_stab_func )
1208	CASE( 0, 1 ) ; rpsi3m = 2.53_wp ! G88 or KC stability functions
1209	CASE( 2 ) ; rpsi3m = 2.38_wp ! Canuto A stability functions
1210	CASE( 3 ) ; rpsi3m = 2.38 ! Canuto B stability functions (caution : constant not identified)
1211	END SELECT
1212	!
1213	CASE( 1 ) ! k-eps
1214	!
1215	IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps'
1216	rpp = 3._wp
1217	rmm = 1.5_wp
1218	rnn = -1._wp
1219	rsc_tke = 1._wp
1220	rsc_psi = 1.3_wp ! Schmidt number for psi
1221	rpsi1 = 1.44_wp
1222	rpsi3p = 1._wp
1223	rpsi2 = 1.92_wp
1224	!
1225	SELECT CASE ( nn_stab_func )
1226	CASE( 0, 1 ) ; rpsi3m = -0.52_wp ! G88 or KC stability functions
1227	CASE( 2 ) ; rpsi3m = -0.629_wp ! Canuto A stability functions
1228	CASE( 3 ) ; rpsi3m = -0.566 ! Canuto B stability functions
1229	END SELECT
1230	!
1231	CASE( 2 ) ! k-omega
1232	!
1233	IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega'
1234	rpp = -1._wp
1235	rmm = 0.5_wp
1236	rnn = -1._wp
1237	rsc_tke = 2._wp
1238	rsc_psi = 2._wp
1239	rpsi1 = 0.555_wp
1240	rpsi3p = 1._wp
1241	rpsi2 = 0.833_wp
1242	!
1243	SELECT CASE ( nn_stab_func )
1244	CASE( 0, 1 ) ; rpsi3m = -0.58_wp ! G88 or KC stability functions
1245	CASE( 2 ) ; rpsi3m = -0.64_wp ! Canuto A stability functions
1246	CASE( 3 ) ; rpsi3m = -0.64_wp ! Canuto B stability functions caution : constant not identified)
1247	END SELECT
1248	!
1249	CASE( 3 ) ! generic
1250	!
1251	IF(lwp) WRITE(numout,*) 'The choosen closure is generic'
1252	rpp = 2._wp
1253	rmm = 1._wp
1254	rnn = -0.67_wp
1255	rsc_tke = 0.8_wp
1256	rsc_psi = 1.07_wp
1257	rpsi1 = 1._wp
1258	rpsi3p = 1._wp
1259	rpsi2 = 1.22_wp
1260	!
1261	SELECT CASE ( nn_stab_func )
1262	CASE( 0, 1 ) ; rpsi3m = 0.1_wp ! G88 or KC stability functions
1263	CASE( 2 ) ; rpsi3m = 0.05_wp ! Canuto A stability functions
1264	CASE( 3 ) ; rpsi3m = 0.05_wp ! Canuto B stability functions caution : constant not identified)
1265	END SELECT
1266	!
1267	END SELECT
1268
1269	!
1270	SELECT CASE ( nn_stab_func ) !* set the parameters of the stability functions
1271	!
1272	CASE ( 0 ) ! Galperin stability functions
1273	!
1274	IF(lwp) WRITE(numout,*) 'Stability functions from Galperin'
1275	rc2 = 0._wp
1276	rc3 = 0._wp
1277	rc_diff = 1._wp
1278	rc0 = 0.5544_wp
1279	rcm_sf = 0.9884_wp
1280	rghmin = -0.28_wp
1281	rgh0 = 0.0233_wp
1282	rghcri = 0.02_wp
1283	!
1284	CASE ( 1 ) ! Kantha-Clayson stability functions
1285	!
1286	IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson'
1287	rc2 = 0.7_wp
1288	rc3 = 0.2_wp
1289	rc_diff = 1._wp
1290	rc0 = 0.5544_wp
1291	rcm_sf = 0.9884_wp
1292	rghmin = -0.28_wp
1293	rgh0 = 0.0233_wp
1294	rghcri = 0.02_wp
1295	!
1296	CASE ( 2 ) ! Canuto A stability functions
1297	!
1298	IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A'
1299	rs0 = 1.5_wp * rl1 * rl5*rl5
1300	rs1 = -rl4(rl6+rl7) + 2._wprl4rl5(rl1-(1._wp/3._wp)rl2-rl3) + 1.5_wprl1rl5rl8
1301	rs2 = -(3._wp/8._wp) * rl1(rl6rl6-rl7*rl7)
1302	rs4 = 2._wp * rl5
1303	rs5 = 2._wp * rl4
1304	rs6 = (2._wp/3._wp) * rl5 * ( 3._wprl3rl3 - rl2rl2 ) - 0.5_wp rl5rl1 (3._wp*rl3-rl2) &
1305	& + 0.75_wp * rl1 * ( rl6 - rl7 )
1306	rd0 = 3._wp * rl5*rl5
1307	rd1 = rl5 * ( 7._wprl4 + 3._wprl8 )
1308	rd2 = rl5rl5 ( 3._wprl3rl3 - rl2rl2 ) - 0.75_wp(rl6rl6 - rl7rl7 )
1309	rd3 = rl4 * ( 4._wprl4 + 3._wprl8)
1310	rd4 = rl4 * ( rl2 * rl6 - 3._wprl3rl7 - rl5(rl2rl2 - rl3rl3 ) ) + rl5rl8 * ( 3._wprl3rl3 - rl2*rl2 )
1311	rd5 = 0.25_wp * ( rl2rl2 - 3._wp rl3rl3 ) ( rl6rl6 - rl7rl7 )
1312	rc0 = 0.5268_wp
1313	rf6 = 8._wp / (rc0**6._wp)
1314	rc_diff = SQRT(2._wp) / (rc0**3._wp)
1315	rcm_sf = 0.7310_wp
1316	rghmin = -0.28_wp
1317	rgh0 = 0.0329_wp
1318	rghcri = 0.03_wp
1319	!
1320	CASE ( 3 ) ! Canuto B stability functions
1321	!
1322	IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B'
1323	rs0 = 1.5_wp * rm1 * rm5*rm5
1324	rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4rm5(rm1-(1._wp/3._wp)rm2-rm3) + 1.5_wp rm1rm5rm8
1325	rs2 = -(3._wp/8._wp) * rm1 * (rm6rm6-rm7rm7 )
1326	rs4 = 2._wp * rm5
1327	rs5 = 2._wp * rm4
1328	rs6 = (2._wp/3._wp) * rm5 * (3._wprm3rm3-rm2rm2) - 0.5_wp rm5rm1(3._wprm3-rm2) + 0.75_wp rm1*(rm6-rm7)
1329	rd0 = 3._wp * rm5*rm5
1330	rd1 = rm5 * (7._wprm4 + 3._wprm8)
1331	rd2 = rm5rm5 (3._wprm3rm3 - rm2rm2) - 0.75_wp (rm6rm6 - rm7rm7)
1332	rd3 = rm4 * ( 4._wprm4 + 3._wprm8 )
1333	rd4 = rm4 * ( rm2rm6 -3._wprm3rm7 - rm5(rm2rm2 - rm3rm3) ) + rm5 * rm8 * ( 3._wprm3rm3 - rm2*rm2 )
1334	rd5 = 0.25_wp * ( rm2rm2 - 3._wprm3rm3 ) ( rm6rm6 - rm7rm7 )
1335	rc0 = 0.5268_wp !! rc0 = 0.5540_wp (Warner ...) to verify !
1336	rf6 = 8._wp / ( rc0**6._wp )
1337	rc_diff = SQRT(2._wp)/(rc0**3.)
1338	rcm_sf = 0.7470_wp
1339	rghmin = -0.28_wp
1340	rgh0 = 0.0444_wp
1341	rghcri = 0.0414_wp
1342	!
1343	END SELECT
1344
1345	! !* Set Schmidt number for psi diffusion in the wave breaking case
1346	! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
1347	! ! or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
1348	IF( ln_sigpsi .AND. ln_crban ) THEN
1349	zcr = SQRT( 1.5_wprsc_tke ) rcm_sf / vkarmn
1350	rsc_psi0 = vkarmnvkarmn / ( rpsi2 rcm_sf*rcm_sf ) &
1351	& * ( rnnrnn - 4._wp/3._wp zcrrnnrmm - 1._wp/3._wp * zcr*rnn &
1352	& + 2._wp/9._wp * rmm * zcrzcr + 4._wp/9._wp zcrzcr rmm*rmm )
1353	ELSE
1354	rsc_psi0 = rsc_psi
1355	ENDIF
1356
1357	! !* Shear free turbulence parameters
1358	!
1359	ra_sf = -4._wp * rnn * SQRT( rsc_tke ) / ( (1._wp+4._wprmm) SQRT( rsc_tke ) &
1360	& - SQRT(rsc_tke + 24._wprsc_psi0rpsi2 ) )
1361	rl_sf = rc0 * SQRT( rc0 / rcm_sf ) &
1362	& * SQRT( ( (1._wp + 4._wprmm + 8._wprmmrmm) rsc_tke &
1363	& + 12._wp * rsc_psi0 * rpsi2 &
1364	& - (1._wp + 4._wprmm) SQRT( rsc_tke(rsc_tke+ 24._wprsc_psi0*rpsi2) ) ) &
1365	& / ( 12._wprnnrnn ) )
1366
1367	!
1368	IF(lwp) THEN !* Control print
1369	WRITE(numout,*)
1370	WRITE(numout,*) 'Limit values'
1371	WRITE(numout,*) '~~~~~~~~~~~~'
1372	WRITE(numout,*) 'Parameter m = ',rmm
1373	WRITE(numout,*) 'Parameter n = ',rnn
1374	WRITE(numout,*) 'Parameter p = ',rpp
1375	WRITE(numout,*) 'rpsi1 = ',rpsi1
1376	WRITE(numout,*) 'rpsi2 = ',rpsi2
1377	WRITE(numout,*) 'rpsi3m = ',rpsi3m
1378	WRITE(numout,*) 'rpsi3p = ',rpsi3p
1379	WRITE(numout,*) 'rsc_tke = ',rsc_tke
1380	WRITE(numout,*) 'rsc_psi = ',rsc_psi
1381	WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0
1382	WRITE(numout,*) 'rc0 = ',rc0
1383	WRITE(numout,*)
1384	WRITE(numout,*) 'Shear free turbulence parameters:'
1385	WRITE(numout,*) 'rcm_sf = ',rcm_sf
1386	WRITE(numout,*) 'ra_sf = ',ra_sf
1387	WRITE(numout,*) 'rl_sf = ',rl_sf
1388	WRITE(numout,*)
1389	ENDIF
1390
1391	! !* Constants initialization
1392	rc02 = rc0 * rc0 ; rc02r = 1. / rc02
1393	rc03 = rc02 * rc0
1394	rc04 = rc03 * rc0
1395	rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp
1396	rsbc_mb = 0.5_wp * (15.8_wprn_crban)*(2._wp/3._wp) ! Surf. bound. cond. from Mellor and Blumberg
1397	rsbc_std = 3.75_wp ! Surf. bound. cond. standard (prod=diss)
1398	rsbc_tke1 = (-rsc_tkern_crban/(rcm_sfra_sfrl_sf))*(2._wp/3._wp) ! k_eps = 53. Dirichlet + Wave breaking
1399	rsbc_tke2 = 0.5_wp / rau0
1400	rsbc_tke3 = rdt * rn_crban ! Neumann + Wave breaking
1401	rsbc_zs = rn_charn / grav ! Charnock formula
1402	rsbc_psi1 = rc0*rpp rsbc_tke1*rmm rl_sf**rnn ! Dirichlet + Wave breaking
1403	rsbc_psi2 = -0.5_wp * rdt * rc0*rpp rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking
1404	rsbc_psi3 = -0.5_wp * rdt * rc0*rpp rl_sf*rnn / rsc_psi (rnn + rmm*ra_sf) ! Neumann + Wave breaking
1405	rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke
1406	rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke
1407
1408	! !* Wall proximity function
1409	zwall (:,:,:) = 1._wp * tmask(:,:,:)
1410
1411	! !* set vertical eddy coef. to the background value
1412	#if defined key_z_first
1413	DO jj = 1, jpj
1414	DO ji = 1, jpi
1415	DO jk = 1, jpk
1416	avt (ji,jj,jk) = avtb(jk) * tmask(ji,jj,jk)
1417	avm (ji,jj,jk) = avmb(jk) * tmask(ji,jj,jk)
1418	avmu(ji,jj,jk) = avmb(jk) * umask(ji,jj,jk)
1419	avmv(ji,jj,jk) = avmb(jk) * vmask(ji,jj,jk)
1420	END DO
1421	END DO
1422	END DO
1423	#else
1424	DO jk = 1, jpk
1425	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
1426	avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
1427	avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
1428	avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
1429	END DO
1430	#endif
1431	!
1432	CALL gls_rst( nit000, 'READ' ) !* read or initialize all required files
1433	!
1434	END SUBROUTINE zdf_gls_init
1435
1436
1437	SUBROUTINE gls_rst( kt, cdrw )
1438	!!---------------------------------------------------------------------
1439	!! * ROUTINE ts_rst *
1440	!!
1441	!! ** Purpose : Read or write TKE file (en) in restart file
1442	!!
1443	!! ** Method : use of IOM library
1444	!! if the restart does not contain TKE, en is either
1445	!! set to rn_emin or recomputed (nn_igls/=0)
1446	!!----------------------------------------------------------------------
1447	INTEGER , INTENT(in) :: kt ! ocean time-step
1448	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1449	!
1450	INTEGER :: jit, jk ! dummy loop indices
1451	INTEGER :: id1, id2, id3, id4, id5, id6
1452	INTEGER :: ji, jj, ikbu, ikbv
1453	REAL(wp):: cbx, cby
1454	!!----------------------------------------------------------------------
1455	!
1456	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
1457	! ! ---------------
1458	IF( ln_rstart ) THEN !* Read the restart file
1459	id1 = iom_varid( numror, 'en' , ldstop = .FALSE. )
1460	id2 = iom_varid( numror, 'avt' , ldstop = .FALSE. )
1461	id3 = iom_varid( numror, 'avm' , ldstop = .FALSE. )
1462	id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
1463	id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
1464	id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. )
1465	!
1466	IF( MIN( id1, id2, id3, id4, id5, id6 ) > 0 ) THEN ! all required arrays exist
1467	CALL iom_get( numror, jpdom_autoglo, 'en' , en )
1468	CALL iom_get( numror, jpdom_autoglo, 'avt' , avt )
1469	CALL iom_get( numror, jpdom_autoglo, 'avm' , avm )
1470	CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu )
1471	CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv )
1472	CALL iom_get( numror, jpdom_autoglo, 'mxln' , mxln )
1473	ELSE
1474	IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
1475	en (:,:,:) = rn_emin
1476	mxln(:,:,:) = 0.001
1477	DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_gls( jit ) ; END DO
1478	ENDIF
1479	ELSE !* Start from rest
1480	IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
1481	en (:,:,:) = rn_emin
1482	mxln(:,:,:) = 0.001
1483	ENDIF
1484	!
1485	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
1486	! ! -------------------
1487	IF(lwp) WRITE(numout,*) '---- gls-rst ----'
1488	CALL iom_rstput( kt, nitrst, numrow, 'en' , en )
1489	CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt )
1490	CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm )
1491	CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu )
1492	CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv )
1493	CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln )
1494	!
1495	ENDIF
1496	!
1497	END SUBROUTINE gls_rst
1498
1499	#else
1500	!!----------------------------------------------------------------------
1501	!! Dummy module : NO TKE scheme
1502	!!----------------------------------------------------------------------
1503	LOGICAL, PUBLIC, PARAMETER :: lk_zdfgls = .FALSE. !: TKE flag
1504	CONTAINS
1505	SUBROUTINE zdf_gls_init ! Empty routine
1506	WRITE(,) 'zdf_gls_init: You should not have seen this print! error?'
1507	END SUBROUTINE zdf_gls_init
1508	SUBROUTINE zdf_gls( kt ) ! Empty routine
1509	WRITE(,) 'zdf_gls: You should not have seen this print! error?', kt
1510	END SUBROUTINE zdf_gls
1511	SUBROUTINE gls_rst( kt, cdrw ) ! Empty routine
1512	INTEGER , INTENT(in) :: kt ! ocean time-step
1513	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1514	WRITE(,) 'gls_rst: You should not have seen this print! error?', kt, cdrw
1515	END SUBROUTINE gls_rst
1516	#endif
1517
1518	!!======================================================================
1519	END MODULE zdfgls
1520

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

Download in other formats: