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 NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdfgls.F90 @ 12616

Last change on this file since 12616 was 12616, checked in by techene, 4 years ago
all: add e3 substitute (sometimes it requires to add ze3t/u/v/w) and limit precompiled files lines to about 130 character, OCE/ASM/asminc.F90, OCE/DOM/domzgr_substitute.h90, OCE/ISF/isfcpl.F90, OCE/SBC/sbcice_cice, OCE/CRS/crsini.F90 : add key_LF
Property svn:keywords set to `Id`
File size: 54.5 KB

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

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

Download in other formats: