Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

zdfgls.F90

Last change on this file was 15145, checked in by smasson, 3 years ago
trunk: pass all sette tests in debug with nn_hls = 2
Property svn:keywords set to `Id`
File size: 66.1 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	!
[15071]	329	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
	330	zdiag(ji,jj,jpk) = 1._wp
	331	!
	332	! Set surface condition on zwall_psi (1 at the bottom)
	333	zwall_psi(ji,jj, 1 ) = zwall_psi(ji,jj,2)
	334	zwall_psi(ji,jj,jpk) = 1._wp
	335	END_2D
[2048]	336	!
	337	! Surface boundary condition on tke
	338	! ---------------------------------
	339	!
[5109]	340	SELECT CASE ( nn_bc_surf )
[2048]	341	!
[14072]	342	CASE ( 0 ) ! Dirichlet boundary condition (set e at k=1 & 2)
[14834]	343	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
	344	! First level
	345	en (ji,jj,1) = MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))rsbc_tke1)*r2_3 )
	346	zd_lw(ji,jj,1) = en(ji,jj,1)
	347	zd_up(ji,jj,1) = 0._wp
	348	zdiag(ji,jj,1) = 1._wp
	349	!
	350	! One level below
	351	en (ji,jj,2) = MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1 &
	352	& * ((zhsro(ji,jj)+gdepw(ji,jj,2,Kmm)) / zhsro(ji,jj) )*(1.5_wpra_sf) )**r2_3 )
	353	zd_lw(ji,jj,2) = 0._wp
	354	zd_up(ji,jj,2) = 0._wp
	355	zdiag(ji,jj,2) = 1._wp
	356	END_2D
	357	!
[14967]	358	IF( ln_isfcav) THEN ! top boundary (ocean cavity)
	359	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
	360	IF( mikt(ji,jj) > 1 )THEN
	361	itop = mikt(ji,jj) ! k top w-point
	362	itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one
	363	! ! mask at the
	364	! ocean surface
	365	! points
	366	z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
	367	!
	368	! Dirichlet condition applied at:
	369	! top level (itop) & Just below it (itopp1)
	370	zd_lw(ji,jj,itop) = 0._wp ; zd_lw(ji,jj,itopp1) = 0._wp
	371	zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp
	372	zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = 1._wp
	373	en (ji,jj,itop) = z_en ; en (ji,jj,itopp1) = z_en
	374	ENDIF
	375	END_2D
	376	ENDIF
[14834]	377	!
[9019]	378	CASE ( 1 ) ! Neumann boundary condition (set d(e)/dz)
[14834]	379	!
	380	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
	381	! Dirichlet conditions at k=1
	382	en (ji,jj,1) = MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))rsbc_tke1)*r2_3 )
	383	zd_lw(ji,jj,1) = en(ji,jj,1)
	384	zd_up(ji,jj,1) = 0._wp
	385	zdiag(ji,jj,1) = 1._wp
	386	!
	387	! at k=2, set de/dz=Fw
	388	!cbr
	389	! zdiag zd_lw not defined/used on the halo
	390	zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
	391	zd_lw(ji,jj,2) = 0._wp
	392	!
	393	zkar (ji,jj) = (rl_sf + (vkarmn-rl_sf)(1.-EXP(-rtransgdept(ji,jj,1,Kmm)/zhsro(ji,jj)) ))
	394	zflxs(ji,jj) = rsbc_tke2 * (1._wp-zice_fra(ji,jj)) * ustar2_surf(ji,jj)*1.5_wp zkar(ji,jj) &
	395	& * ( ( zhsro(ji,jj)+gdept(ji,jj,1,Kmm) ) / zhsro(ji,jj) )*(1.5_wpra_sf)
[12377]	396	!!gm why not : * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )*(1.5_wpra_sf)
[14834]	397	en(ji,jj,2) = en(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm)
	398	END_2D
	399	!
[14967]	400	IF( ln_isfcav) THEN ! top boundary (ocean cavity)
	401	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
	402	IF( mikt(ji,jj) > 1 )THEN
	403	itop = mikt(ji,jj) ! k top w-point
	404	itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one
	405	! ! mask at the
	406	! ocean surface
	407	! points
	408	z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
	409	!
	410	! Bottom level Dirichlet condition:
	411	! Bottom level (ibot) & Just above it (ibotm1)
	412	! Dirichlet ! Neumann
	413	zd_lw(ji,jj,itop) = 0._wp ! ! Remove zd_up from zdiag
	414	zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1)
	415	zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp
	416	en (ji,jj,itop) = z_en
	417	ENDIF
	418	END_2D
	419	ENDIF
[14834]	420	!
[2048]	421	END SELECT
	422
	423	! Bottom boundary condition on tke
	424	! --------------------------------
	425	!
[5109]	426	SELECT CASE ( nn_bc_bot )
[2048]	427	!
[14072]	428	CASE ( 0 ) ! Dirichlet
[9019]	429	! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
[2397]	430	! ! Balance between the production and the dissipation terms
[14834]	431	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
[9019]	432	!!gm This means that bottom and ocean w-level above have a specified "en" value. Sure ????
	433	!! With thick deep ocean level thickness, this may be quite large, no ???
	434	!! in particular in ocean cavities where top stratification can be large...
[12377]	435	ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point
	436	ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1
	437	!
	438	z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin )
	439	!
[14072]	440	! Dirichlet condition applied at:
	441	! Bottom level (ibot) & Just above it (ibotm1)
[12377]	442	zd_lw(ji,jj,ibot) = 0._wp ; zd_lw(ji,jj,ibotm1) = 0._wp
	443	zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp
	444	zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = 1._wp
	445	en (ji,jj,ibot) = z_en ; en (ji,jj,ibotm1) = z_en
	446	END_2D
[2397]	447	!
[14834]	448	! NOTE: ctl_stop with ln_isfcav when using GLS
[9019]	449	IF( ln_isfcav) THEN ! top boundary (ocean cavity)
[14834]	450	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
[12377]	451	itop = mikt(ji,jj) ! k top w-point
	452	itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one
	453	! ! mask at the ocean surface points
	454	z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
	455	!
[9019]	456	!!gm TO BE VERIFIED !!!
[14072]	457	! Dirichlet condition applied at:
	458	! top level (itop) & Just below it (itopp1)
[12377]	459	zd_lw(ji,jj,itop) = 0._wp ; zd_lw(ji,jj,itopp1) = 0._wp
	460	zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp
	461	zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = 1._wp
	462	en (ji,jj,itop) = z_en ; en (ji,jj,itopp1) = z_en
	463	END_2D
[9019]	464	ENDIF
	465	!
[2048]	466	CASE ( 1 ) ! Neumman boundary condition
[14072]	467	!
[14834]	468	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
[12377]	469	ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point
	470	ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1
	471	!
	472	z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin )
	473	!
	474	! Bottom level Dirichlet condition:
[14072]	475	! Bottom level (ibot) & Just above it (ibotm1)
[12377]	476	! Dirichlet ! Neumann
	477	zd_lw(ji,jj,ibot) = 0._wp ! ! Remove zd_up from zdiag
	478	zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1)
	479	zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp
[14156]	480	en (ji,jj,ibot) = z_en
[12377]	481	END_2D
[14834]	482	! NOTE: ctl_stop with ln_isfcav when using GLS
[12377]	483	IF( ln_isfcav) THEN ! top boundary (ocean cavity)
[14834]	484	DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
[12377]	485	itop = mikt(ji,jj) ! k top w-point
	486	itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one
	487	! ! mask at the ocean surface points
	488	z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
[2397]	489	!
	490	! Bottom level Dirichlet condition:
[14072]	491	! Bottom level (ibot) & Just above it (ibotm1)
[9019]	492	! Dirichlet ! Neumann
[12377]	493	zd_lw(ji,jj,itop) = 0._wp ! ! Remove zd_up from zdiag
	494	zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1)
	495	zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp
[14156]	496	en (ji,jj,itop) = z_en
[12377]	497	END_2D
[9019]	498	ENDIF
[2397]	499	!
[2048]	500	END SELECT
	501
	502	! Matrix inversion (en prescribed at surface and the bottom)
	503	! ----------------------------------------------------------
	504	!
[15071]	505	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
[12377]	506	zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
	507	END_3D
[15071]	508	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
[12377]	509	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)
	510	END_3D
[14834]	511	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]	512	en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
	513	END_3D
[14072]	514	! ! set the minimum value of tke
[14834]	515	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
	516	en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin )
	517	END_3D
[5109]	518
[2048]	519	!!----------------------------------------!!
	520	!! Solve prognostic equation for psi !!
	521	!!----------------------------------------!!
	522
	523	! Set psi to previous time step value
	524	!
	525	SELECT CASE ( nn_clos )
	526	!
	527	CASE( 0 ) ! k-kl (Mellor-Yamada)
[14834]	528	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
[12377]	529	psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
	530	END_3D
[2397]	531	!
[2048]	532	CASE( 1 ) ! k-eps
[14834]	533	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
[12377]	534	psi(ji,jj,jk) = eps(ji,jj,jk)
	535	END_3D
[2397]	536	!
[2048]	537	CASE( 2 ) ! k-w
[14834]	538	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
[12377]	539	psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
	540	END_3D
[2397]	541	!
	542	CASE( 3 ) ! generic
[14834]	543	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
[14072]	544	psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn
[12377]	545	END_3D
[2397]	546	!
[2048]	547	END SELECT
	548	!
	549	! Now gls (output in psi)
	550	! -------------------------------
	551	! Resolution of a tridiagonal linear system by a "methode de chasse"
	552	! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ).
[9019]	553	! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
[2048]	554	! Warning : after this step, en : right hand side of the matrix
	555
[14834]	556	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
[12377]	557	!
	558	! psi / k
[14072]	559	zratio = psi(ji,jj,jk) / eb(ji,jj,jk)
[12377]	560	!
	561	! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable)
	562	zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )
	563	!
	564	rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p
	565	!
	566	! shear prod. - stratif. destruction
	567	prod = rpsi1 * zratio * p_sh2(ji,jj,jk)
	568	!
	569	! stratif. destruction
	570	buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) )
	571	!
	572	! shear prod. - stratif. destruction
	573	diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
	574	!
	575	zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
	576	!
	577	zesh2 = zdir * ( prod + buoy ) + (1._wp - zdir ) * prod ! production term
	578	zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
[14072]	579	!
[12377]	580	! building the matrix
	581	zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
	582	! ! lower diagonal
[13237]	583	zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) &
	584	& / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) )
[12377]	585	! ! upper diagonal
[13237]	586	zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) &
	587	& / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk,Kmm) )
[12377]	588	! ! diagonal
[12489]	589	zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rn_Dt * zdiss * wmask(ji,jj,jk)
[12377]	590	! ! right hand side in psi
[12489]	591	psi(ji,jj,jk) = psi(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk)
[12377]	592	END_3D
[2048]	593	!
[15071]	594	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
	595	zdiag(ji,jj,jpk) = 1._wp
	596	END_2D
	597
[2048]	598	! Surface boundary condition on psi
	599	! ---------------------------------
	600	!
[5109]	601	SELECT CASE ( nn_bc_surf )
[2048]	602	!
	603	CASE ( 0 ) ! Dirichlet boundary conditions
[9019]	604	!
[14834]	605	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
	606	! Surface value
	607	zdep (ji,jj) = zhsro(ji,jj) * rl_sf ! Cosmetic
	608	psi (ji,jj,1) = rc0*rpp en(ji,jj,1)*rmm zdep(ji,jj)*rnn tmask(ji,jj,1)
	609	zd_lw(ji,jj,1) = psi(ji,jj,1)
	610	zd_up(ji,jj,1) = 0._wp
	611	zdiag(ji,jj,1) = 1._wp
	612	!
	613	! One level below
	614	zkar (ji,jj) = (rl_sf + (vkarmn-rl_sf)(1._wp-EXP(-rtransgdepw(ji,jj,2,Kmm)/zhsro(ji,jj) )))
	615	zdep (ji,jj) = (zhsro(ji,jj) + gdepw(ji,jj,2,Kmm)) * zkar(ji,jj)
	616	psi (ji,jj,2) = rc0*rpp en(ji,jj,2)*rmm zdep(ji,jj)*rnn tmask(ji,jj,1)
	617	zd_lw(ji,jj,2) = 0._wp
	618	zd_up(ji,jj,2) = 0._wp
	619	zdiag(ji,jj,2) = 1._wp
	620	END_2D
[9019]	621	!
[2048]	622	CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz
[9019]	623	!
[14834]	624	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
	625	! Surface value: Dirichlet
	626	zdep (ji,jj) = zhsro(ji,jj) * rl_sf
	627	psi (ji,jj,1) = rc0*rpp en(ji,jj,1)*rmm zdep(ji,jj)*rnn tmask(ji,jj,1)
	628	zd_lw(ji,jj,1) = psi(ji,jj,1)
	629	zd_up(ji,jj,1) = 0._wp
	630	zdiag(ji,jj,1) = 1._wp
	631	!
	632	! Neumann condition at k=2, zdiag zd_lw not defined/used on the halo
[13546]	633	zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
	634	zd_lw(ji,jj,2) = 0._wp
[14834]	635	!
	636	! Set psi vertical flux at the surface:
	637	zkar (ji,jj) = rl_sf + (vkarmn-rl_sf)(1._wp-EXP(-rtransgdept(ji,jj,1,Kmm)/zhsro(ji,jj) )) ! Lengh scale slope
	638	zdep (ji,jj) = ((zhsro(ji,jj) + gdept(ji,jj,1,Kmm)) / zhsro(ji,jj))*(rmmra_sf)
	639	zflxs(ji,jj) = (rnn + (1._wp-zice_fra(ji,jj))rsbc_tke1 (rnn + rmmra_sf) zdep(ji,jj)) &
	640	& (1._wp + (1._wp-zice_fra(ji,jj))rsbc_tke1zdep(ji,jj))(2._wprmm/3._wp-1_wp)
	641	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)) * &
	642	& ustar2_surf(ji,jj)*rmm zkar(ji,jj)*rnn (zhsro(ji,jj) + gdept(ji,jj,1,Kmm))**(rnn-1.)
	643	zflxs(ji,jj) = zdep(ji,jj) * zflxs(ji,jj)
	644	psi (ji,jj,2) = psi(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm)
[13546]	645	END_2D
[9019]	646	!
[2048]	647	END SELECT
	648
	649	! Bottom boundary condition on psi
	650	! --------------------------------
	651	!
[9019]	652	!!gm should be done for ISF (top boundary cond.)
	653	!!gm so, totally new staff needed ===>>> think about that !
	654	!
	655	SELECT CASE ( nn_bc_bot ) ! bottom boundary
[2048]	656	!
[14072]	657	CASE ( 0 ) ! Dirichlet
[9019]	658	! ! en(ibot) = u^2 / Co2 and hmxl_n(ibot) = vkarmn r_z0_bot
[2397]	659	! ! Balance between the production and the dissipation terms
[14834]	660	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
[12377]	661	ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point
	662	ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1
	663	zdep(ji,jj) = vkarmn * r_z0_bot
	664	psi (ji,jj,ibot) = rc0*rpp en(ji,jj,ibot)*rmm zdep(ji,jj)**rnn
	665	zd_lw(ji,jj,ibot) = 0._wp
	666	zd_up(ji,jj,ibot) = 0._wp
	667	zdiag(ji,jj,ibot) = 1._wp
	668	!
	669	! Just above last level, Dirichlet condition again (GOTM like)
	670	zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t(ji,jj,ibotm1,Kmm) )
	671	psi (ji,jj,ibotm1) = rc0*rpp en(ji,jj,ibot )*rmm zdep(ji,jj)**rnn
	672	zd_lw(ji,jj,ibotm1) = 0._wp
	673	zd_up(ji,jj,ibotm1) = 0._wp
	674	zdiag(ji,jj,ibotm1) = 1._wp
	675	END_2D
[2397]	676	!
[14967]	677	IF( ln_isfcav) THEN ! top boundary (ocean cavity)
	678	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
	679	IF ( mikt(ji,jj) > 1 ) THEN
	680	itop = mikt(ji,jj) ! k top w-point
	681	itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one
	682	!
	683	zdep(ji,jj) = vkarmn * r_z0_top
	684	psi (ji,jj,itop) = rc0*rpp en(ji,jj,itop)*rmm zdep(ji,jj)**rnn
	685	zd_lw(ji,jj,itop) = 0._wp
	686	zd_up(ji,jj,itop) = 0._wp
	687	zdiag(ji,jj,itop) = 1._wp
	688	!
	689	! Just above last level, Dirichlet condition again (GOTM like)
	690	zdep(ji,jj) = vkarmn * ( r_z0_top + e3t(ji,jj,itopp1,Kmm) )
	691	psi (ji,jj,itopp1) = rc0*rpp en(ji,jj,itop )*rmm zdep(ji,jj)**rnn
	692	zd_lw(ji,jj,itopp1) = 0._wp
	693	zd_up(ji,jj,itopp1) = 0._wp
	694	zdiag(ji,jj,itopp1) = 1._wp
	695	END IF
	696	END_2D
	697	END IF
	698	!
[2048]	699	CASE ( 1 ) ! Neumman boundary condition
[14072]	700	!
[14834]	701	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
[12377]	702	ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point
	703	ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1
	704	!
	705	! Bottom level Dirichlet condition:
	706	zdep(ji,jj) = vkarmn * r_z0_bot
	707	psi (ji,jj,ibot) = rc0*rpp en(ji,jj,ibot)*rmm zdep(ji,jj)**rnn
	708	!
	709	zd_lw(ji,jj,ibot) = 0._wp
	710	zd_up(ji,jj,ibot) = 0._wp
	711	zdiag(ji,jj,ibot) = 1._wp
	712	!
	713	! Just above last level: Neumann condition with flux injection
	714	zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag
	715	zd_up(ji,jj,ibotm1) = 0.
	716	!
	717	! Set psi vertical flux at the bottom:
	718	zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t(ji,jj,ibotm1,Kmm)
	719	zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) ) &
	720	& * (0.5_wp(en(ji,jj,ibot)+en(ji,jj,ibotm1)))rmm zdep(ji,jj)**(rnn-1._wp)
	721	psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w(ji,jj,ibotm1,Kmm)
	722	END_2D
[2397]	723	!
[14967]	724	IF( ln_isfcav) THEN ! top boundary (ocean cavity)
	725	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
	726	IF ( mikt(ji,jj) > 1 ) THEN
	727	itop = mikt(ji,jj) ! k top w-point
	728	itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one
	729	!
	730	! Bottom level Dirichlet condition:
	731	zdep(ji,jj) = vkarmn * r_z0_top
	732	psi (ji,jj,itop) = rc0*rpp en(ji,jj,itop)*rmm zdep(ji,jj)**rnn
	733	!
	734	zd_lw(ji,jj,itop) = 0._wp
	735	zd_up(ji,jj,itop) = 0._wp
	736	zdiag(ji,jj,itop) = 1._wp
	737	!
	738	! Just below cavity level: Neumann condition with flux
	739	! injection
	740	zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) ! Remove zd_up from zdiag
	741	zd_up(ji,jj,itopp1) = 0._wp
	742	!
	743	! Set psi vertical flux below cavity:
	744	zdep(ji,jj) = r_z0_top + 0.5_wp*e3t(ji,jj,itopp1,Kmm)
	745	zflxb = rsbc_psi2 * ( p_avm(ji,jj,itop) + p_avm(ji,jj,itopp1)) &
	746	& * (0.5_wp(en(ji,jj,itop)+en(ji,jj,itopp1)))rmm zdep(ji,jj)**(rnn-1._wp)
	747	psi(ji,jj,itopp1) = psi(ji,jj,itopp1) + zflxb / e3w(ji,jj,itopp1,Kmm)
	748	END IF
	749	END_2D
	750	END IF
	751
	752	!
[2048]	753	END SELECT
	754
	755	! Matrix inversion
	756	! ----------------
	757	!
[14834]	758	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]	759	zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
	760	END_3D
[14834]	761	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]	762	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)
	763	END_3D
[14834]	764	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]	765	psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
	766	END_3D
[2048]	767
	768	! Set dissipation
	769	!----------------
	770
	771	SELECT CASE ( nn_clos )
	772	!
	773	CASE( 0 ) ! k-kl (Mellor-Yamada)
[14834]	774	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
[12377]	775	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)
	776	END_3D
[2397]	777	!
[2048]	778	CASE( 1 ) ! k-eps
[14834]	779	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
[12377]	780	eps(ji,jj,jk) = psi(ji,jj,jk)
	781	END_3D
[2397]	782	!
[2048]	783	CASE( 2 ) ! k-w
[14834]	784	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
[14072]	785	eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)
[12377]	786	END_3D
[2397]	787	!
	788	CASE( 3 ) ! generic
	789	zcoef = rc0**( 3._wp + rpp/rnn )
	790	zex1 = ( 1.5_wp + rmm/rnn )
	791	zex2 = -1._wp / rnn
[14834]	792	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
[12377]	793	eps(ji,jj,jk) = zcoef * en(ji,jj,jk)*zex1 psi(ji,jj,jk)**zex2
	794	END_3D
[2397]	795	!
[2048]	796	END SELECT
	797
	798	! Limit dissipation rate under stable stratification
	799	! --------------------------------------------------
[14834]	800	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]	801	! limitation
	802	eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin )
	803	hmxl_n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
	804	END_3D
[14834]	805	IF( ln_length_lim ) THEN ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
	806	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
	807	zrn2 = MAX( rn2(ji,jj,jk), rsmall )
	808	hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) )
	809	END_3D
	810	ENDIF
[2048]	811
	812	!
	813	! Stability function and vertical viscosity and diffusivity
	814	! ---------------------------------------------------------
	815	!
	816	SELECT CASE ( nn_stab_func )
	817	!
	818	CASE ( 0 , 1 ) ! Galperin or Kantha-Clayson stability functions
[14834]	819	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
[12377]	820	! zcof = l²/q²
	821	zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
	822	! Gh = -N²l²/q²
	823	gh = - rn2(ji,jj,jk) * zcof
	824	gh = MIN( gh, rgh0 )
	825	gh = MAX( gh, rghmin )
	826	! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
	827	sh = ra2( 1._wp-6._wpra1/rb1 ) / ( 1.-3.ra2gh(6.ra1+rb2*( 1._wp-rc3 ) ) )
	828	sm = ( rb1*(-1._wp/3._wp) + ( 18._wpra1ra1 + 9._wpra1ra2(1._wp-rc2) )shgh ) / (1._wp-9._wpra1ra2*gh)
	829	!
	830	! Store stability function in zstt and zstm
	831	zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
	832	zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
	833	END_3D
[2397]	834	!
[2048]	835	CASE ( 2, 3 ) ! Canuto stability functions
[14834]	836	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
[12377]	837	! zcof = l²/q²
	838	zcof = hmxl_b(ji,jj,jk)hmxl_b(ji,jj,jk) / ( 2._wp eb(ji,jj,jk) )
	839	! Gh = -N²l²/q²
	840	gh = - rn2(ji,jj,jk) * zcof
	841	gh = MIN( gh, rgh0 )
	842	gh = MAX( gh, rghmin )
	843	gh = gh * rf6
	844	! Gm = M²l²/q² Shear number
	845	shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall )
	846	gm = MAX( shr * zcof , 1.e-10 )
	847	gm = gm * rf6
	848	gm = MIN ( (rd0 - rd1gh + rd3ghgh) / (rd2-rd4gh) , gm )
	849	! Stability functions from Canuto
	850	rcff = rd0 - rd1gh +rd2gm + rd3ghgh - rd4ghgm + rd5gmgm
	851	sm = (rs0 - rs1gh + rs2gm) / rcff
	852	sh = (rs4 - rs5gh + rs6gm) / rcff
	853	!
	854	! Store stability function in zstt and zstm
	855	zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
	856	zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
	857	END_3D
[2397]	858	!
[2048]	859	END SELECT
	860
	861	! Boundary conditions on stability functions for momentum (Neumann):
	862	! Lines below are useless if GOTM style Dirichlet conditions are used
[5109]	863
[15071]	864	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
	865	zstm(ji,jj,1) = zstm(ji,jj,2)
	866	zstm(ji,jj,jpk) = 0. ! default value, in case jpk > mbkt(ji,jj)+1
	867	! ! Not needed but avoid a bug when looking for undefined values (-fpe0)
	868	END_2D
[14834]	869	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! update bottom with good values
[12377]	870	zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
	871	END_2D
[10342]	872
[14834]	873	zstt(:,:, 1) = wmask(A2D(nn_hls), 1) ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
	874	zstt(:,:,jpk) = wmask(A2D(nn_hls),jpk) ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
[10342]	875
[9019]	876	!!gm should be done for ISF (top boundary cond.)
	877	!!gm so, totally new staff needed!!gm
[2048]	878
	879	! Compute diffusivities/viscosities
	880	! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
[10342]	881	! -> yes BUT p_avm(:,:1) and p_avm(:,:jpk) are used when we compute zd_lw(:,:2) and zd_up(:,:jpkm1). These values are
	882	! later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk)
	883	! for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0)
[14834]	884	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
[12377]	885	zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
	886	zavt = zsqen * zstt(ji,jj,jk)
	887	zavm = zsqen * zstm(ji,jj,jk)
	888	p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
	889	p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom
	890	END_3D
[14834]	891	p_avt(A2D(nn_hls),1) = 0._wp
[2048]	892	!
[12377]	893	IF(sn_cfctl%l_prtctl) THEN
[9440]	894	CALL prt_ctl( tab3d_1=en , clinfo1=' gls - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk)
	895	CALL prt_ctl( tab3d_1=p_avm, clinfo1=' gls - m: ', kdim=jpk )
[2048]	896	ENDIF
	897	!
	898	END SUBROUTINE zdf_gls
	899
[2329]	900
[2048]	901	SUBROUTINE zdf_gls_init
	902	!!----------------------------------------------------------------------
	903	!! * ROUTINE zdf_gls_init *
[14072]	904	!!
	905	!! ** Purpose : Initialization of the vertical eddy diffivity and
[9019]	906	!! viscosity computed using a GLS turbulent closure scheme
[2048]	907	!!
	908	!! ** Method : Read the namzdf_gls namelist and check the parameters
	909	!!
	910	!! ** input : Namlist namzdf_gls
	911	!!
	912	!! ** Action : Increase by 1 the nstop flag is setting problem encounter
	913	!!
	914	!!----------------------------------------------------------------------
[2329]	915	INTEGER :: jk ! dummy loop indices
[4147]	916	INTEGER :: ios ! Local integer output status for namelist read
[2329]	917	REAL(wp):: zcr ! local scalar
[2048]	918	!!
[13472]	919	NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
	920	& rn_clim_galp, ln_sigpsi, rn_hsro, rn_hsri, &
[14966]	921	& nn_mxlice, rn_crban, rn_charn, rn_frac_hs, &
[13472]	922	& nn_bc_surf, nn_bc_bot, nn_z0_met, nn_z0_ice, &
[2048]	923	& nn_stab_func, nn_clos
	924	!!----------------------------------------------------------
[3294]	925	!
[4147]	926	READ ( numnam_ref, namzdf_gls, IOSTAT = ios, ERR = 901)
[11536]	927	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in reference namelist' )
[2048]	928
[4147]	929	READ ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 )
[11536]	930	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist' )
[4624]	931	IF(lwm) WRITE ( numond, namzdf_gls )
[4147]	932
[2397]	933	IF(lwp) THEN !* Control print
[2048]	934	WRITE(numout,*)
[9019]	935	WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme'
[2048]	936	WRITE(numout,*) '~~~~~~~~~~~~'
[2397]	937	WRITE(numout,*) ' Namelist namzdf_gls : set gls mixing parameters'
[5109]	938	WRITE(numout,*) ' minimum value of en rn_emin = ', rn_emin
	939	WRITE(numout,*) ' minimum value of eps rn_epsmin = ', rn_epsmin
	940	WRITE(numout,*) ' Limit dissipation rate under stable stratif. ln_length_lim = ', ln_length_lim
	941	WRITE(numout,*) ' Galperin limit (Standard: 0.53, Holt: 0.26) rn_clim_galp = ', rn_clim_galp
	942	WRITE(numout,*) ' TKE Surface boundary condition nn_bc_surf = ', nn_bc_surf
	943	WRITE(numout,*) ' TKE Bottom boundary condition nn_bc_bot = ', nn_bc_bot
	944	WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi
[2397]	945	WRITE(numout,*) ' Craig and Banner coefficient rn_crban = ', rn_crban
	946	WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn
[5109]	947	WRITE(numout,*) ' Surface roughness formula nn_z0_met = ', nn_z0_met
[13472]	948	WRITE(numout,*) ' surface wave breaking under ice nn_z0_ice = ', nn_z0_ice
	949	SELECT CASE( nn_z0_ice )
	950	CASE( 0 ) ; WRITE(numout,*) ' ==>>> no impact of ice cover on surface wave breaking'
	951	CASE( 1 ) ; WRITE(numout,) ' ==>>> roughness uses rn_hsri and is weigthed by 1-TANH( fr_i(:,:) 10 )'
	952	CASE( 2 ) ; WRITE(numout,*) ' ==>>> roughness uses rn_hsri and is weighted by 1-fr_i(:,:)'
	953	CASE( 3 ) ; WRITE(numout,) ' ==>>> roughness uses rn_hsri and is weighted by 1-MIN( 1, 4 fr_i(:,:) )'
	954	CASE DEFAULT
	955	CALL ctl_stop( 'zdf_gls_init: wrong value for nn_z0_ice, should be 0,1,2, or 3')
	956	END SELECT
[5109]	957	WRITE(numout,*) ' Wave height frac. (used if nn_z0_met=2) rn_frac_hs = ', rn_frac_hs
[2397]	958	WRITE(numout,*) ' Stability functions nn_stab_func = ', nn_stab_func
	959	WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos
[5109]	960	WRITE(numout,*) ' Surface roughness (m) rn_hsro = ', rn_hsro
[14966]	961	WRITE(numout,*) ' type of scaling under sea-ice nn_mxlice = ', nn_mxlice
	962	IF( nn_mxlice == 1 ) &
	963	WRITE(numout,*) ' Ice-ocean roughness (used if nn_z0_ice/=0) rn_hsri = ', rn_hsri
	964	SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice
	965	CASE( 0 ) ; WRITE(numout,*) ' ==>>> No scaling under sea-ice'
	966	CASE( 1 ) ; WRITE(numout,*) ' ==>>> scaling with constant sea-ice thickness'
	967	CASE( 2 ) ; WRITE(numout,*) ' ==>>> scaling with mean sea-ice thickness'
	968	CASE( 3 ) ; WRITE(numout,*) ' ==>>> scaling with max sea-ice thickness'
	969	CASE DEFAULT
	970	CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 ')
	971	END SELECT
[9019]	972	WRITE(numout,*)
[2048]	973	ENDIF
	974
[9019]	975	! !* allocate GLS arrays
[2715]	976	IF( zdf_gls_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
	977
[2397]	978	! !* Check of some namelist values
[13283]	979	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' )
	980	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' )
	981	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' )
	982	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' )
	983	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' )
	984	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]	985
[2715]	986	SELECT CASE ( nn_clos ) !* set the parameters for the chosen closure
[2048]	987	!
[2715]	988	CASE( 0 ) ! k-kl (Mellor-Yamada)
[2397]	989	!
[9019]	990	IF(lwp) WRITE(numout,*) ' ==>> k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)'
	991	IF(lwp) WRITE(numout,*)
[2397]	992	rpp = 0._wp
	993	rmm = 1._wp
	994	rnn = 1._wp
	995	rsc_tke = 1.96_wp
	996	rsc_psi = 1.96_wp
	997	rpsi1 = 0.9_wp
	998	rpsi3p = 1._wp
	999	rpsi2 = 0.5_wp
	1000	!
[2048]	1001	SELECT CASE ( nn_stab_func )
[2397]	1002	CASE( 0, 1 ) ; rpsi3m = 2.53_wp ! G88 or KC stability functions
[5109]	1003	CASE( 2 ) ; rpsi3m = 2.62_wp ! Canuto A stability functions
[2397]	1004	CASE( 3 ) ; rpsi3m = 2.38 ! Canuto B stability functions (caution : constant not identified)
	1005	END SELECT
[2048]	1006	!
[2715]	1007	CASE( 1 ) ! k-eps
[2397]	1008	!
[9019]	1009	IF(lwp) WRITE(numout,*) ' ==>> k-eps closure chosen'
	1010	IF(lwp) WRITE(numout,*)
[2397]	1011	rpp = 3._wp
	1012	rmm = 1.5_wp
	1013	rnn = -1._wp
	1014	rsc_tke = 1._wp
[5109]	1015	rsc_psi = 1.2_wp ! Schmidt number for psi
[2397]	1016	rpsi1 = 1.44_wp
	1017	rpsi3p = 1._wp
	1018	rpsi2 = 1.92_wp
	1019	!
	1020	SELECT CASE ( nn_stab_func )
	1021	CASE( 0, 1 ) ; rpsi3m = -0.52_wp ! G88 or KC stability functions
	1022	CASE( 2 ) ; rpsi3m = -0.629_wp ! Canuto A stability functions
	1023	CASE( 3 ) ; rpsi3m = -0.566 ! Canuto B stability functions
[2048]	1024	END SELECT
[2397]	1025	!
[2715]	1026	CASE( 2 ) ! k-omega
[2397]	1027	!
[9019]	1028	IF(lwp) WRITE(numout,*) ' ==>> k-omega closure chosen'
	1029	IF(lwp) WRITE(numout,*)
[2397]	1030	rpp = -1._wp
	1031	rmm = 0.5_wp
	1032	rnn = -1._wp
	1033	rsc_tke = 2._wp
	1034	rsc_psi = 2._wp
	1035	rpsi1 = 0.555_wp
	1036	rpsi3p = 1._wp
	1037	rpsi2 = 0.833_wp
	1038	!
	1039	SELECT CASE ( nn_stab_func )
	1040	CASE( 0, 1 ) ; rpsi3m = -0.58_wp ! G88 or KC stability functions
	1041	CASE( 2 ) ; rpsi3m = -0.64_wp ! Canuto A stability functions
	1042	CASE( 3 ) ; rpsi3m = -0.64_wp ! Canuto B stability functions caution : constant not identified)
	1043	END SELECT
	1044	!
[2715]	1045	CASE( 3 ) ! generic
[2397]	1046	!
[9019]	1047	IF(lwp) WRITE(numout,*) ' ==>> generic closure chosen'
	1048	IF(lwp) WRITE(numout,*)
[2397]	1049	rpp = 2._wp
	1050	rmm = 1._wp
	1051	rnn = -0.67_wp
	1052	rsc_tke = 0.8_wp
	1053	rsc_psi = 1.07_wp
	1054	rpsi1 = 1._wp
	1055	rpsi3p = 1._wp
	1056	rpsi2 = 1.22_wp
	1057	!
	1058	SELECT CASE ( nn_stab_func )
	1059	CASE( 0, 1 ) ; rpsi3m = 0.1_wp ! G88 or KC stability functions
	1060	CASE( 2 ) ; rpsi3m = 0.05_wp ! Canuto A stability functions
	1061	CASE( 3 ) ; rpsi3m = 0.05_wp ! Canuto B stability functions caution : constant not identified)
	1062	END SELECT
	1063	!
[2048]	1064	END SELECT
	1065
	1066	!
[2715]	1067	SELECT CASE ( nn_stab_func ) !* set the parameters of the stability functions
[2048]	1068	!
[2715]	1069	CASE ( 0 ) ! Galperin stability functions
[2397]	1070	!
[9019]	1071	IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Galperin'
[2397]	1072	rc2 = 0._wp
	1073	rc3 = 0._wp
	1074	rc_diff = 1._wp
	1075	rc0 = 0.5544_wp
	1076	rcm_sf = 0.9884_wp
	1077	rghmin = -0.28_wp
	1078	rgh0 = 0.0233_wp
	1079	rghcri = 0.02_wp
	1080	!
[2715]	1081	CASE ( 1 ) ! Kantha-Clayson stability functions
[2397]	1082	!
[9019]	1083	IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Kantha-Clayson'
[2397]	1084	rc2 = 0.7_wp
	1085	rc3 = 0.2_wp
	1086	rc_diff = 1._wp
	1087	rc0 = 0.5544_wp
	1088	rcm_sf = 0.9884_wp
	1089	rghmin = -0.28_wp
	1090	rgh0 = 0.0233_wp
	1091	rghcri = 0.02_wp
	1092	!
[2715]	1093	CASE ( 2 ) ! Canuto A stability functions
[2397]	1094	!
[9019]	1095	IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Canuto A'
[2397]	1096	rs0 = 1.5_wp * rl1 * rl5*rl5
	1097	rs1 = -rl4(rl6+rl7) + 2._wprl4rl5(rl1-(1._wp/3._wp)rl2-rl3) + 1.5_wprl1rl5rl8
	1098	rs2 = -(3._wp/8._wp) * rl1(rl6rl6-rl7*rl7)
	1099	rs4 = 2._wp * rl5
	1100	rs5 = 2._wp * rl4
	1101	rs6 = (2._wp/3._wp) * rl5 * ( 3._wprl3rl3 - rl2rl2 ) - 0.5_wp rl5rl1 (3._wp*rl3-rl2) &
	1102	& + 0.75_wp * rl1 * ( rl6 - rl7 )
	1103	rd0 = 3._wp * rl5*rl5
	1104	rd1 = rl5 * ( 7._wprl4 + 3._wprl8 )
	1105	rd2 = rl5rl5 ( 3._wprl3rl3 - rl2rl2 ) - 0.75_wp(rl6rl6 - rl7rl7 )
	1106	rd3 = rl4 * ( 4._wprl4 + 3._wprl8)
	1107	rd4 = rl4 * ( rl2 * rl6 - 3._wprl3rl7 - rl5(rl2rl2 - rl3rl3 ) ) + rl5rl8 * ( 3._wprl3rl3 - rl2*rl2 )
	1108	rd5 = 0.25_wp * ( rl2rl2 - 3._wp rl3rl3 ) ( rl6rl6 - rl7rl7 )
	1109	rc0 = 0.5268_wp
	1110	rf6 = 8._wp / (rc0**6._wp)
	1111	rc_diff = SQRT(2._wp) / (rc0**3._wp)
	1112	rcm_sf = 0.7310_wp
	1113	rghmin = -0.28_wp
	1114	rgh0 = 0.0329_wp
	1115	rghcri = 0.03_wp
	1116	!
[2715]	1117	CASE ( 3 ) ! Canuto B stability functions
[2397]	1118	!
[9019]	1119	IF(lwp) WRITE(numout,*) ' ==>> Stability functions from Canuto B'
[2397]	1120	rs0 = 1.5_wp * rm1 * rm5*rm5
	1121	rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4rm5(rm1-(1._wp/3._wp)rm2-rm3) + 1.5_wp rm1rm5rm8
	1122	rs2 = -(3._wp/8._wp) * rm1 * (rm6rm6-rm7rm7 )
	1123	rs4 = 2._wp * rm5
	1124	rs5 = 2._wp * rm4
	1125	rs6 = (2._wp/3._wp) * rm5 * (3._wprm3rm3-rm2rm2) - 0.5_wp rm5rm1(3._wprm3-rm2) + 0.75_wp rm1*(rm6-rm7)
	1126	rd0 = 3._wp * rm5*rm5
	1127	rd1 = rm5 * (7._wprm4 + 3._wprm8)
	1128	rd2 = rm5rm5 (3._wprm3rm3 - rm2rm2) - 0.75_wp (rm6rm6 - rm7rm7)
	1129	rd3 = rm4 * ( 4._wprm4 + 3._wprm8 )
	1130	rd4 = rm4 * ( rm2rm6 -3._wprm3rm7 - rm5(rm2rm2 - rm3rm3) ) + rm5 * rm8 * ( 3._wprm3rm3 - rm2*rm2 )
	1131	rd5 = 0.25_wp * ( rm2rm2 - 3._wprm3rm3 ) ( rm6rm6 - rm7rm7 )
	1132	rc0 = 0.5268_wp !! rc0 = 0.5540_wp (Warner ...) to verify !
	1133	rf6 = 8._wp / ( rc0**6._wp )
	1134	rc_diff = SQRT(2._wp)/(rc0**3.)
	1135	rcm_sf = 0.7470_wp
	1136	rghmin = -0.28_wp
	1137	rgh0 = 0.0444_wp
	1138	rghcri = 0.0414_wp
	1139	!
[2048]	1140	END SELECT
[14072]	1141
[2715]	1142	! !* Set Schmidt number for psi diffusion in the wave breaking case
	1143	! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
	1144	! ! or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
[5109]	1145	IF( ln_sigpsi ) THEN
[14072]	1146	ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf
[5109]	1147	! Verification: retrieve Burchard (2001) results by uncomenting the line below:
	1148	! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work
	1149	! ra_sf = -SQRT(2./3.rc03./rn_cm_sfrn_sc_tke)/vkarmn
	1150	rsc_psi0 = rsc_tke/(24.rpsi2)(-1.+(4.rnn + ra_sf(1.+4.rmm))2./(ra_sf*2.))
[2048]	1151	ELSE
[2299]	1152	rsc_psi0 = rsc_psi
[2048]	1153	ENDIF
[14072]	1154
[2715]	1155	! !* Shear free turbulence parameters
[2048]	1156	!
[5109]	1157	ra_sf = -4._wprnnSQRT(rsc_tke) / ( (1._wp+4._wprmm)SQRT(rsc_tke) &
	1158	& - SQRT(rsc_tke + 24._wprsc_psi0rpsi2 ) )
[2048]	1159
[5109]	1160	IF ( rn_crban==0._wp ) THEN
	1161	rl_sf = vkarmn
	1162	ELSE
[9019]	1163	rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wprmm + 8._wprmm*2_wp) rsc_tke &
	1164	& + 12._wprsc_psi0rpsi2 - (1._wp + 4._wp*rmm) &
	1165	& SQRT(rsc_tke(rsc_tke &
	1166	& + 24._wprsc_psi0rpsi2)) ) &
	1167	& /(12._wprnn*2.) )
[5109]	1168	ENDIF
	1169
[2048]	1170	!
[2715]	1171	IF(lwp) THEN !* Control print
[2048]	1172	WRITE(numout,*)
[9019]	1173	WRITE(numout,*) ' Limit values :'
	1174	WRITE(numout,*) ' Parameter m = ', rmm
	1175	WRITE(numout,*) ' Parameter n = ', rnn
	1176	WRITE(numout,*) ' Parameter p = ', rpp
	1177	WRITE(numout,*) ' rpsi1 = ', rpsi1
	1178	WRITE(numout,*) ' rpsi2 = ', rpsi2
	1179	WRITE(numout,*) ' rpsi3m = ', rpsi3m
	1180	WRITE(numout,*) ' rpsi3p = ', rpsi3p
	1181	WRITE(numout,*) ' rsc_tke = ', rsc_tke
	1182	WRITE(numout,*) ' rsc_psi = ', rsc_psi
	1183	WRITE(numout,*) ' rsc_psi0 = ', rsc_psi0
	1184	WRITE(numout,*) ' rc0 = ', rc0
[2048]	1185	WRITE(numout,*)
[9019]	1186	WRITE(numout,*) ' Shear free turbulence parameters:'
	1187	WRITE(numout,*) ' rcm_sf = ', rcm_sf
	1188	WRITE(numout,*) ' ra_sf = ', ra_sf
	1189	WRITE(numout,*) ' rl_sf = ', rl_sf
[2048]	1190	ENDIF
	1191
[2715]	1192	! !* Constants initialization
[2397]	1193	rc02 = rc0 * rc0 ; rc02r = 1. / rc02
	1194	rc03 = rc02 * rc0
	1195	rc04 = rc03 * rc0
[5109]	1196	rsbc_tke1 = -3._wp/2._wprn_crbanra_sf*rl_sf ! Dirichlet + Wave breaking
[14072]	1197	rsbc_tke2 = rn_Dt * rn_crban / rl_sf ! Neumann + Wave breaking
[5109]	1198	zcr = MAX(rsmall, rsbc_tke1*(1./(-ra_sf3._wp/2._wp))-1._wp )
[14072]	1199	rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer
[5109]	1200	rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness
[14072]	1201	rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness
[12489]	1202	rsbc_psi1 = -0.5_wp * rn_Dt * rc0*(rpp-2._wprmm) / rsc_psi
[14072]	1203	rsbc_psi2 = -0.5_wp * rn_Dt * rc0*rpp rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking
[9019]	1204	!
[12489]	1205	rfact_tke = -0.5_wp / rsc_tke * rn_Dt ! Cst used for the Diffusion term of tke
	1206	rfact_psi = -0.5_wp / rsc_psi * rn_Dt ! Cst used for the Diffusion term of tke
[9019]	1207	!
[2397]	1208	! !* Wall proximity function
[9019]	1209	!!gm tmask or wmask ????
	1210	zwall(:,:,:) = 1._wp * tmask(:,:,:)
[2048]	1211
[14072]	1212	! !* read or initialize all required files
[9019]	1213	CALL gls_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, hmxl_n)
[2048]	1214	!
	1215	END SUBROUTINE zdf_gls_init
	1216
[2329]	1217
[2048]	1218	SUBROUTINE gls_rst( kt, cdrw )
[2452]	1219	!!---------------------------------------------------------------------
[9124]	1220	!! * ROUTINE gls_rst *
[14072]	1221	!!
[2452]	1222	!! ** Purpose : Read or write TKE file (en) in restart file
	1223	!!
	1224	!! ** Method : use of IOM library
[14072]	1225	!! if the restart does not contain TKE, en is either
[2452]	1226	!! set to rn_emin or recomputed (nn_igls/=0)
	1227	!!----------------------------------------------------------------------
[9019]	1228	USE zdf_oce , ONLY : en, avt_k, avm_k ! ocean vertical physics
	1229	!!
	1230	INTEGER , INTENT(in) :: kt ! ocean time-step
	1231	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
[2452]	1232	!
	1233	INTEGER :: jit, jk ! dummy loop indices
[9019]	1234	INTEGER :: id1, id2, id3, id4
[2452]	1235	INTEGER :: ji, jj, ikbu, ikbv
	1236	REAL(wp):: cbx, cby
	1237	!!----------------------------------------------------------------------
	1238	!
[14072]	1239	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
[2452]	1240	! ! ---------------
	1241	IF( ln_rstart ) THEN !* Read the restart file
[9019]	1242	id1 = iom_varid( numror, 'en' , ldstop = .FALSE. )
	1243	id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. )
	1244	id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. )
	1245	id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. )
[2452]	1246	!
[9019]	1247	IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! all required arrays exist
[15145]	1248	CALL iom_get( numror, jpdom_auto, 'en' , en , kfill = jpfillcopy ) ! we devide by en -> must be != 0.
[13970]	1249	CALL iom_get( numror, jpdom_auto, 'avt_k' , avt_k )
	1250	CALL iom_get( numror, jpdom_auto, 'avm_k' , avm_k )
[15145]	1251	CALL iom_get( numror, jpdom_auto, 'hmxl_n', hmxl_n, kfill = jpfillcopy ) ! we devide by hmxl_n -> must be != 0.
[14072]	1252	ELSE
[9019]	1253	IF(lwp) WRITE(numout,*)
	1254	IF(lwp) WRITE(numout,*) ' ==>> previous run without GLS scheme, set en and hmxl_n to background values'
	1255	en (:,:,:) = rn_emin
	1256	hmxl_n(:,:,:) = 0.05_wp
	1257	! avt_k, avm_k already set to the background value in zdf_phy_init
[2452]	1258	ENDIF
	1259	ELSE !* Start from rest
[9019]	1260	IF(lwp) WRITE(numout,*)
	1261	IF(lwp) WRITE(numout,*) ' ==>> start from rest, set en and hmxl_n by background values'
	1262	en (:,:,:) = rn_emin
	1263	hmxl_n(:,:,:) = 0.05_wp
	1264	! avt_k, avm_k already set to the background value in zdf_phy_init
[2452]	1265	ENDIF
	1266	!
	1267	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
	1268	! ! -------------------
	1269	IF(lwp) WRITE(numout,*) '---- gls-rst ----'
[13970]	1270	CALL iom_rstput( kt, nitrst, numrow, 'en' , en )
	1271	CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k )
	1272	CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k )
	1273	CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n )
[2452]	1274	!
	1275	ENDIF
	1276	!
[2048]	1277	END SUBROUTINE gls_rst
	1278
	1279	!!======================================================================
	1280	END MODULE zdfgls

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

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

Context Navigation

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

Download in other formats: