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 branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 @ 7953

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

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

Download in other formats: