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/trunk/src/OCE/ZDF – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/ZDF/zdfgls.F90 @ 13497

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

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

Download in other formats: