Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

zdfgls.F90 @ 15608

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

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/ZDF/zdfgls.F90 @ 15608

Download in other formats: