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 @ 14971

Last change on this file since 14971 was 14967, checked in by cbricaud, 3 years ago
introduce limit conditions under ice caps in GLS; see #2593 and #2604
Property svn:keywords set to `Id`
File size: 65.8 KB

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

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

Download in other formats: