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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

zdfgls.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfgls.F90 @ 13108

Last change on this file since 13108 was 11822, checked in by acc, 5 years ago
Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Sette tested updates to branch to align with trunk changes between 10721 and 11740. Sette tests are passing but results differ from branch before these changes (except for GYRE_PISCES and VORTEX) and branch results already differed from trunk because of algorithmic fixes. Will need more checks to confirm correctness.
Property svn:keywords set to `Id`
File size: 59.3 KB

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

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

Download in other formats: