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.

domvvl.F90 in NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM/domvvl.F90 @ 10164

Last change on this file since 10164 was 10164, checked in by jchanut, 6 years ago
ztilde update, #2126: Add higher order filtering + alternative scale factor decomposition
Property svn:keywords set to `Id`
File size: 120.8 KB

Rev	Line
[592]	1	MODULE domvvl
	2	!!======================================================================
	3	!! * MODULE domvvl *
	4	!! Ocean :
	5	!!======================================================================
[1438]	6	!! History : 2.0 ! 2006-06 (B. Levier, L. Marie) original code
	7	!! 3.1 ! 2009-02 (G. Madec, M. Leclair, R. Benshila) pure z* coordinate
[9019]	8	!! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates
[5120]	9	!! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability
[10116]	10	!! 4.0 ! 2018-09 (J. Chanut) improve z_tilde robustness
[592]	11	!!----------------------------------------------------------------------
[5836]	12
[592]	13	!!----------------------------------------------------------------------
[4292]	14	!! dom_vvl_init : define initial vertical scale factors, depths and column thickness
	15	!! dom_vvl_sf_nxt : Compute next vertical scale factors
	16	!! dom_vvl_sf_swp : Swap vertical scale factors and update the vertical grid
	17	!! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
	18	!! dom_vvl_rst : read/write restart file
	19	!! dom_vvl_ctl : Check the vvl options
	20	!!----------------------------------------------------------------------
[592]	21	USE oce ! ocean dynamics and tracers
[6140]	22	USE phycst ! physical constant
[592]	23	USE dom_oce ! ocean space and time domain
[4292]	24	USE sbc_oce ! ocean surface boundary condition
[6152]	25	USE wet_dry ! wetting and drying
[7646]	26	USE usrdef_istate ! user defined initial state (wad only)
[6140]	27	USE restart ! ocean restart
	28	!
[592]	29	USE in_out_manager ! I/O manager
[4292]	30	USE iom ! I/O manager library
[592]	31	USE lib_mpp ! distributed memory computing library
	32	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[3294]	33	USE timing ! Timing
[10116]	34	USE bdy_oce ! ocean open boundary conditions
	35	USE sbcrnf ! river runoff
	36	USE dynspg_ts, ONLY: un_adv, vn_adv
[592]	37
	38	IMPLICIT NONE
	39	PRIVATE
	40
[4292]	41	PUBLIC dom_vvl_init ! called by domain.F90
	42	PUBLIC dom_vvl_sf_nxt ! called by step.F90
	43	PUBLIC dom_vvl_sf_swp ! called by step.F90
	44	PUBLIC dom_vvl_interpol ! called by dynnxt.F90
[592]	45
[5836]	46	! !!* Namelist nam_vvl
	47	LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate
	48	LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate
	49	LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate
[10116]	50	LOGICAL :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate
	51	LOGICAL :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate
	52	LOGICAL :: ln_vvl_zstar_on_shelf = .FALSE. ! revert to zstar on shelves
	53	LOGICAL :: ln_vvl_adv_fct = .FALSE. ! Centred thickness advection
	54	LOGICAL :: ln_vvl_adv_cn2 = .TRUE. ! FCT thickness advection
	55	LOGICAL :: ln_vvl_dbg = .FALSE. ! debug control prints
	56	LOGICAL :: ln_vvl_ramp = .FALSE. ! Ramp on interfaces displacement
	57	LOGICAL :: ln_vvl_lap = .FALSE. ! Laplacian thickness diffusion
	58	LOGICAL :: ln_vvl_blp = .FALSE. ! Bilaplacian thickness diffusion
	59	LOGICAL :: ln_vvl_regrid = .FALSE. ! ensure layer separation
	60	LOGICAL :: ll_shorizd = .FALSE. ! Use "shelf horizon depths"
	61	LOGICAL :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer
[5836]	62	! ! conservation: not used yet
[10164]	63	INTEGER :: nn_filt_order=1
[10116]	64	REAL(wp) :: rn_ahe3_lap ! thickness diffusion coefficient (Laplacian)
	65	REAL(wp) :: rn_ahe3_blp ! thickness diffusion coefficient (Bilaplacian)
	66	REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days]
	67	REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days]
	68	REAL(wp) :: rn_day_ramp ! Duration of linear ramp [days]
	69	REAL(wp) :: hsmall=0.01_wp ! small thickness [m]
[592]	70
[5836]	71	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport
[10164]	72	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: un_lf, vn_lf, hdivn_lf ! low frequency fluxes and divergence
[5836]	73	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors
[10164]	74	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a ! baroclinic scale factors
[10116]	75	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: tildemask ! mask tilde tendency
	76	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! restoring period for scale factors
	77	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! restoring period for low freq. divergence
	78	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: hsm, dsm !
	79	INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: i_int_bot
[1438]	80
[592]	81	!! * Substitutions
	82	# include "vectopt_loop_substitute.h90"
	83	!!----------------------------------------------------------------------
[10126]	84	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]	85	!! $Id$
[10126]	86	!! Software governed by the CeCILL license (see ./LICENSE)
[592]	87	!!----------------------------------------------------------------------
[4292]	88	CONTAINS
	89
[2715]	90	INTEGER FUNCTION dom_vvl_alloc()
	91	!!----------------------------------------------------------------------
[4292]	92	!! * FUNCTION dom_vvl_alloc *
[2715]	93	!!----------------------------------------------------------------------
[6140]	94	IF( ln_vvl_zstar ) dom_vvl_alloc = 0
[4292]	95	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
[4338]	96	ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , &
[10164]	97	& un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , &
[10116]	98	& tildemask(jpi,jpj) , hsm(jpi,jpj) , dsm(jpi,jpj) , i_int_bot(jpi,jpj), STAT = dom_vvl_alloc )
[4292]	99	IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc )
	100	IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
[6140]	101	un_td = 0._wp
	102	vn_td = 0._wp
[4292]	103	ENDIF
	104	IF( ln_vvl_ztilde ) THEN
[10164]	105	ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdivn_lf(jpi,jpj,jpk,nn_filt_order), &
	106	& un_lf(jpi,jpj,jpk,nn_filt_order), vn_lf(jpi,jpj,jpk,nn_filt_order), STAT= dom_vvl_alloc )
[4292]	107	IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc )
	108	IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
	109	ENDIF
[6140]	110	!
[2715]	111	END FUNCTION dom_vvl_alloc
	112
	113
[4292]	114	SUBROUTINE dom_vvl_init
[592]	115	!!----------------------------------------------------------------------
[4292]	116	!! * ROUTINE dom_vvl_init *
[592]	117	!!
[4292]	118	!! ** Purpose : Initialization of all scale factors, depths
	119	!! and water column heights
	120	!!
	121	!! ** Method : - use restart file and/or initialize
	122	!! - interpolate scale factors
	123	!!
[6140]	124	!! ** Action : - e3t_(n/b) and tilde_e3t_(n/b)
	125	!! - Regrid: e3(u/v)_n
	126	!! e3(u/v)_b
	127	!! e3w_n
	128	!! e3(u/v)w_b
	129	!! e3(u/v)w_n
	130	!! gdept_n, gdepw_n and gde3w_n
[4292]	131	!! - h(t/u/v)_0
	132	!! - frq_rst_e3t and frq_rst_hdv
	133	!!
	134	!! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling.
[592]	135	!!----------------------------------------------------------------------
[6140]	136	INTEGER :: ji, jj, jk
[4292]	137	INTEGER :: ii0, ii1, ij0, ij1
[10116]	138	REAL(wp):: zcoef, zwgt, ztmp, zhmin, zhmax
[592]	139	!!----------------------------------------------------------------------
[6140]	140	!
[4292]	141	IF(lwp) WRITE(numout,*)
	142	IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
	143	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
[6140]	144	!
	145	CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer)
	146	!
	147	! ! Allocate module arrays
[4292]	148	IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
[6140]	149	!
[10116]	150	! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdivn_lf
[4292]	151	CALL dom_vvl_rst( nit000, 'READ' )
[7753]	152	e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all
[6140]	153	!
	154	! !== Set of all other vertical scale factors ==! (now and before)
	155	! ! Horizontal interpolation of e3t
	156	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U
	157	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
	158	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V
	159	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
[10060]	160	CALL dom_vvl_interpol( e3t_n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F
[6140]	161	! ! Vertical interpolation of e3t,u,v
	162	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W
	163	CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W' )
	164	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW
	165	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
	166	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW
	167	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
[9449]	168
	169	! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
	170	e3t_a(:,:,:) = e3t_n(:,:,:)
	171	e3u_a(:,:,:) = e3u_n(:,:,:)
	172	e3v_a(:,:,:) = e3v_n(:,:,:)
[6140]	173	!
	174	! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep)
[7753]	175	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration)
	176	gdepw_n(:,:,1) = 0.0_wp
	177	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg
	178	gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
	179	gdepw_b(:,:,1) = 0.0_wp
[6140]	180	DO jk = 2, jpk ! vertical sum
[5120]	181	DO jj = 1,jpj
	182	DO ji = 1,jpi
[6140]	183	! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
	184	! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
	185	! ! 0.5 where jk = mikt
	186	!!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ??
	187	zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
	188	gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
	189	gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) &
	190	& + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk))
	191	gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
	192	gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
	193	gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) &
	194	& + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk))
[4990]	195	END DO
	196	END DO
[592]	197	END DO
[6140]	198	!
	199	! !== thickness of the water column !! (ocean portion only)
[7753]	200	ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) ....
	201	hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
	202	hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
	203	hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
	204	hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
[6140]	205	DO jk = 2, jpkm1
[7753]	206	ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
	207	hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
	208	hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
	209	hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
	210	hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
[4370]	211	END DO
[6140]	212	!
	213	! !== inverse of water column thickness ==! (u- and v- points)
[7753]	214	r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF
	215	r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
	216	r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
	217	r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
	218
[6140]	219	! !== z_tilde coordinate case ==! (Restoring frequencies)
[10116]	220	tildemask(:,:) = 1._wp
	221
[4292]	222	IF( ln_vvl_ztilde ) THEN
[6140]	223	!!gm : idea: add here a READ in a file of custumized restoring frequency
[10116]	224	! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings
	225	frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp )
	226	frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
[6140]	227	!
[10116]	228	IF( ln_vvl_ztilde_as_zstar ) THEN
	229	! Ignore namelist settings and use these next two to emulate z-star using z-tilde
	230	frq_rst_e3t(:,:) = 0.0_wp
	231	frq_rst_hdv(:,:) = 1.0_wp / rdt
	232	tildemask(:,:) = 0._wp
[4292]	233	ENDIF
[10116]	234
	235	IF ( ln_vvl_zstar_at_eqtor ) THEN
[4292]	236	DO jj = 1, jpj
	237	DO ji = 1, jpi
[10116]	238	!!gm case \|gphi\| >= 6 degrees is useless initialized just above by default
[4292]	239	IF( ABS(gphit(ji,jj)) >= 6.) THEN
	240	! values outside the equatorial band and transition zone (ztilde)
	241	frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp )
[10116]	242	! frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
	243
	244	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
[4292]	245	! values inside the equatorial band (ztilde as zstar)
	246	frq_rst_e3t(ji,jj) = 0.0_wp
[10116]	247	! frq_rst_hdv(ji,jj) = 1.0_wp / rdt
	248	tildemask(ji,jj) = 0._wp
	249	ELSE
	250	! values in the transition band (linearly vary from ztilde to ztilde as zstar values)
[4292]	251	frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp &
	252	& * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
	253	& * 180._wp / 3.5_wp ) )
[10116]	254	! frq_rst_hdv(ji,jj) = (1.0_wp / rdt) &
	255	! & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp &
	256	! & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
	257	! & * 180._wp / 3.5_wp ) )
	258	tildemask(ji,jj) = 0.5_wp * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
	259	& * 180._wp / 3.5_wp ) )
[4292]	260	ENDIF
	261	END DO
	262	END DO
	263	ENDIF
[10116]	264	!
	265	IF ( ln_vvl_zstar_on_shelf ) THEN
[10164]	266	zhmin = 50._wp
[10116]	267	zhmax = 100._wp
	268	DO jj = 1, jpj
	269	DO ji = 1, jpi
	270	zwgt = 1._wp
	271	IF(( ht_0(ji,jj)>zhmin).AND.(ht_0(ji,jj) <=zhmax)) THEN
	272	zwgt = (ht_0(ji,jj)-zhmin)/(zhmax-zhmin)
	273	ELSEIF ( ht_0(ji,jj)<=zhmin) THEN
	274	zwgt = 0._wp
	275	ENDIF
	276	frq_rst_e3t(ji,jj) = MIN(frq_rst_e3t(ji,jj), frq_rst_e3t(ji,jj)*zwgt)
	277	tildemask(ji,jj) = MIN(tildemask(ji,jj), zwgt)
	278	END DO
	279	END DO
	280	ENDIF
	281	!
	282	ztmp = MAXVAL( frq_rst_hdv(:,:) )
	283	IF( lk_mpp ) CALL mpp_max( ztmp ) ! max over the global domain
	284	!
	285	IF ( (ztmp*rdt) > 1._wp) CALL ctl_stop( 'dom_vvl_init: rn_lf_cuttoff is too small' )
	286	!
[4292]	287	ENDIF
[10164]	288
	289	IF( ln_vvl_layer ) THEN
	290	IF ( ln_vvl_zstar_on_shelf ) THEN
	291	zhmin = 100._wp
	292	zhmax = 150._wp
	293	DO jj = 1, jpj
	294	DO ji = 1, jpi
	295	zwgt = 1._wp
	296	IF(( ht_0(ji,jj)>zhmin).AND.(ht_0(ji,jj) <=zhmax)) THEN
	297	zwgt = (ht_0(ji,jj)-zhmin)/(zhmax-zhmin)
	298	ELSEIF ( ht_0(ji,jj)<=zhmin) THEN
	299	zwgt = 0._wp
	300	ENDIF
	301	tildemask(ji,jj) = MIN(tildemask(ji,jj), zwgt)
	302	END DO
	303	END DO
	304	ENDIF
	305	IF ( ln_vvl_zstar_at_eqtor ) THEN
	306	DO jj = 1, jpj
	307	DO ji = 1, jpi
	308	!!gm case \|gphi\| >= 6 degrees is useless initialized just above by default
	309	IF( ABS(gphit(ji,jj)) >= 6.) THEN
	310	! values outside the equatorial band and transition zone (ztilde)
	311
	312	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
	313	! values inside the equatorial band (ztilde as zstar)
	314	tildemask(ji,jj) = 0._wp
	315	ELSE
	316	tildemask(ji,jj) = 0.5_wp * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
	317	& * 180._wp / 3.5_wp ) )
	318	ENDIF
	319	END DO
	320	END DO
	321	ENDIF
	322	ENDIF
[6140]	323	!
[9367]	324	IF(lwxios) THEN
	325	! define variables in restart file when writing with XIOS
	326	CALL iom_set_rstw_var_active('e3t_b')
	327	CALL iom_set_rstw_var_active('e3t_n')
	328	! ! ----------------------- !
	329	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
	330	! ! ----------------------- !
	331	CALL iom_set_rstw_var_active('tilde_e3t_b')
	332	CALL iom_set_rstw_var_active('tilde_e3t_n')
	333	END IF
	334	! ! -------------!
	335	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
	336	! ! ------------ !
[10116]	337	CALL iom_set_rstw_var_active('un_lf')
	338	CALL iom_set_rstw_var_active('vn_lf')
	339	CALL iom_set_rstw_var_active('hdivn_lf')
[9367]	340	ENDIF
	341	!
	342	ENDIF
	343	!
[4292]	344	END SUBROUTINE dom_vvl_init
	345
	346
[4338]	347	SUBROUTINE dom_vvl_sf_nxt( kt, kcall )
[4292]	348	!!----------------------------------------------------------------------
	349	!! * ROUTINE dom_vvl_sf_nxt *
	350	!!
	351	!! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt,
	352	!! tranxt and dynspg routines
	353	!!
	354	!! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness.
	355	!! - z_tilde_case: after scale factor increment =
	356	!! high frequency part of horizontal divergence
	357	!! + retsoring towards the background grid
	358	!! + thickness difusion
	359	!! Then repartition of ssh INCREMENT proportionnaly
	360	!! to the "baroclinic" level thickness.
	361	!!
[10116]	362	!! ** Action : - hdivn_lf : restoring towards full baroclinic divergence in z_tilde case
[4292]	363	!! - tilde_e3t_a: after increment of vertical scale factor
	364	!! in z_tilde case
[6140]	365	!! - e3(t/u/v)_a
[4292]	366	!!
	367	!! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling.
	368	!!----------------------------------------------------------------------
[6140]	369	INTEGER, INTENT( in ) :: kt ! time step
	370	INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence
	371	!
[10116]	372	LOGICAL :: ll_do_bclinic ! local logical
[10164]	373	INTEGER :: ji, jj, jk, jo ! dummy loop indices
[10116]	374	INTEGER :: ib, ib_bdy, ip, jp ! " " "
[6140]	375	INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers
[10116]	376	INTEGER :: ncall
	377	REAL(wp) :: z2dt , z_tmin, z_tmax! local scalars
	378	REAL(wp) :: zalpha, zwgt ! temporary scalars
[10164]	379	REAL(wp) :: zdu, zdv, zramp, zmet
[10116]	380	REAL(wp) :: ztra, zbtr, ztout, ztin, zfac, zmsku, zmskv
[9019]	381	REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv
[10116]	382	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, ztu, ztv
[4292]	383	!!----------------------------------------------------------------------
[6140]	384	!
	385	IF( ln_linssh ) RETURN ! No calculation in linear free surface
	386	!
[9019]	387	IF( ln_timing ) CALL timing_start('dom_vvl_sf_nxt')
[6140]	388	!
	389	IF( kt == nit000 ) THEN
[4292]	390	IF(lwp) WRITE(numout,*)
	391	IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
	392	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
	393	ENDIF
	394
[10164]	395	zmet = 1._wp
	396
[4338]	397	ll_do_bclinic = .TRUE.
[10116]	398	ncall = 1
	399
[4338]	400	IF( PRESENT(kcall) ) THEN
[10116]	401	! comment line below if tilda coordinate has to be computed at each call
[10164]	402	IF (kcall == 2 .AND. (ln_vvl_ztilde.OR.ln_vvl_layer) ) ll_do_bclinic = .FALSE.
[10116]	403	ncall = kcall
[4338]	404	ENDIF
	405
[10116]	406	IF( neuler == 0 .AND. kt == nit000 ) THEN
	407	z2dt = rdt
	408	ELSE
	409	z2dt = 2.0_wp * rdt
	410	ENDIF
	411
[4292]	412	! ******************************* !
[10116]	413	! After scale factors at t-points !
[4292]	414	! ******************************* !
[4338]	415	! ! --------------------------------------------- !
[6140]	416	! ! z_star coordinate and barotropic z-tilde part !
[4338]	417	! ! --------------------------------------------- !
[6140]	418	!
[7753]	419	z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
[4338]	420	DO jk = 1, jpkm1
[7753]	421	! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
	422	e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
[4338]	423	END DO
[6140]	424	!
[10164]	425	IF ((ln_vvl_ztilde.OR.ln_vvl_layer).AND.(zmet==1._wp)) THEN
	426	DO jk = 1, jpkm1
	427	e3t_a(:,:,jk) = e3t_a(:,:,jk) - tilde_e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
	428	END DO
	429	ENDIF
	430
[10116]	431	IF( (ln_vvl_ztilde.OR.ln_vvl_layer) .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !
[10164]	432	! ! ------baroclinic part------ !
[10116]	433	tilde_e3t_a(:,:,:) = 0.0_wp ! tilde_e3t_a used to store tendency terms
	434	un_td(:,:,:) = 0.0_wp ! Transport corrections
	435	vn_td(:,:,:) = 0.0_wp
[4292]	436
[10116]	437	zhdiv(:,:) = 0.
[4292]	438	DO jk = 1, jpkm1
[7753]	439	zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
[592]	440	END DO
[10116]	441	zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) )
[2528]	442
[10116]	443	ze3t(:,:,:) = 0._wp
	444	IF( ln_rnf ) THEN
	445	CALL sbc_rnf_div( ze3t ) ! runoffs
	446	DO jk=1,jpkm1
	447	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ze3t(:,:,jk) * e3t_n(:,:,jk)
	448	END DO
	449	ENDIF
	450
	451	! Thickness advection:
	452	! --------------------
	453	! Set advection velocities and source term
	454	IF ( ln_vvl_ztilde ) THEN
	455	IF ( ncall==1 ) THEN
	456	zalpha = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
[4292]	457	DO jk = 1, jpkm1
[10116]	458	ztu(:,:,jk) = un(:,:,jk) * e3u_n(:,:,jk) * e2u(:,:)
	459	ztv(:,:,jk) = vn(:,:,jk) * e3v_n(:,:,jk) * e1v(:,:)
[10164]	460	ze3t(:,:,jk) = -tilde_e3t_a(:,:,jk) - (e3t_n(:,:,jk)-zmettilde_e3t_n(:,:,jk)) zhdiv(:,:)
[4292]	461	END DO
[10116]	462	!
[10164]	463	un_lf(:,:,:,nn_filt_order) = un_lf(:,:,:,nn_filt_order) * (1._wp - zalpha) + zalpha * ztu(:,:,:)
	464	vn_lf(:,:,:,nn_filt_order) = vn_lf(:,:,:,nn_filt_order) * (1._wp - zalpha) + zalpha * ztv(:,:,:)
	465	hdivn_lf(:,:,:,nn_filt_order) = hdivn_lf(:,:,:,nn_filt_order) * (1._wp - zalpha) + zalpha * ze3t(:,:,:)
	466
	467	DO jo = nn_filt_order-1,1,-1
	468	un_lf(:,:,:,jo) = un_lf(:,:,:,jo) * (1._wp - zalpha) + zalpha * un_lf(:,:,:,jo+1)
	469	vn_lf(:,:,:,jo) = vn_lf(:,:,:,jo) * (1._wp - zalpha) + zalpha * vn_lf(:,:,:,jo+1)
	470	hdivn_lf(:,:,:,jo) = hdivn_lf(:,:,:,jo) * (1._wp - zalpha) + zalpha * hdivn_lf(:,:,:,jo+1)
	471	END DO
[4292]	472	ENDIF
[10116]	473	!
	474	DO jk = 1, jpkm1
[10164]	475	ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk,1)/e3u_n(:,:,jk)r1_e2u(:,:))umask(:,:,jk)
	476	ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk,1)/e3v_n(:,:,jk)r1_e1v(:,:))vmask(:,:,jk)
	477	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk,1) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
[10116]	478	END DO
	479	!
	480	ELSEIF ( ln_vvl_layer ) THEN
	481	!
[4292]	482	DO jk = 1, jpkm1
[10116]	483	ztu(:,:,jk) = un(:,:,jk)
	484	ztv(:,:,jk) = vn(:,:,jk)
[4292]	485	END DO
[10116]	486	!
	487	ENDIF
	488	!
	489	! Block fluxes through small layers:
[10164]	490	! DO jk=1,jpkm1
	491	! DO ji = 1, jpi
	492	! DO jj= 1, jpj
	493	! zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, e3u_n(ji,jj,jk) - hsmall) )
	494	! un_td(ji,jj,jk) = un_td(ji,jj,jk) - (1. - zmsku) * un(ji,jj,jk) * e3u_n(ji,jj,jk) * e2u(ji,jj)
	495	! ztu(ji,jj,jk) = zmsku * ztu(ji,jj,jk)
	496	! IF ( ln_vvl_ztilde ) un_lf(ji,jj,jk) = zmsku * un_lf(ji,jj,jk)
	497	! !
	498	! zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, e3v_n(ji,jj,jk) - hsmall) )
	499	! vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - (1. - zmskv) * vn(ji,jj,jk) * e3v_n(ji,jj,jk) * e1v(ji,jj)
	500	! ztv(ji,jj,jk) = zmskv * ztv(ji,jj,jk)
	501	! IF ( ln_vvl_ztilde ) vn_lf(ji,jj,jk) = zmskv * vn_lf(ji,jj,jk)
	502	! END DO
	503	! END DO
	504	! END DO
[10116]	505	!
	506	! Do advection
	507	IF (ln_vvl_adv_fct) THEN
	508	CALL dom_vvl_adv_fct( kt, tilde_e3t_a, ztu, ztv )
	509	!
	510	ELSEIF (ln_vvl_adv_cn2) THEN
[4292]	511	DO jk = 1, jpkm1
[10116]	512	DO jj = 2, jpjm1
	513	DO ji = fs_2, fs_jpim1 ! vector opt.
	514	tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) &
	515	& -( e2u(ji,jj)e3u_n(ji,jj,jk) ztu(ji,jj,jk) - e2u(ji-1,jj )e3u_n(ji-1,jj ,jk) ztu(ji-1,jj ,jk) &
	516	& + e1v(ji,jj)e3v_n(ji,jj,jk) ztv(ji,jj,jk) - e1v(ji ,jj-1)e3v_n(ji ,jj-1,jk) ztv(ji ,jj-1,jk) ) &
	517	& / ( e1t(ji,jj) * e2t(ji,jj) )
	518	END DO
	519	END DO
[4292]	520	END DO
[6140]	521	ENDIF
[10116]	522	!
	523	! Thickness anomaly diffusion:
[10164]	524	! ----------------------------
[10116]	525	ztu(:,:,:) = 0.0_wp
	526	ztv(:,:,:) = 0.0_wp
[4292]	527
[10116]	528	ze3t(:,:,1) = 0.e0
	529	DO jj = 1, jpj
	530	DO ji = 1, jpi
	531	DO jk = 2, jpk
	532	ze3t(ji,jj,jk) = ze3t(ji,jj,jk-1) + tilde_e3t_b(ji,jj,jk-1) * tmask(ji,jj,jk-1)
	533	END DO
[4292]	534	END DO
[10116]	535	END DO
	536
	537	IF ( ln_vvl_blp ) THEN ! Bilaplacian
	538	DO jk = 1, jpkm1
	539	DO jj = 1, jpjm1 ! First derivative (gradient)
	540	DO ji = 1, fs_jpim1 ! vector opt.
	541	ztu(ji,jj,jk) = umask(ji,jj,jk) * e2_e1u(ji,jj) &
	542	& * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj ,jk) )
	543	ztv(ji,jj,jk) = vmask(ji,jj,jk) * e1_e2v(ji,jj) &
	544	& * ( ze3t(ji,jj,jk) - ze3t(ji ,jj+1,jk) )
	545	END DO
	546	END DO
	547
	548	DO jj = 2, jpjm1 ! Second derivative (divergence) time the eddy diffusivity coefficient
	549	DO ji = fs_2, fs_jpim1 ! vector opt.
	550	zht(ji,jj) = rn_ahe3_blp * r1_e1e2t(ji,jj) * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) &
	551	& + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )
	552	END DO
	553	END DO
	554
	555	! Open boundary conditions:
	556	IF ( ln_bdy ) THEN
	557	DO ib_bdy=1, nb_bdy
	558	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
	559	ji = idx_bdy(ib_bdy)%nbi(ib,1)
	560	jj = idx_bdy(ib_bdy)%nbj(ib,1)
	561	zht(ji,jj) = 0._wp
	562	END DO
	563	END DO
	564	END IF
	565
	566	CALL lbc_lnk( zht, 'T', 1. ) ! Lateral boundary conditions (unchanged sgn)
	567
	568	DO jj = 1, jpjm1 ! third derivative (gradient)
	569	DO ji = 1, fs_jpim1 ! vector opt.
	570	ztu(ji,jj,jk) = umask(ji,jj,jk) * e2_e1u(ji,jj) * ( zht(ji+1,jj ) - zht(ji,jj) )
	571	ztv(ji,jj,jk) = vmask(ji,jj,jk) * e1_e2v(ji,jj) * ( zht(ji ,jj+1) - zht(ji,jj) )
	572	END DO
	573	END DO
	574	END DO
[6140]	575	ENDIF
[4292]	576
[10116]	577	IF ( ln_vvl_lap ) THEN ! Laplacian
	578	DO jk = 1, jpkm1 ! First derivative (gradient)
	579	DO jj = 1, jpjm1
	580	DO ji = 1, fs_jpim1 ! vector opt.
	581	zdu = rn_ahe3_lap * umask(ji,jj,jk) * e2_e1u(ji,jj) &
	582	& * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj ,jk) )
	583	zdv = rn_ahe3_lap * vmask(ji,jj,jk) * e1_e2v(ji,jj) &
	584	& * ( ze3t(ji,jj,jk) - ze3t(ji ,jj+1,jk) )
	585	ztu(ji,jj,jk) = ztu(ji,jj,jk) + zdu
	586	ztv(ji,jj,jk) = ztv(ji,jj,jk) + zdv
	587	END DO
[4292]	588	END DO
	589	END DO
[10116]	590	ENDIF
	591
	592	! divergence of diffusive fluxes
	593	DO jk = 1, jpkm1
[4292]	594	DO jj = 2, jpjm1
	595	DO ji = fs_2, fs_jpim1 ! vector opt.
[10116]	596	un_td(ji,jj,jk) = un_td(ji,jj,jk) + ztu(ji,jj,jk+1) - ztu(ji,jj,jk )
	597	vn_td(ji,jj,jk) = vn_td(ji,jj,jk) + ztv(ji,jj,jk+1) - ztv(ji,jj,jk )
	598	tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( ztu(ji-1,jj ,jk+1) - ztu(ji,jj,jk+1) &
	599	& +ztv(ji ,jj-1,jk+1) - ztv(ji,jj,jk+1) &
	600	& -ztu(ji-1,jj ,jk ) + ztu(ji,jj,jk ) &
	601	& -ztv(ji ,jj-1,jk ) + ztv(ji,jj,jk ) &
[5836]	602	& ) * r1_e1e2t(ji,jj)
[4292]	603	END DO
	604	END DO
	605	END DO
	606
[10128]	607	CALL lbc_lnk_multi( un_td, 'U', -1., vn_td, 'V', -1. ) !* local domain boundaries
[10116]	608	!
	609	CALL dom_vvl_ups_cor( kt, tilde_e3t_a, un_td, vn_td )
	610
	611	! IF ( ln_vvl_ztilde ) THEN
	612	! ztu(:,:,:) = -un_lf(:,:,:)
	613	! ztv(:,:,:) = -vn_lf(:,:,:)
	614	! CALL dom_vvl_ups_cor( kt, tilde_e3t_a, ztu, ztv )
	615	! ENDIF
	616	!
	617	! Remove "external thickness" tendency:
	618	DO jk = 1, jpkm1
[10164]	619	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + (e3t_n(:,:,jk)-zmettilde_e3t_n(:,:,jk)) zhdiv(:,:) * tmask(:,:,jk)
[10116]	620	END DO
	621
[4292]	622	! Leapfrog time stepping
	623	! ~~~~~~~~~~~~~~~~~~~~~~
[10116]	624	zramp = 1._wp
	625	IF ((.NOT.ln_rstart).AND.ln_vvl_ramp) zramp = MIN(MAX( ((kt-nit000)rdt)/(rn_day_ramprday),0._wp),1._wp)
[4292]	626
[10116]	627	DO jk=1,jpkm1
	628	tilde_e3t_a(:,:,jk) = tilde_e3t_b(:,:,jk) + z2dt * tmask(:,:,jk) * tilde_e3t_a(:,:,jk) &
	629	& * tildemask(:,:) * zramp
[4292]	630	END DO
[10116]	631
	632	! Ensure layer separation:
	633	! ~~~~~~~~~~~~~~~~~~~~~~~~
	634	IF ( ln_vvl_regrid ) CALL dom_vvl_regrid( kt )
	635
	636	! Boundary conditions:
	637	! ~~~~~~~~~~~~~~~~~~~~
	638	IF ( ln_bdy ) THEN
	639	DO ib_bdy = 1, nb_bdy
	640	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
	641	!! DO ib = 1, idx_bdy(ib_bdy)%nblen(1)
	642	ji = idx_bdy(ib_bdy)%nbi(ib,1)
	643	jj = idx_bdy(ib_bdy)%nbj(ib,1)
	644	zwgt = idx_bdy(ib_bdy)%nbw(ib,1)
	645	ip = bdytmask(ji+1,jj ) - bdytmask(ji-1,jj )
	646	jp = bdytmask(ji ,jj+1) - bdytmask(ji ,jj-1)
	647	DO jk = 1, jpkm1
	648	tilde_e3t_a(ji,jj,jk) = 0.e0
	649	!! tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) * (1._wp - zwgt)
	650	!! tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji+ip,jj+jp,jk) * tmask(ji+ip,jj+jp,jk)
	651	END DO
	652	END DO
	653	END DO
[4292]	654	ENDIF
	655
[10116]	656	CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )
	657
	658	ENDIF
	659
	660	IF( ln_vvl_ztilde.AND.( ncall==1)) THEN
	661	zalpha = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
[4338]	662	!
[10116]	663	! divergence of diffusive fluxes
[4338]	664	DO jk = 1, jpkm1
[10116]	665	DO jj = 2, jpjm1
	666	DO ji = fs_2, fs_jpim1 ! vector opt.
	667	ze3t(ji,jj,jk) = ( un_td(ji,jj,jk) - un_td(ji-1,jj ,jk) &
	668	& + vn_td(ji,jj,jk) - vn_td(ji ,jj-1,jk) &
	669	& ) / ( e1t(ji,jj) * e2t(ji,jj) )
	670	END DO
	671	END DO
[4338]	672	END DO
[10116]	673	CALL lbc_lnk( ze3t(:,:,:), 'T', 1. )
[10164]	674
	675	DO jo = nn_filt_order,1,-1
	676	hdivn_lf(:,:,:,jo) = hdivn_lf(:,:,:,jo) + zalpha*(nn_filt_order-jo+1) ze3t(:,:,:)
	677	END DO
[4292]	678	ENDIF
	679
[4338]	680	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate !
[10116]	681	! ! ---baroclinic part--------- !
[10164]	682
	683	IF ( (ncall==2).AND.(.NOT.ll_do_bclinic) ) THEN
	684
	685	DO jk = 1, jpkm1
	686	ztu(:,:,jk) = (un_adv(:,:)r1_hu_n(:,:) - un_b(:,:) ) e3u_n(:,:,jk) * e2u(:,:) * umask(:,:,jk)
	687	ztv(:,:,jk) = (vn_adv(:,:)r1_hv_n(:,:) - vn_b(:,:) ) e3v_n(:,:,jk) * e1v(:,:) * vmask(:,:,jk)
	688	DO jj = 2, jpjm1
	689	DO ji = fs_2, fs_jpim1 ! vector opt.
	690	ze3t(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk) &
	691	& + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) &
	692	& ) / ( e1t(ji,jj) * e2t(ji,jj) )
	693	END DO
	694	END DO
	695	END DO
	696	!
	697	zhdiv(:,:) = 0.
	698	DO jk = 1, jpkm1
	699	zhdiv(:,:) = zhdiv(:,:) + ze3t(:,:,jk) * tmask(:,:,jk)
	700	END DO
	701	zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) )
	702	!
	703	DO jk = 1, jpkm1
	704	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - z2dt * (ze3t(:,:,jk) &
	705	& - zhdiv(:,:)(e3t_n(:,:,jk)-zmettilde_e3t_n(:,:,jk))*tmask(:,:,jk))
	706	END DO
	707	CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )
	708	ENDIF
	709	!
[4338]	710	DO jk = 1, jpkm1
[10164]	711	e3t_a(:,:,jk) = e3t_a(:,:,jk) + tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
[4338]	712	END DO
[10164]	713	! IF( ln_vvl_ztilde ) THEN ! Relax barotropic component:
	714	! DO jk = 1, jpkm1
	715	! e3t_a(:,:,jk) = e3t_a(:,:,jk) &
	716	! & - z2dt * frq_rst_e3t(:,:) * (e3t_b(:,:,jk) - tilde_e3t_b(:,:,jk) &
	717	! & - e3t_0(:,:,jk) * (ht_0(:,:) + sshb(:,:))/ (ht_0(:,:)*tmask(:,:,1) + 1._wp - tmask(:,:,1)))
	718	! END DO
	719	! ENDIF
[4338]	720	ENDIF
	721
[10116]	722	IF( ln_vvl_dbg.AND.(ln_vvl_ztilde .OR. ln_vvl_layer) ) THEN ! - ML - test: control prints for debuging
[4292]	723	!
[10116]	724	zht(:,:) = 0.0_wp
	725	DO jk = 1, jpkm1
	726	zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
	727	END DO
	728	IF( lwp ) WRITE(numout, *) 'kt = ', kt
	729	IF( lwp ) WRITE(numout, *) 'ncall = ', ncall
	730	IF( lwp ) WRITE(numout, *) 'll_do_bclinic', ll_do_bclinic
[4292]	731	IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
[10116]	732	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ), mask = tmask(:,:,1) == 1.e0 )
[4292]	733	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
	734	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
	735	END IF
	736	!
[10116]	737	z_tmin = MINVAL( e3t_n(:,:,:) )
	738	IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain
	739	IF( lwp ) WRITE(numout, *) kt,' MINVAL(e3t_n) =', z_tmin
	740	IF ( z_tmin .LE. 0._wp ) THEN
	741	IF( lk_mpp ) THEN
	742	CALL mpp_minloc(e3t_n(:,:,:), tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
	743	ELSE
	744	ijk_min = MINLOC( e3t_n(:,:,:) )
	745	ijk_min(1) = ijk_min(1) + nimpp - 1
	746	ijk_min(2) = ijk_min(2) + njmpp - 1
	747	ENDIF
	748	IF (lwp) THEN
[10164]	749	ji = ijk_min(1) ; jj = ijk_min(2) ; jk = ijk_min(3)
[10116]	750	WRITE(numout, *) 'Negative scale factor, e3t_n =', z_tmin
	751	WRITE(numout, *) 'at i, j, k=', ijk_min
	752	CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
	753	ENDIF
	754	ENDIF
	755	!
	756	z_tmin = MINVAL( e3u_n(:,:,:))
	757	IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain
	758	IF( lwp ) WRITE(numout, *) kt,' MINVAL(e3u_n) =', z_tmin
	759	IF ( z_tmin .LE. 0._wp ) THEN
	760	IF( lk_mpp ) THEN
	761	CALL mpp_minloc(e3u_n(:,:,:), umask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
	762	ELSE
	763	ijk_min = MINLOC( e3u_n(:,:,:) )
	764	ijk_min(1) = ijk_min(1) + nimpp - 1
	765	ijk_min(2) = ijk_min(2) + njmpp - 1
	766	ENDIF
	767	IF (lwp) THEN
	768	WRITE(numout, *) 'Negative scale factor, e3u_n =', z_tmin
	769	WRITE(numout, *) 'at i, j, k=', ijk_min
	770	CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
	771	ENDIF
	772	ENDIF
	773	!
	774	z_tmin = MINVAL( e3v_n(:,:,:) )
	775	IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain
	776	IF( lwp ) WRITE(numout, *) kt,' MINVAL(e3v_n) =', z_tmin
	777	IF ( z_tmin .LE. 0._wp ) THEN
	778	IF( lk_mpp ) THEN
	779	CALL mpp_minloc(e3v_n(:,:,:), vmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
	780	ELSE
	781	ijk_min = MINLOC( e3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0 )
	782	ijk_min(1) = ijk_min(1) + nimpp - 1
	783	ijk_min(2) = ijk_min(2) + njmpp - 1
	784	ENDIF
	785	IF (lwp) THEN
	786	WRITE(numout, *) 'Negative scale factor, e3v_n =', z_tmin
	787	WRITE(numout, *) 'at i, j, k=', ijk_min
	788	CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
	789	ENDIF
	790	ENDIF
	791	!
	792	z_tmin = MINVAL( e3f_n(:,:,:))
	793	IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain
	794	IF( lwp ) WRITE(numout, *) kt,' MINVAL(e3f_n) =', z_tmin
	795	IF ( z_tmin .LE. 0._wp ) THEN
	796	IF( lk_mpp ) THEN
	797	CALL mpp_minloc(e3f_n(:,:,:), fmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
	798	ELSE
	799	ijk_min = MINLOC( e3f_n(:,:,:) )
	800	ijk_min(1) = ijk_min(1) + nimpp - 1
	801	ijk_min(2) = ijk_min(2) + njmpp - 1
	802	ENDIF
	803	IF (lwp) THEN
	804	WRITE(numout, *) 'Negative scale factor, e3f_n =', z_tmin
	805	WRITE(numout, *) 'at i, j, k=', ijk_min
	806	CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
	807	ENDIF
	808	ENDIF
	809	!
[7753]	810	zht(:,:) = 0.0_wp
[4292]	811	DO jk = 1, jpkm1
[7753]	812	zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
[4292]	813	END DO
	814	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
	815	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
[10116]	816	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn -SUM(e3t_n))) =', z_tmax
[4292]	817	!
[7753]	818	zht(:,:) = 0.0_wp
[4292]	819	DO jk = 1, jpkm1
[10116]	820	zht(:,:) = zht(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
[4292]	821	END DO
[10116]	822	zwu(:,:) = 0._wp
	823	DO jj = 1, jpjm1
	824	DO ji = 1, fs_jpim1 ! vector opt.
	825	zwu(ji,jj) = 0.5_wp * umask(ji,jj,1) * r1_e1e2u(ji,jj) &
	826	& * ( e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj) )
	827	END DO
	828	END DO
	829	CALL lbc_lnk( zwu(:,:), 'U', 1._wp )
	830	z_tmax = MAXVAL( ssumask(:,:) * ssumask(:,:) * ABS( hu_0(:,:) + zwu(:,:) - zht(:,:) ) )
[4292]	831	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
[10116]	832	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hu_0+sshu_n-SUM(e3u_n))) =', z_tmax
[4292]	833	!
[7753]	834	zht(:,:) = 0.0_wp
[4292]	835	DO jk = 1, jpkm1
[10116]	836	zht(:,:) = zht(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
[4292]	837	END DO
[10116]	838	zwv(:,:) = 0._wp
	839	DO jj = 1, jpjm1
	840	DO ji = 1, fs_jpim1 ! vector opt.
	841	zwv(ji,jj) = 0.5_wp * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &
	842	& * ( e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1) )
	843	END DO
	844	END DO
	845	CALL lbc_lnk( zwv(:,:), 'V', 1._wp )
	846	z_tmax = MAXVAL( ssvmask(:,:) * ssvmask(:,:) * ABS( hv_0(:,:) + zwv(:,:) - zht(:,:) ) )
[4292]	847	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
[10116]	848	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hv_0+sshv_n-SUM(e3v_n))) =', z_tmax
[4292]	849	!
[10116]	850	zht(:,:) = 0.0_wp
	851	DO jk = 1, jpkm1
	852	DO jj = 1, jpjm1
	853	zht(:,jj) = zht(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk)*umask(:,jj+1,jk)
	854	END DO
	855	END DO
	856	zwu(:,:) = 0._wp
	857	DO jj = 1, jpjm1
	858	DO ji = 1, fs_jpim1 ! vector opt.
	859	zwu(ji,jj) = 0.25_wp * umask(ji,jj,1) * umask(ji,jj+1,1) * r1_e1e2f(ji,jj) &
	860	& * ( e1e2t(ji ,jj) * sshn(ji ,jj) + e1e2t(ji ,jj+1) * sshn(ji ,jj+1) &
	861	& + e1e2t(ji+1,jj) * sshn(ji+1,jj) + e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1) )
	862	END DO
	863	END DO
	864	CALL lbc_lnk( zht(:,:), 'F', 1._wp )
	865	CALL lbc_lnk( zwu(:,:), 'F', 1._wp )
	866	z_tmax = MAXVAL( fmask(:,:,1) * fmask(:,:,1) * ABS( hf_0(:,:) + zwu(:,:) - zht(:,:) ) )
[4292]	867	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
[10116]	868	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hf_0+sshf_n-SUM(e3f_n))) =', z_tmax
[4292]	869	!
	870	END IF
	871
	872	! *********************************** !
	873	! After scale factors at u- v- points !
	874	! *********************************** !
	875
[6140]	876	CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
	877	CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
[4292]	878
[4370]	879	! *********************************** !
	880	! After depths at u- v points !
	881	! *********************************** !
	882
[7753]	883	hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
	884	hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
[6140]	885	DO jk = 2, jpkm1
[7753]	886	hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
	887	hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
[4370]	888	END DO
	889	! ! Inverse of the local depth
[6140]	890	!!gm BUG ? don't understand the use of umask_i here .....
[7753]	891	r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
	892	r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
[6140]	893	!
[9019]	894	IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt')
[6140]	895	!
[4292]	896	END SUBROUTINE dom_vvl_sf_nxt
	897
	898
	899	SUBROUTINE dom_vvl_sf_swp( kt )
[3294]	900	!!----------------------------------------------------------------------
[4292]	901	!! * ROUTINE dom_vvl_sf_swp *
[3294]	902	!!
[4292]	903	!! ** Purpose : compute time filter and swap of scale factors
	904	!! compute all depths and related variables for next time step
	905	!! write outputs and restart file
[3294]	906	!!
[4292]	907	!! ** Method : - swap of e3t with trick for volume/tracer conservation
	908	!! - reconstruct scale factor at other grid points (interpolate)
	909	!! - recompute depths and water height fields
	910	!!
[6140]	911	!! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
[4292]	912	!! - Recompute:
[6140]	913	!! e3(u/v)_b
	914	!! e3w_n
	915	!! e3(u/v)w_b
	916	!! e3(u/v)w_n
	917	!! gdept_n, gdepw_n and gde3w_n
[4292]	918	!! h(u/v) and h(u/v)r
	919	!!
	920	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
	921	!! Leclair, M., and G. Madec, 2011, Ocean Modelling.
[3294]	922	!!----------------------------------------------------------------------
[6140]	923	INTEGER, INTENT( in ) :: kt ! time step
	924	!
	925	INTEGER :: ji, jj, jk ! dummy loop indices
	926	REAL(wp) :: zcoef ! local scalar
[3294]	927	!!----------------------------------------------------------------------
[6140]	928	!
	929	IF( ln_linssh ) RETURN ! No calculation in linear free surface
	930	!
[9019]	931	IF( ln_timing ) CALL timing_start('dom_vvl_sf_swp')
[3294]	932	!
[4292]	933	IF( kt == nit000 ) THEN
	934	IF(lwp) WRITE(numout,*)
	935	IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
	936	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step'
[3294]	937	ENDIF
[4292]	938	!
	939	! Time filter and swap of scale factors
	940	! =====================================
[6140]	941	! - ML - e3(t/u/v)_b are allready computed in dynnxt.
[4292]	942	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
	943	IF( neuler == 0 .AND. kt == nit000 ) THEN
[7753]	944	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
[4292]	945	ELSE
[7753]	946	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &
	947	& + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
[4292]	948	ENDIF
[7753]	949	tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
[4292]	950	ENDIF
[7753]	951	gdept_b(:,:,:) = gdept_n(:,:,:)
	952	gdepw_b(:,:,:) = gdepw_n(:,:,:)
[4488]	953
[7753]	954	e3t_n(:,:,:) = e3t_a(:,:,:)
	955	e3u_n(:,:,:) = e3u_a(:,:,:)
	956	e3v_n(:,:,:) = e3v_a(:,:,:)
	957
[4292]	958	! Compute all missing vertical scale factor and depths
	959	! ====================================================
	960	! Horizontal scale factor interpolations
	961	! --------------------------------------
[6140]	962	! - ML - e3u_b and e3v_b are allready computed in dynnxt
[4370]	963	! - JC - hu_b, hv_b, hur_b, hvr_b also
[6140]	964
[10060]	965	CALL dom_vvl_interpol( e3t_n(:,:,:), e3f_n(:,:,:), 'F' )
[6140]	966
[4292]	967	! Vertical scale factor interpolations
[6140]	968	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' )
	969	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
	970	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
	971	CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b(:,:,:), 'W' )
	972	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
	973	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
[5120]	974
[6140]	975	! t- and w- points depth (set the isf depth as it is in the initial step)
[7753]	976	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
	977	gdepw_n(:,:,1) = 0.0_wp
	978	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
[5120]	979	DO jk = 2, jpk
	980	DO jj = 1,jpj
	981	DO ji = 1,jpi
	982	! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
	983	! 1 for jk = mikt
	984	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
[6140]	985	gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
	986	gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) &
	987	& + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) )
	988	gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
[4990]	989	END DO
	990	END DO
[4292]	991	END DO
[5120]	992
[6140]	993	! Local depth and Inverse of the local depth of the water
	994	! -------------------------------------------------------
[7753]	995	hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:)
	996	hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:)
	997	!
	998	ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
[6140]	999	DO jk = 2, jpkm1
[7753]	1000	ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
[4370]	1001	END DO
[7753]	1002
[4292]	1003	! write restart file
	1004	! ==================
[9019]	1005	IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
[4292]	1006	!
[9019]	1007	IF( ln_timing ) CALL timing_stop('dom_vvl_sf_swp')
[6140]	1008	!
[4292]	1009	END SUBROUTINE dom_vvl_sf_swp
	1010
	1011
	1012	SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
	1013	!!---------------------------------------------------------------------
	1014	!! * ROUTINE dom_vvl__interpol *
	1015	!!
	1016	!! ** Purpose : interpolate scale factors from one grid point to another
	1017	!!
	1018	!! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0)
	1019	!! - horizontal interpolation: grid cell surface averaging
	1020	!! - vertical interpolation: simple averaging
	1021	!!----------------------------------------------------------------------
[5836]	1022	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3_in ! input e3 to be interpolated
	1023	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3_out ! output interpolated e3
	1024	CHARACTER(LEN=*) , INTENT(in ) :: pout ! grid point of out scale factors
	1025	! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'
	1026	!
[10060]	1027	INTEGER :: ji, jj, jk, jkbot ! dummy loop indices
	1028	INTEGER :: nmet ! horizontal interpolation method
[9023]	1029	REAL(wp) :: zlnwd ! =1./0. when ln_wd_il = T/F
[10060]	1030	REAL(wp) :: ztap, zsmall ! Parameters defining minimum thicknesses UVF-points
	1031	REAL(wp) :: zmin
	1032	REAL(wp) :: zdo, zup ! Lower and upper interfaces depths anomalies
	1033	REAL(wp), DIMENSION(jpi,jpj) :: zs ! Surface interface depth anomaly
	1034	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw ! Interface depth anomaly
[4292]	1035	!!----------------------------------------------------------------------
[5836]	1036	!
[10116]	1037	! nmet = 0 ! Original method (Surely wrong)
	1038	nmet = 2 ! Interface interpolation
	1039	! nmet = 2 ! Internal interfaces interpolation only, spread barotropic increment
[10060]	1040	! Note that we kept surface weighted interpolation for barotropic increment to be compliant
	1041	! with what is done in surface pressure module.
	1042	!
[9023]	1043	IF(ln_wd_il) THEN
[6152]	1044	zlnwd = 1.0_wp
	1045	ELSE
	1046	zlnwd = 0.0_wp
	1047	END IF
	1048	!
[10116]	1049	ztap = 0._wp ! Minimum fraction of T-point thickness at cell interfaces
	1050	zsmall = 1.e-8_wp ! Minimum thickness at U or V points (m)
[10060]	1051	!
	1052	IF ( (nmet==1).OR.(nmet==2) ) THEN
	1053	SELECT CASE ( pout )
	1054	!
	1055	CASE( 'U', 'V', 'F' )
	1056	! Compute interface depth anomaly at T-points
	1057	!
	1058	zw(:,:,:) = 0._wp
	1059	!
	1060	DO jk=2,jpk
	1061	zw(:,:,jk) = zw(:,:,jk-1) + pe3_in(:,:,jk-1)*tmask(:,:,jk-1)
	1062	END DO
	1063	! Interface depth anomalies:
	1064	DO jk=1,jpkm1
	1065	zw(:,:,jk) = zw(:,:,jk) - zw(:,:,jpk) + ht_0(:,:)
	1066	END DO
	1067	zw(:,:,jpk) = ht_0(:,:)
	1068	!
	1069	IF (nmet==2) THEN ! Consider "internal" interfaces only
	1070	zs(:,:) = - zw(:,:,1) ! Save surface anomaly (ssh)
	1071	!
	1072	DO jj = 1, jpj
	1073	DO ji = 1, jpi
	1074	DO jk=1,jpk
	1075	zw(ji,jj,jk) = (zw(ji,jj,jk) + zs(ji,jj)) &
	1076	& * ht_0(ji,jj) / (ht_0(ji,jj) + zs(ji,jj) + 1._wp - tmask(ji,jj,1)) &
	1077	& * tmask(ji,jj,jk)
	1078	END DO
	1079	END DO
	1080	END DO
	1081	ENDIF
[10116]	1082	zw(:,:,:) = (zw(:,:,:) - gdepw_0(:,:,:))*tmask(:,:,:)
[10060]	1083	!
	1084	END SELECT
	1085	END IF
	1086	!
	1087	pe3_out(:,:,:) = 0.0_wp
	1088	!
[5836]	1089	SELECT CASE ( pout ) !== type of interpolation ==!
[4292]	1090	!
[5836]	1091	CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean
[10060]	1092	IF (nmet==0) THEN
	1093	DO jk = 1, jpk
	1094	DO jj = 1, jpjm1
	1095	DO ji = 1, fs_jpim1 ! vector opt.
	1096	pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) &
	1097	& * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	1098	& + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
	1099	END DO
	1100	END DO
	1101	END DO
	1102	ELSE
[4292]	1103	DO jj = 1, jpjm1
	1104	DO ji = 1, fs_jpim1 ! vector opt.
[10060]	1105	! Correction at last level:
	1106	jkbot = mbku(ji,jj)
[10116]	1107	zdo = 0._wp
[10060]	1108	DO jk=jkbot,1,-1
[10116]	1109	zup = 0.5_wp * ( e1e2t(ji ,jj)*zw(ji ,jj,jk) &
	1110	& + e1e2t(ji+1,jj)zw(ji+1,jj,jk) ) r1_e1e2u(ji,jj)
[10060]	1111	!
	1112	! If there is a step, taper bottom interface:
[10116]	1113	! IF ((hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ).AND.(jk==jkbot)) THEN
	1114	! IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN
	1115	! zmin = ztap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk))
	1116	! ELSE
	1117	! zmin = ztap * (zw(ji ,jj,jk+1)-zw(ji ,jj,jk))
	1118	! ENDIF
	1119	! zup = MIN(zup, zdo-zmin)
	1120	! ENDIF
	1121	zup = MIN(zup, zdo+e3u_0(ji,jj,jk)-zsmall)
	1122	pe3_out(ji,jj,jk) = (zdo - zup) * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )
[10060]	1123	zdo = zup
	1124	END DO
[4292]	1125	END DO
[2528]	1126	END DO
[10060]	1127	END IF
	1128	!
	1129	IF (nmet==2) THEN ! Spread sea level anomaly
	1130	DO jj = 1, jpjm1
	1131	DO ji = 1, fs_jpim1 ! vector opt.
	1132	DO jk=1,jpk
	1133	pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) &
	1134	& + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) ) &
[10116]	1135	& / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) ) &
[10060]	1136	& * 0.5_wp * r1_e1e2u(ji,jj) &
	1137	& * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) &
	1138	& * ( e1e2t(ji,jj)zs(ji,jj) + e1e2t(ji+1,jj)zs(ji+1,jj) )
	1139	END DO
	1140	END DO
	1141	END DO
	1142	!
	1143	ENDIF
	1144	!
[4990]	1145	CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
[7753]	1146	pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
[5836]	1147	!
	1148	CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean
[10060]	1149	IF (nmet==0) THEN
	1150	DO jk = 1, jpk
	1151	DO jj = 1, jpjm1
	1152	DO ji = 1, fs_jpim1 ! vector opt.
	1153	pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) &
	1154	& * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
	1155	& + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
	1156	END DO
	1157	END DO
	1158	END DO
	1159	ELSE
[4292]	1160	DO jj = 1, jpjm1
	1161	DO ji = 1, fs_jpim1 ! vector opt.
[10060]	1162	! Correction at last level:
	1163	jkbot = mbkv(ji,jj)
[10116]	1164	zdo = 0._wp
[10060]	1165	DO jk=jkbot,1,-1
[10116]	1166	zup = 0.5_wp * ( e1e2t(ji,jj ) * zw(ji,jj ,jk) &
	1167	& + e1e2t(ji,jj+1) * zw(ji,jj+1,jk) ) * r1_e1e2v(ji,jj)
[10060]	1168	!
	1169	! If there is a step, taper bottom interface:
[10116]	1170	! IF ((hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ).AND.(jk==jkbot)) THEN
	1171	! IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN
	1172	! zmin = ztap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk))
	1173	! ELSE
	1174	! zmin = ztap * (zw(ji ,jj,jk+1)-zw(ji ,jj,jk))
	1175	! ENDIF
	1176	! zup = MIN(zup, zdo-zmin)
	1177	! ENDIF
	1178	zup = MIN(zup, zdo + e3v_0(ji,jj,jk) - zsmall)
	1179	pe3_out(ji,jj,jk) = (zdo - zup) * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )
[10060]	1180	zdo = zup
	1181	END DO
[4292]	1182	END DO
	1183	END DO
[10060]	1184	END IF
	1185	!
	1186	IF (nmet==2) THEN ! Spread sea level anomaly
	1187	DO jj = 1, jpjm1
	1188	DO ji = 1, fs_jpim1 ! vector opt.
	1189	DO jk=1,jpk
	1190	pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) &
	1191	& + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) ) &
[10116]	1192	& / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) ) &
[10060]	1193	& * 0.5_wp * r1_e1e2v(ji,jj) &
	1194	* ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) &
	1195	& * ( e1e2t(ji,jj)zs(ji,jj) + e1e2t(ji,jj+1)zs(ji,jj+1) )
	1196	END DO
	1197	END DO
	1198	END DO
	1199	!
	1200	ENDIF
	1201	!
[4990]	1202	CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
[7753]	1203	pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
[5836]	1204	!
[10060]	1205	CASE( 'F' ) !* from T-point to F-point : hor. surface weighted mean
	1206	IF (nmet==0) THEN
	1207	DO jk=1,jpk
	1208	DO jj = 1, jpjm1
	1209	DO ji = 1, fs_jpim1 ! vector opt.
	1210	pe3_out(ji,jj,jk) = 0.25_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
	1211	& * r1_e1e2f(ji,jj) &
	1212	& * ( e1e2t(ji ,jj ) * ( pe3_in(ji ,jj ,jk)-e3t_0(ji ,jj ,jk) ) &
	1213	& + e1e2t(ji ,jj+1) * ( pe3_in(ji ,jj+1,jk)-e3t_0(ji ,jj+1,jk) ) &
	1214	& + e1e2t(ji+1,jj ) * ( pe3_in(ji+1,jj ,jk)-e3t_0(ji+1,jj ,jk) ) &
	1215	& + e1e2t(ji+1,jj+1) * ( pe3_in(ji+1,jj+1,jk)-e3t_0(ji+1,jj+1,jk) ) )
	1216	END DO
	1217	END DO
	1218	END DO
	1219	ELSE
[4292]	1220	DO jj = 1, jpjm1
	1221	DO ji = 1, fs_jpim1 ! vector opt.
[10060]	1222	! bottom correction:
	1223	jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1))
[10116]	1224	zdo = 0._wp
[10060]	1225	DO jk=jkbot,1,-1
[10116]	1226	zup = 0.25_wp * ( e1e2t(ji ,jj ) * zw(ji ,jj ,jk) &
	1227	& + e1e2t(ji+1,jj ) * zw(ji+1,jj ,jk) &
	1228	& + e1e2t(ji ,jj+1) * zw(ji ,jj+1,jk) &
	1229	& + e1e2t(ji+1,jj+1) * zw(ji+1,jj+1,jk) ) * r1_e1e2f(ji,jj)
[10060]	1230	!
	1231	! If there is a step, taper bottom interface:
[10116]	1232	! IF ((hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj ) + hu_0(ji,jj+1) ) ).AND.(jk==jkbot)) THEN
	1233	! IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN
	1234	! IF ( ht_0(ji+1,jj+1) < ht_0(ji ,jj+1) ) THEN
	1235	! zmin = ztap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk))
	1236	! ELSE
	1237	! zmin = ztap * (zw(ji ,jj+1,jk+1)-zw(ji ,jj+1,jk))
	1238	! ENDIF
	1239	! ELSE
	1240	! IF ( ht_0(ji+1,jj ) < ht_0(ji ,jj ) ) THEN
	1241	! zmin = ztap * (zw(ji+1,jj ,jk+1)-zw(ji+1,jj ,jk))
	1242	! ELSE
	1243	! zmin = ztap * (zw(ji ,jj ,jk+1)-zw(ji ,jj ,jk))
	1244	! ENDIF
	1245	! ENDIF
	1246	! zup = MIN(zup, zdo-zmin)
	1247	! ENDIF
	1248	zup = MIN(zup, zdo+e3f_0(ji,jj,jk)-zsmall)
[10060]	1249	!
[10116]	1250	pe3_out(ji,jj,jk) = ( zdo - zup ) &
[10060]	1251	& ( umask(ji,jj,jk) umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd )
	1252	zdo = zup
	1253	END DO
[4292]	1254	END DO
	1255	END DO
[10060]	1256	END IF
	1257	!
	1258	IF (nmet==2) THEN ! Spread sea level anomaly
	1259	!
	1260	DO jj = 1, jpjm1
	1261	DO ji = 1, fs_jpim1 ! vector opt.
	1262	DO jk=1,jpk
	1263	pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) &
	1264	& + ( pe3_out(ji,jj,jk) + e3f_0(ji,jj,jk) ) &
	1265	& / ( hf_0(ji,jj) + 1._wp - umask(ji,jj,1)*umask(ji,jj+1,1) ) &
	1266	& * 0.25_wp * r1_e1e2f(ji,jj) &
	1267	& * ( umask(ji,jj,jk)umask(ji,jj+1,jk)(1.0_wp - zlnwd) + zlnwd )&
	1268	& * ( e1e2t(ji ,jj)zs(ji ,jj) + e1e2t(ji ,jj+1)zs(ji ,jj+1) &
	1269	& +e1e2t(ji+1,jj)zs(ji+1,jj) + e1e2t(ji+1,jj+1)zs(ji+1,jj+1) )
	1270	END DO
	1271	END DO
	1272	END DO
	1273	END IF
	1274	!
[4990]	1275	CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
[7753]	1276	pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
[5836]	1277	!
	1278	CASE( 'W' ) !* from T- to W-point : vertical simple mean
	1279	!
[7753]	1280	pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
[5836]	1281	! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
	1282	!!gm BUG? use here wmask in case of ISF ? to be checked
[4292]	1283	DO jk = 2, jpk
[7753]	1284	pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
	1285	& * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) &
	1286	& + 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
	1287	& * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) )
[4292]	1288	END DO
[5836]	1289	!
	1290	CASE( 'UW' ) !* from U- to UW-point : vertical simple mean
	1291	!
[7753]	1292	pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
[4292]	1293	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
[5836]	1294	!!gm BUG? use here wumask in case of ISF ? to be checked
[4292]	1295	DO jk = 2, jpk
[7753]	1296	pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
	1297	& * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) &
	1298	& + 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
	1299	& * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) )
[4292]	1300	END DO
[5836]	1301	!
	1302	CASE( 'VW' ) !* from V- to VW-point : vertical simple mean
	1303	!
[7753]	1304	pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
[4292]	1305	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
[5836]	1306	!!gm BUG? use here wvmask in case of ISF ? to be checked
[4292]	1307	DO jk = 2, jpk
[7753]	1308	pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
	1309	& * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) &
	1310	& + 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
	1311	& * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) )
[4292]	1312	END DO
	1313	END SELECT
	1314	!
	1315	END SUBROUTINE dom_vvl_interpol
	1316
[5836]	1317
[4292]	1318	SUBROUTINE dom_vvl_rst( kt, cdrw )
	1319	!!---------------------------------------------------------------------
	1320	!! * ROUTINE dom_vvl_rst *
	1321	!!
	1322	!! ** Purpose : Read or write VVL file in restart file
	1323	!!
	1324	!! ** Method : use of IOM library
	1325	!! if the restart does not contain vertical scale factors,
	1326	!! they are set to the _0 values
	1327	!! if the restart does not contain vertical scale factors increments (z_tilde),
	1328	!! they are set to 0.
	1329	!!----------------------------------------------------------------------
	1330	INTEGER , INTENT(in) :: kt ! ocean time-step
	1331	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
[5836]	1332	!
[6152]	1333	INTEGER :: ji, jj, jk
[10164]	1334	INTEGER :: id1, id2, id3, id4, id5, id6, id7 ! local integers
[4292]	1335	!!----------------------------------------------------------------------
	1336	!
	1337	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
	1338	! ! ===============
	1339	IF( ln_rstart ) THEN !* Read the restart file
	1340	CALL rst_read_open ! open the restart file if necessary
[9367]	1341	CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn, ldxios = lrxios )
[4366]	1342	!
[6140]	1343	id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
	1344	id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
[4292]	1345	id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
	1346	id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
[10116]	1347	id5 = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. )
[10164]	1348	id6 = iom_varid( numror, 'un_lf', ldstop = .FALSE. )
	1349	id7 = iom_varid( numror, 'vn_lf', ldstop = .FALSE. )
[4292]	1350	! ! --------- !
	1351	! ! all cases !
	1352	! ! --------- !
	1353	IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist
[9367]	1354	CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
	1355	CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
[4990]	1356	! needed to restart if land processor not computed
[6140]	1357	IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
[4990]	1358	WHERE ( tmask(:,:,:) == 0.0_wp )
[6140]	1359	e3t_n(:,:,:) = e3t_0(:,:,:)
	1360	e3t_b(:,:,:) = e3t_0(:,:,:)
[4990]	1361	END WHERE
[4292]	1362	IF( neuler == 0 ) THEN
[6140]	1363	e3t_b(:,:,:) = e3t_n(:,:,:)
[4292]	1364	ENDIF
	1365	ELSE IF( id1 > 0 ) THEN
[6140]	1366	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
	1367	IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
[4990]	1368	IF(lwp) write(numout,*) 'neuler is forced to 0'
[9367]	1369	CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
[6140]	1370	e3t_n(:,:,:) = e3t_b(:,:,:)
[4990]	1371	neuler = 0
	1372	ELSE IF( id2 > 0 ) THEN
[6140]	1373	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
	1374	IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
[4490]	1375	IF(lwp) write(numout,*) 'neuler is forced to 0'
[9367]	1376	CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
[6140]	1377	e3t_b(:,:,:) = e3t_n(:,:,:)
[4490]	1378	neuler = 0
	1379	ELSE
[6140]	1380	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
[4490]	1381	IF(lwp) write(numout,*) 'Compute scale factor from sshn'
	1382	IF(lwp) write(numout,*) 'neuler is forced to 0'
[6140]	1383	DO jk = 1, jpk
	1384	e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
	1385	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
	1386	& + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
[4490]	1387	END DO
[6140]	1388	e3t_b(:,:,:) = e3t_n(:,:,:)
[4490]	1389	neuler = 0
[4292]	1390	ENDIF
	1391	! ! ----------- !
	1392	IF( ln_vvl_zstar ) THEN ! z_star case !
	1393	! ! ----------- !
	1394	IF( MIN( id3, id4 ) > 0 ) THEN
	1395	CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
	1396	ENDIF
	1397	! ! ----------------------- !
	1398	ELSE ! z_tilde and layer cases !
	1399	! ! ----------------------- !
	1400	IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist
[9367]	1401	CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )
	1402	CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )
[4292]	1403	ELSE ! one at least array is missing
	1404	tilde_e3t_b(:,:,:) = 0.0_wp
	1405	tilde_e3t_n(:,:,:) = 0.0_wp
	1406	ENDIF
	1407	! ! ------------ !
	1408	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
	1409	! ! ------------ !
[10164]	1410	IF( MIN(id5, id6, id7) > 0 ) THEN ! required arrays exist
	1411	CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:,1), ldxios = lrxios )
	1412	CALL iom_get( numror, jpdom_autoglo, 'un_lf', un_lf(:,:,:,1), ldxios = lrxios )
	1413	CALL iom_get( numror, jpdom_autoglo, 'vn_lf', vn_lf(:,:,:,1), ldxios = lrxios )
[4292]	1414	ELSE ! array is missing
[10164]	1415	hdivn_lf(:,:,:,:) = 0.0_wp
	1416	un_lf(:,:,:,:) = 0.0_wp
	1417	vn_lf(:,:,:,:) = 0.0_wp
[4292]	1418	ENDIF
	1419	ENDIF
	1420	ENDIF
	1421	!
	1422	ELSE !* Initialize at "rest"
[7646]	1423	!
[9023]	1424
	1425	IF( ll_wd ) THEN ! MJB ll_wd edits start here - these are essential
	1426	!
	1427	IF( cn_cfg == 'wad' ) THEN
	1428	! Wetting and drying test case
	1429	CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb )
	1430	tsn (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones
	1431	sshn (:,:) = sshb(:,:)
	1432	un (:,:,:) = ub (:,:,:)
	1433	vn (:,:,:) = vb (:,:,:)
	1434	ELSE
	1435	! if not test case
	1436	sshn(:,:) = -ssh_ref
	1437	sshb(:,:) = -ssh_ref
	1438
	1439	DO jj = 1, jpj
	1440	DO ji = 1, jpi
	1441	IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth
	1442
	1443	sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
	1444	sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
	1445	ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
	1446	ENDIF
	1447	ENDDO
	1448	ENDDO
	1449	ENDIF !If test case else
	1450
	1451	! Adjust vertical metrics for all wad
[7646]	1452	DO jk = 1, jpk
[9023]	1453	e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
[7646]	1454	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
[9023]	1455	& + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
[7646]	1456	END DO
	1457	e3t_b(:,:,:) = e3t_n(:,:,:)
[9023]	1458
	1459	DO ji = 1, jpi
	1460	DO jj = 1, jpj
	1461	IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
	1462	CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
	1463	ENDIF
	1464	END DO
	1465	END DO
[7646]	1466	!
	1467	ELSE
	1468	!
[9059]	1469	! Just to read set ssh in fact, called latter once vertical grid
	1470	! is set up:
[9065]	1471	! CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb )
	1472	! !
	1473	! DO jk=1,jpk
	1474	! e3t_b(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &
	1475	! & / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
	1476	! END DO
	1477	! e3t_n(:,:,:) = e3t_b(:,:,:)
	1478	sshn(:,:)=0._wp
	1479	e3t_n(:,:,:)=e3t_0(:,:,:)
	1480	e3t_b(:,:,:)=e3t_0(:,:,:)
[7646]	1481	!
[9023]	1482	END IF ! end of ll_wd edits
[6152]	1483
[4292]	1484	IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
[7646]	1485	tilde_e3t_b(:,:,:) = 0._wp
	1486	tilde_e3t_n(:,:,:) = 0._wp
[10116]	1487	IF( ln_vvl_ztilde ) THEN
[10164]	1488	hdivn_lf(:,:,:,:) = 0._wp
	1489	un_lf(:,:,:,:) = 0._wp
	1490	vn_lf(:,:,:,:) = 0._wp
[10116]	1491	ENDIF
[4292]	1492	END IF
	1493	ENDIF
[5836]	1494	!
[4292]	1495	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
	1496	! ! ===================
	1497	IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
[9367]	1498	IF( lwxios ) CALL iom_swap( cwxios_context )
[4292]	1499	! ! --------- !
	1500	! ! all cases !
	1501	! ! --------- !
[9367]	1502	CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios )
	1503	CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )
[4292]	1504	! ! ----------------------- !
	1505	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
	1506	! ! ----------------------- !
[9367]	1507	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
	1508	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
[4292]	1509	END IF
	1510	! ! -------------!
	1511	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
	1512	! ! ------------ !
[10164]	1513	CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:,1), ldxios = lwxios)
	1514	CALL iom_rstput( kt, nitrst, numrow, 'un_lf', un_lf(:,:,:,1), ldxios = lwxios)
	1515	CALL iom_rstput( kt, nitrst, numrow, 'vn_lf', vn_lf(:,:,:,1), ldxios = lwxios)
[4292]	1516	ENDIF
[5836]	1517	!
[9367]	1518	IF( lwxios ) CALL iom_swap( cxios_context )
[4292]	1519	ENDIF
[5836]	1520	!
[4292]	1521	END SUBROUTINE dom_vvl_rst
	1522
	1523
	1524	SUBROUTINE dom_vvl_ctl
	1525	!!---------------------------------------------------------------------
	1526	!! * ROUTINE dom_vvl_ctl *
	1527	!!
	1528	!! ** Purpose : Control the consistency between namelist options
	1529	!! for vertical coordinate
	1530	!!----------------------------------------------------------------------
[5836]	1531	INTEGER :: ioptio, ios
[10116]	1532
	1533	NAMELIST/nam_vvl/ ln_vvl_zstar , ln_vvl_ztilde , &
	1534	& ln_vvl_layer , ln_vvl_ztilde_as_zstar , &
	1535	& ln_vvl_zstar_at_eqtor , ln_vvl_zstar_on_shelf , &
	1536	& ln_vvl_adv_cn2 , ln_vvl_adv_fct , &
	1537	& ln_vvl_lap , ln_vvl_blp , &
	1538	& rn_ahe3_lap , rn_ahe3_blp , &
	1539	& rn_rst_e3t , rn_lf_cutoff , &
	1540	& ln_vvl_regrid , &
	1541	& ln_vvl_ramp , rn_day_ramp , &
	1542	& ln_vvl_dbg ! not yet implemented: ln_vvl_kepe
	1543	!!----------------------------------------------------------------------
[5836]	1544	!
[4294]	1545	REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist :
	1546	READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
[9168]	1547	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
[4294]	1548	REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run
	1549	READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
[9168]	1550	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
[4624]	1551	IF(lwm) WRITE ( numond, nam_vvl )
[5836]	1552	!
[4292]	1553	IF(lwp) THEN ! Namelist print
	1554	WRITE(numout,*)
	1555	WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
	1556	WRITE(numout,*) '~~~~~~~~~~~'
[10116]	1557	WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate'
	1558	WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar
	1559	WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde
	1560	WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer
	1561	WRITE(numout,*) ' ztilde as zstar ln_vvl_ztilde_as_zstar = ', ln_vvl_ztilde_as_zstar
	1562	! WRITE(numout,*) ' Namelist nam_vvl : chose kinetic-to-potential energy conservation'
	1563	! WRITE(numout,*) ' ln_vvl_kepe = ', ln_vvl_kepe
[4292]	1564	WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor
[10116]	1565	WRITE(numout,*) ' ztilde on shelves ln_vvl_zstar_on_shelf = ', ln_vvl_zstar_on_shelf
	1566	WRITE(numout,*) ' Namelist nam_vvl : thickness advection scheme'
	1567	WRITE(numout,*) ' 2nd order ln_vvl_adv_cn2 = ', ln_vvl_adv_cn2
	1568	WRITE(numout,*) ' 2nd order FCT ln_vvl_adv_fct = ', ln_vvl_adv_fct
	1569	WRITE(numout,*) ' Namelist nam_vvl : thickness diffusion scheme'
	1570	WRITE(numout,*) ' Laplacian ln_vvl_lap = ', ln_vvl_lap
	1571	WRITE(numout,*) ' Bilaplacian ln_vvl_blp = ', ln_vvl_blp
	1572	WRITE(numout,*) ' Laplacian coefficient rn_ahe3_lap = ', rn_ahe3_lap
	1573	WRITE(numout,*) ' Bilaplacian coefficient rn_ahe3_blp = ', rn_ahe3_blp
	1574	WRITE(numout,*) ' Namelist nam_vvl : layers regriding'
	1575	WRITE(numout,*) ' ln_vvl_regrid = ', ln_vvl_regrid
	1576	WRITE(numout,*) ' Namelist nam_vvl : linear ramp at startup'
	1577	WRITE(numout,*) ' ln_vvl_ramp = ', ln_vvl_ramp
	1578	WRITE(numout,*) ' rn_day_ramp = ', rn_day_ramp
[4292]	1579	IF( ln_vvl_ztilde_as_zstar ) THEN
[10116]	1580	WRITE(numout,*) ' ztilde running in zstar emulation mode; '
	1581	WRITE(numout,*) ' ignoring namelist timescale parameters and using:'
	1582	WRITE(numout,*) ' hard-wired : z-tilde to zstar restoration timescale (days)'
	1583	WRITE(numout,*) ' rn_rst_e3t = 0.0'
	1584	WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
	1585	WRITE(numout,*) ' rn_lf_cutoff = 1.0/rdt'
[4292]	1586	ELSE
[10116]	1587	WRITE(numout,*) ' Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
	1588	WRITE(numout,*) ' rn_rst_e3t = ', rn_rst_e3t
	1589	WRITE(numout,*) ' Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
	1590	WRITE(numout,*) ' rn_lf_cutoff = ', rn_lf_cutoff
[4292]	1591	ENDIF
[10116]	1592	WRITE(numout,*) ' Namelist nam_vvl : debug prints'
	1593	WRITE(numout,*) ' ln_vvl_dbg = ', ln_vvl_dbg
[4292]	1594	ENDIF
[5836]	1595	!
[10116]	1596	IF ( ln_vvl_ztilde.OR.ln_vvl_layer ) THEN
	1597	ioptio = 0 ! Choose one advection scheme at most
	1598	IF( ln_vvl_adv_cn2 ) ioptio = ioptio + 1
	1599	IF( ln_vvl_adv_fct ) ioptio = ioptio + 1
	1600	IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE thickness advection scheme in namelist nam_vvl' )
	1601	ENDIF
	1602	!
[4292]	1603	ioptio = 0 ! Parameter control
[5836]	1604	IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
	1605	IF( ln_vvl_zstar ) ioptio = ioptio + 1
	1606	IF( ln_vvl_ztilde ) ioptio = ioptio + 1
	1607	IF( ln_vvl_layer ) ioptio = ioptio + 1
	1608	!
[4292]	1609	IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
[6140]	1610	IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
[5836]	1611	!
[4292]	1612	IF(lwp) THEN ! Print the choice
	1613	WRITE(numout,*)
[9168]	1614	IF( ln_vvl_zstar ) WRITE(numout,*) ' ==>>> zstar vertical coordinate is used'
	1615	IF( ln_vvl_ztilde ) WRITE(numout,*) ' ==>>> ztilde vertical coordinate is used'
	1616	IF( ln_vvl_layer ) WRITE(numout,*) ' ==>>> layer vertical coordinate is used'
	1617	IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) ' ==>>> to emulate a zstar coordinate'
[4292]	1618	ENDIF
[5836]	1619	!
[10116]	1620	! Use of "shelf horizon depths" should be allowed with s-z coordinates, but we restrict it to zco and zps
	1621	! for the time being
	1622	IF ( ln_sco ) THEN
	1623	ll_shorizd=.FALSE.
	1624	ELSE
	1625	ll_shorizd=.TRUE.
	1626	ENDIF
	1627	!
[4486]	1628	#if defined key_agrif
[9190]	1629	IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) ) CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
[4486]	1630	#endif
[5836]	1631	!
[4292]	1632	END SUBROUTINE dom_vvl_ctl
	1633
[10116]	1634	SUBROUTINE dom_vvl_regrid( kt )
	1635	!!----------------------------------------------------------------------
	1636	!! * ROUTINE dom_vvl_regrid *
	1637	!!
	1638	!! ** Purpose : Ensure "well-behaved" vertical grid
	1639	!!
	1640	!! ** Method : More or less adapted from references below.
	1641	!!regrid
	1642	!! ** Action : Ensure that thickness are above a given value, spaced enough
	1643	!! and revert to Eulerian coordinates near the bottom.
	1644	!!
	1645	!! References : Bleck, R. and S. Benjamin, 1993: Regional Weather Prediction
	1646	!! with a Model Combining Terrain-following and Isentropic
	1647	!! coordinates. Part I: Model Description. Monthly Weather Rev.,
	1648	!! 121, 1770-1785.
	1649	!! Toy, M., 2011: Incorporating Condensational Heating into a
	1650	!! Nonhydrostatic Atmospheric Model Based on a Hybrid Isentropic-
	1651	!! Sigma Vertical Coordinate. Monthly Weather Rev., 139, 2940-2954.
	1652	!!----------------------------------------------------------------------
	1653	!! * Arguments
	1654	INTEGER, INTENT( in ) :: kt ! time step
	1655
	1656	!! * Local declarations
	1657	INTEGER :: ji, jj, jk ! dummy loop indices
	1658	LOGICAL :: ll_chk_bot2top, ll_chk_top2bot, ll_lapdiff_cond
	1659	LOGICAL :: ll_zdiff_cond, ll_blpdiff_cond
	1660	INTEGER :: jkbot
	1661	REAL(wp) :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd
	1662	REAL(wp) :: zufim1, zufi, zvfjm1, zvfj, dzmin_int, dzmin_surf
	1663	REAL(wp) :: zh_new, zh_old, zh_bef, ztmp, ztmp1, z2dt, zh_up, zh_dwn
	1664	REAL(wp) :: zeu2, zev2, zfrch_stp, zfrch_rel, zfrac_bot, zscal_bot
	1665	REAL(wp) :: zhdiff, zhdiff2, zvdiff, zhlim, zhlim2, zvlim
	1666	REAL(wp), DIMENSION(jpi,jpj) :: zdw, zwu, zwv
	1667	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwdw, zwdw_b
	1668	!!----------------------------------------------------------------------
	1669
	1670	IF( ln_timing ) CALL timing_start('dom_vvl_regrid')
	1671	!
	1672	!
	1673	! Some user defined parameters below:
	1674	ll_chk_bot2top = .TRUE.
	1675	ll_chk_top2bot = .TRUE.
[10164]	1676	dzmin_int = 1.0_wp ! Absolute minimum depth in the interior (in meters)
[10116]	1677	dzmin_surf = 1.0_wp ! Absolute minimum depth at the surface (in meters)
[10164]	1678	zfrch_stp = 5._wp ! Maximum fractionnal thickness change in one time step (<= 1.)
[10116]	1679	zfrch_rel = 0.4_wp ! Maximum relative thickness change in the vertical (<= 1.)
	1680	zfrac_bot = 0.05_wp ! Fraction of bottom level allowed to change
	1681	zscal_bot = 2.0_wp ! Depth lengthscale
	1682	ll_zdiff_cond = .TRUE. ! Conditionnal vertical diffusion of interfaces
	1683	zvdiff = 0.2_wp ! m
	1684	zvlim = 0.5_wp ! max d2h/dh
	1685	ll_lapdiff_cond = .TRUE. ! Conditionnal Laplacian diffusion of interfaces
[10164]	1686	zhdiff = 0.01_wp ! ad.
[10116]	1687	zhlim = 0.03_wp ! ad. max lap(z)*e1
[10164]	1688	ll_blpdiff_cond = .TRUE. ! Conditionnal Bilaplacian diffusion of interfaces
[10116]	1689	zhdiff2 = 0.2_wp ! ad.
	1690	zhlim2 = 0.01_wp ! ad. max bilap(z)e1*3
	1691	! ---------------------------------------------------------------------------------------
	1692	!
	1693	! Set arrays determining maximum vertical displacement at the bottom:
	1694	!--------------------------------------------------------------------
	1695	IF ( kt==nit000 ) THEN
	1696	DO jj = 2, jpjm1
	1697	DO ji = 2, jpim1
	1698	jk = MIN(mbkt(ji,jj), mbkt(ji+1,jj), mbkt(ji-1,jj), mbkt(ji,jj+1), mbkt(ji,jj-1))
	1699	jk = MIN(jk,mbkt(ji-1,jj-1), mbkt(ji-1,jj+1), mbkt(ji+1,jj+1), mbkt(ji+1,jj-1))
	1700	i_int_bot(ji,jj) = jk
	1701	END DO
	1702	END DO
	1703	dsm(:,:) = REAL( i_int_bot(:,:), wp ) ; CALL lbc_lnk(dsm(:,:),'T',1.)
	1704	i_int_bot(:,:) = MAX( INT( dsm(:,:) ), 1 )
	1705
	1706	DO jj = 2, jpjm1
	1707	DO ji = 2, jpim1
	1708	zdw(ji,jj) = MAX(ABS(ht_0(ji,jj)-ht_0(ji+1,jj))*umask(ji ,jj,1), &
	1709	& ABS(ht_0(ji,jj)-ht_0(ji-1,jj))*umask(ji-1,jj,1), &
	1710	& ABS(ht_0(ji,jj)-ht_0(ji,jj+1))*vmask(ji,jj ,1), &
	1711	& ABS(ht_0(ji,jj)-ht_0(ji,jj-1))*vmask(ji,jj-1,1) )
	1712	zdw(ji,jj) = MAX(zscal_bot * zdw(ji,jj), rsmall )
	1713	END DO
	1714	END DO
	1715	CALL lbc_lnk( zdw(:,:), 'T', 1. )
	1716
	1717	DO jj = 2, jpjm1
	1718	DO ji = 2, jpim1
	1719	dsm(ji,jj) = 1._wp/16._wp * ( zdw(ji-1,jj-1) + zdw(ji+1,jj-1) &
	1720	& + zdw(ji-1,jj+1) + zdw(ji+1,jj+1) &
	1721	& + 2._wp*( zdw(ji ,jj-1) + zdw(ji-1,jj ) &
	1722	& + zdw(ji+1,jj ) + zdw(ji ,jj+1) ) &
	1723	& + 4._wp* zdw(ji ,jj ) )
	1724	END DO
	1725	END DO
	1726
	1727	CALL lbc_lnk( dsm(:,:), 'T', 1. )
	1728
	1729	IF (ln_zps) THEN
	1730	DO jj = 1, jpj
	1731	DO ji = 1, jpi
	1732	jk = i_int_bot(ji,jj)
	1733	hsm(ji,jj) = zfrac_bot * e3w_1d(jk)
[10128]	1734	dsm(ji,jj) = MAX(dsm(ji,jj), 0.05_wp*ht_0(ji,jj))
[10116]	1735	END DO
	1736	END DO
	1737	ELSE
	1738	DO jj = 1, jpj
	1739	DO ji = 1, jpi
	1740	jk = i_int_bot(ji,jj)
	1741	hsm(ji,jj) = zfrac_bot * e3w_0(ji,jj,jk)
[10128]	1742	dsm(ji,jj) = MAX(dsm(ji,jj), 0.05_wp*ht_0(ji,jj))
[10116]	1743	END DO
	1744	END DO
	1745	ENDIF
	1746	END IF
	1747
	1748	! Provisionnal interface depths:
	1749	!-------------------------------
	1750	zwdw(:,:,1) = 0.e0
	1751	DO jj = 1, jpj
	1752	DO ji = 1, jpi
	1753	DO jk = 2, jpk
	1754	zwdw(ji,jj,jk) = zwdw(ji,jj,jk-1) + &
	1755	& (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
	1756	END DO
	1757	END DO
	1758	END DO
	1759	!
	1760	! Conditionnal horizontal Laplacian diffusion:
	1761	!---------------------------------------------
	1762	IF ( ll_lapdiff_cond ) THEN
	1763	!
	1764	zwdw_b(:,:,1) = 0._wp
	1765	DO jj = 1, jpj
	1766	DO ji = 1, jpi
	1767	DO jk=2,jpk
	1768	zwdw_b(ji,jj,jk) = zwdw_b(ji,jj,jk-1) + &
	1769	& (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
	1770	END DO
	1771	END DO
	1772	END DO
	1773	!
	1774	DO jk = 2, jpkm1
	1775	zwu(:,:) = 0._wp
	1776	zwv(:,:) = 0._wp
	1777	DO jj = 1, jpjm1
	1778	DO ji = 1, fs_jpim1 ! vector opt.
	1779	zwu(ji,jj) = umask(ji,jj,jk) * e2_e1u(ji,jj) &
	1780	& * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj ,jk) )
	1781	zwv(ji,jj) = vmask(ji,jj,jk) * e1_e2v(ji,jj) &
	1782	& * ( zwdw_b(ji,jj,jk) - zwdw_b(ji ,jj+1,jk) )
	1783	END DO
	1784	END DO
	1785	DO jj = 2, jpjm1
	1786	DO ji = fs_2, fs_jpim1 ! vector opt.
	1787	ztmp1 = ( zwu(ji-1,jj ) - zwu(ji,jj) &
	1788	& + zwv(ji ,jj-1) - zwv(ji,jj) ) * r1_e1e2t(ji,jj)
	1789	zh2 = MAX(abs(ztmp1)-zhlim*SQRT(r1_e1e2t(ji,jj)), 0._wp)
	1790	ztmp = SIGN(zh2, ztmp1)
	1791	zeu2 = zhdiff * e1e2t(ji,jj)e1e2t(ji,jj)/(e1t(ji,jj)e1t(ji,jj) + e2t(ji,jj)*e2t(ji,jj))
	1792	zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp * tmask(ji,jj,jk)
	1793	END DO
	1794	END DO
	1795	END DO
	1796	!
	1797	ENDIF
	1798
	1799	! Conditionnal horizontal Bilaplacian diffusion:
	1800	!-----------------------------------------------
	1801	IF ( ll_blpdiff_cond ) THEN
	1802	!
	1803	zwdw_b(:,:,1) = 0._wp
	1804	DO jj = 1, jpj
	1805	DO ji = 1, jpi
	1806	DO jk = 2,jpkm1
	1807	zwdw_b(ji,jj,jk) = zwdw_b(ji,jj,jk-1) + &
	1808	& (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
	1809	END DO
	1810	END DO
	1811	END DO
	1812	!
	1813	DO jk = 2, jpkm1
	1814	zwu(:,:) = 0._wp
	1815	zwv(:,:) = 0._wp
	1816	DO jj = 1, jpjm1
	1817	DO ji = 1, fs_jpim1 ! vector opt.
	1818	zwu(ji,jj) = umask(ji,jj,jk) * e2_e1u(ji,jj) &
	1819	& * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj ,jk) )
	1820	zwv(ji,jj) = vmask(ji,jj,jk) * e1_e2v(ji,jj) &
	1821	& * ( zwdw_b(ji,jj,jk) - zwdw_b(ji ,jj+1,jk) )
	1822	END DO
	1823	END DO
	1824	DO jj = 2, jpjm1
	1825	DO ji = fs_2, fs_jpim1 ! vector opt.
	1826	zwdw_b(ji,jj,jk) = -( (zwu(ji-1,jj ) - zwu(ji,jj)) &
	1827	& + (zwv(ji ,jj-1) - zwv(ji,jj)) ) * r1_e1e2t(ji,jj)
	1828	END DO
	1829	END DO
	1830	END DO
	1831	!
	1832	CALL lbc_lnk( zwdw_b(:,:,:), 'T', 1. )
	1833	!
	1834	DO jk = 2, jpkm1
	1835	zwu(:,:) = 0._wp
	1836	zwv(:,:) = 0._wp
	1837	DO jj = 1, jpjm1
	1838	DO ji = 1, fs_jpim1 ! vector opt.
	1839	zwu(ji,jj) = umask(ji,jj,jk) * e2_e1u(ji,jj) &
	1840	& * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj ,jk) )
	1841	zwv(ji,jj) = vmask(ji,jj,jk) * e1_e2v(ji,jj) &
	1842	& * ( zwdw_b(ji,jj,jk) - zwdw_b(ji ,jj+1,jk) )
	1843	END DO
	1844	END DO
	1845	DO jj = 2, jpjm1
	1846	DO ji = fs_2, fs_jpim1 ! vector opt.
	1847	ztmp1 = ( (zwu(ji-1,jj ) - zwv(ji,jj)) &
	1848	& + (zwv(ji ,jj-1) - zwv(ji,jj)) ) * r1_e1e2t(ji,jj)
	1849	zh2 = MAX(abs(ztmp1)-zhlim2SQRT(r1_e1e2t(ji,jj))r1_e1e2t(ji,jj), 0._wp)
	1850	ztmp = SIGN(zh2, ztmp1)
	1851	zeu2 = zhdiff2 * e1e2t(ji,jj)*e1e2t(ji,jj) / 16._wp
	1852	zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp * tmask(ji,jj,jk)
	1853	END DO
	1854	END DO
	1855	END DO
	1856	!
	1857	ENDIF
	1858
	1859	! Conditionnal vertical diffusion:
	1860	!---------------------------------
	1861	IF ( ll_zdiff_cond ) THEN
	1862	DO jk = 2, jpkm1
	1863	DO jj = 2, jpjm1
	1864	DO ji = fs_2, fs_jpim1 ! vector opt.
	1865	ztmp = -( (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))*tmask(ji,jj,jk-1) &
	1866	-(tilde_e3t_b(ji,jj,jk )+e3t_0(ji,jj,jk ))*tmask(ji,jj,jk ) )
	1867	ztmp1 = 0.5_wp * ( tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) &
	1868	& +tilde_e3t_b(ji,jj,jk ) + e3t_0(ji,jj,jk ) )
	1869	zh2 = MAX(abs(ztmp)-zvlim*ztmp1, 0._wp)
	1870	ztmp = SIGN(zh2, ztmp)
	1871	IF ((jk==mbkt(ji,jj)).AND.(ln_zps)) ztmp=0.e0
	1872	zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zvdiff * ztmp * tmask(ji,jj,jk)
	1873	END DO
	1874	END DO
	1875	END DO
	1876	ENDIF
	1877	!
	1878	! Check grid from the bottom to the surface
	1879	!------------------------------------------
	1880	IF ( ll_chk_bot2top ) THEN
	1881	DO jj = 2, jpjm1
	1882	DO ji = 2, jpim1
	1883	jkbot = mbkt(ji,jj)
	1884	DO jk = jkbot,2,-1
	1885	!
	1886	zh_0 = e3t_0(ji,jj,jk)
	1887	zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1))
	1888	zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
	1889	zh_min = MIN(zh_0/3._wp, dzmin_int)
[10164]	1890	! zh_min = MAX(zh_min, zh_min-e3t_a(ji,jj,jk)+e3t_0(ji,jj,jk))
[10116]	1891	!
	1892	! Set maximum and minimum vertical excursions
	1893	ztmph = hsm(ji,jj)
	1894	ztmpd = dsm(ji,jj)
	1895	zh2 = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)))/ztmpd)
	1896	! zh2 = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd)
	1897	zh2 = MAX(zh2,0.001_wp) ! Extend tolerance a bit for stability reasons (to be explored)
	1898	zdiff = cush_max(gdepw_0(ji,jj,jk)-zwdw(ji,jj,jk), zh2 )
	1899	zwdw(ji,jj,jk) = MAX(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) - zdiff)
	1900	zdiff = cush_max(zwdw(ji,jj,jk)-gdepw_0(ji,jj,jk), zh2 )
	1901	zwdw(ji,jj,jk) = MIN(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) + zdiff)
	1902	!
	1903	! New layer thickness:
	1904	zh_new = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
	1905	!
	1906	! Ensure minimum layer thickness:
	1907	! zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
[10164]	1908	! zh_new = cush(zh_new, zh_min)
	1909	zh_new = MAX(zh_new, zh_min)
[10116]	1910	!
	1911	! Final flux:
	1912	zdiff = (zh_new - zh_old) * tmask(ji,jj,jk)
	1913	!
	1914	! Limit thickness change in 1 time step:
[10164]	1915	! ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
	1916	! zdiff = SIGN(ztmp, zh_new - zh_old)
[10116]	1917	zh_new = zdiff + zh_old
	1918	!
	1919	zwdw(ji,jj,jk) = zwdw(ji,jj,jk+1) - zh_new
	1920	END DO
	1921	END DO
	1922	END DO
	1923	END IF
	1924	!
	1925	! Check grid from the surface to the bottom
	1926	!------------------------------------------
	1927	IF ( ll_chk_top2bot ) THEN
	1928	DO jj = 2, jpjm1
	1929	DO ji = 2, jpim1
	1930	jkbot = mbkt(ji,jj)
	1931	DO jk = 1, jkbot-1
	1932	!
	1933	zh_0 = e3t_0(ji,jj,jk)
	1934	zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk+1) + e3t_0(ji,jj,jk+1))
	1935	zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
	1936	zh_min = MIN(zh_0/3._wp, dzmin_int)
[10164]	1937	! zh_min = MAX(zh_min, zh_min-e3t_a(ji,jj,jk)+e3t_0(ji,jj,jk))
[10116]	1938	!
[10164]	1939	! zwdw(ji,jj,jk+1) = MAX(zwdw(ji,jj,jk+1), REAL(jk)*dzmin_surf)
[10116]	1940	!
	1941	! New layer thickness:
	1942	zh_new = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
	1943	!
	1944	! Ensure minimum layer thickness:
	1945	! zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
[10164]	1946	! zh_new = cush(zh_new, zh_min)
	1947	zh_new = MAX(zh_new, zh_min)
[10116]	1948	!
	1949	! Final flux:
	1950	zdiff = (zh_new -zh_old) * tmask(ji,jj,jk)
	1951	!
	1952	! Limit flux:
	1953	! ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
	1954	! zdiff = SIGN(ztmp, zh_new - zh_old)
	1955	zh_new = zdiff + zh_old
	1956	!
	1957	zwdw(ji,jj,jk+1) = zwdw(ji,jj,jk) + zh_new
	1958	END DO
	1959	!
	1960	END DO
	1961	END DO
	1962	ENDIF
	1963	!
	1964	DO jj = 2, jpjm1
	1965	DO ji = 2, jpim1
	1966	DO jk = 1, jpkm1
	1967	tilde_e3t_a(ji,jj,jk) = (zwdw(ji,jj,jk+1)-zwdw(ji,jj,jk)-e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
	1968	END DO
	1969	END DO
	1970	END DO
	1971	!
	1972	!
	1973	IF( ln_timing ) CALL timing_stop('dom_vvl_regrid')
	1974	!
	1975	END SUBROUTINE dom_vvl_regrid
	1976
	1977	FUNCTION cush(hin, hmin) RESULT(hout)
	1978	!!----------------------------------------------------------------------
	1979	!! * FUNCTION cush *
	1980	!!
	1981	!! ** Purpose :
	1982	!!
	1983	!! ** Method :
	1984	!!
	1985	!!----------------------------------------------------------------------
	1986	IMPLICIT NONE
	1987	REAL(wp), INTENT(in) :: hin, hmin
	1988	REAL(wp) :: hout, zx, zh_cri
	1989	!!----------------------------------------------------------------------
	1990	zh_cri = 3._wp * hmin
	1991	!
	1992	IF ( hin<=0._wp ) THEN
	1993	hout = hmin
	1994	!
	1995	ELSEIF ( (hin>0._wp).AND.(hin<=zh_cri) ) THEN
	1996	zx = hin/zh_cri
	1997	hout = hmin * (1._wp + zx + zx*zx)
	1998	!
	1999	ELSEIF ( hin>zh_cri ) THEN
	2000	hout = hin
	2001	!
	2002	ENDIF
	2003	!
	2004	END FUNCTION cush
	2005
	2006	FUNCTION cush_max(hin, hmax) RESULT(hout)
	2007	!!----------------------------------------------------------------------
	2008	!! * FUNCTION cush *
	2009	!!
	2010	!! ** Purpose :
	2011	!!
	2012	!! ** Method :
	2013	!!
	2014	!!----------------------------------------------------------------------
	2015	IMPLICIT NONE
	2016	REAL(wp), INTENT(in) :: hin, hmax
	2017	REAL(wp) :: hout, hmin, zx, zh_cri
	2018	!!----------------------------------------------------------------------
	2019	hmin = 0.1_wp * hmax
	2020	zh_cri = 3._wp * hmin
	2021	!
	2022	IF ( (hin>=(hmax-zh_cri)).AND.(hin<=(hmax-hmin))) THEN
	2023	zx = (hmax-hin)/zh_cri
	2024	hout = hmax - hmin * (1._wp + zx + zx*zx)
	2025	!
	2026	ELSEIF ( hin>(hmax-zh_cri) ) THEN
	2027	hout = hmax - hmin
	2028	!
	2029	ELSE
	2030	hout = hin
	2031	!
	2032	ENDIF
	2033	!
	2034	END FUNCTION cush_max
	2035
	2036	SUBROUTINE dom_vvl_adv_fct( kt, pta, uin, vin )
	2037	!!----------------------------------------------------------------------
	2038	!! * ROUTINE dom_vvl_adv_fct *
	2039	!!
	2040	!! ** Purpose : Do thickness advection
	2041	!!
	2042	!! ** Method : FCT scheme to ensure positivity
	2043	!!
	2044	!! ** Action : - Update pta thickness tendency and diffusive fluxes
	2045	!! - this is the total trend, hence it does include sea level motions
	2046	!!----------------------------------------------------------------------
	2047	!
	2048	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	2049	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pta ! thickness baroclinic trend
	2050	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: uin, vin ! input velocities
	2051	!
	2052	INTEGER :: ji, jj, jk, ib, ib_bdy ! dummy loop indices
	2053	INTEGER :: ikbu, ikbv, ibot
	2054	REAL(wp) :: z2dtt, zbtr, ztra ! local scalar
	2055	REAL(wp) :: zdi, zdj, zmin ! - -
	2056	REAL(wp) :: zfp_ui, zfp_vj ! - -
	2057	REAL(wp) :: zfm_ui, zfm_vj ! - -
	2058	REAL(wp) :: zfp_hi, zfp_hj ! - -
	2059	REAL(wp) :: zfm_hi, zfm_hj ! - -
	2060	REAL(wp) :: ztout, ztin, zfac ! - -
	2061	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy, zwi
	2062	!!----------------------------------------------------------------------
	2063	!
	2064	IF( ln_timing ) CALL timing_start('dom_vvl_adv_fct')
	2065	!
	2066	!
	2067	! 1. Initializations
	2068	! ------------------
	2069	!
	2070	IF( neuler == 0 .AND. kt == nit000 ) THEN
	2071	z2dtt = rdt
	2072	ELSE
	2073	z2dtt = 2.0_wp * rdt
	2074	ENDIF
	2075	!
	2076	zwi(:,:,:) = 0.e0
	2077	zwx(:,:,:) = 0.e0
	2078	zwy(:,:,:) = 0.e0
	2079	!
	2080	!
	2081	! 2. upstream advection with initial mass fluxes & intermediate update
	2082	! --------------------------------------------------------------------
	2083	IF ( ll_shorizd ) THEN
	2084	DO jk = 1, jpkm1
	2085	DO jj = 1, jpjm1
	2086	DO ji = 1, fs_jpim1 ! vector opt.
	2087	!
	2088	zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji ,jj ,jk), 0._wp)
	2089	zfp_hi = MIN(zfp_hi, e3t_b(ji ,jj ,jk))
	2090	zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) )
	2091	!
	2092	zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj ,jk), 0._wp)
	2093	zfm_hi = MIN(zfm_hi, e3t_b(ji+1,jj ,jk))
	2094	zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) )
	2095	!
	2096	zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji ,jj ,jk), 0._wp)
	2097	zfp_hj = MIN(zfp_hj, e3t_b(ji ,jj ,jk))
	2098	zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) )
	2099	!
	2100	zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji ,jj+1,jk), 0._wp)
	2101	zfm_hj = MIN(zfm_hj, e3t_b(ji ,jj+1,jk))
	2102	zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) )
	2103	!
	2104	zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) )
	2105	zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) )
	2106	zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) )
	2107	zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) )
	2108	zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
	2109	zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
	2110	END DO
	2111	END DO
	2112	END DO
	2113	ELSE
	2114	DO jk = 1, jpkm1
	2115	DO jj = 1, jpjm1
	2116	DO ji = 1, fs_jpim1 ! vector opt.
	2117	!
	2118	zfp_hi = e3t_b(ji ,jj ,jk)
	2119	zfm_hi = e3t_b(ji+1,jj ,jk)
	2120	zfp_hj = e3t_b(ji ,jj ,jk)
	2121	zfm_hj = e3t_b(ji ,jj+1,jk)
	2122	!
	2123	zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) )
	2124	zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) )
	2125	zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) )
	2126	zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) )
	2127	zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
	2128	zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
	2129	END DO
	2130	END DO
	2131	END DO
	2132	ENDIF
	2133
	2134	! total advective trend
	2135	DO jk = 1, jpkm1
	2136	DO jj = 2, jpjm1
	2137	DO ji = fs_2, fs_jpim1 ! vector opt.
	2138	zbtr = r1_e1e2t(ji,jj)
	2139	! total intermediate advective trends
	2140	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	2141	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) )
	2142	!
	2143	! update and guess with monotonic sheme
	2144	pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
	2145	zwi(ji,jj,jk) = (e3t_b(ji,jj,jk) + z2dtt * ztra ) * tmask(ji,jj,jk)
	2146	END DO
	2147	END DO
	2148	END DO
	2149
	2150	CALL lbc_lnk( zwi, 'T', 1. )
	2151
	2152	IF ( ln_bdy ) THEN
	2153	DO ib_bdy=1, nb_bdy
	2154	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
	2155	ji = idx_bdy(ib_bdy)%nbi(ib,1)
	2156	jj = idx_bdy(ib_bdy)%nbj(ib,1)
	2157	DO jk = 1, jpkm1
	2158	zwi(ji,jj,jk) = e3t_a(ji,jj,jk)
	2159	END DO
	2160	END DO
	2161	END DO
	2162	ENDIF
	2163
[10164]	2164	! IF ( ln_vvl_dbg ) THEN
	2165	! zmin = MINVAL( zwi(:,:,:), mask = tmask(:,:,:) == 1.e0 )
	2166	! IF( lk_mpp ) CALL mpp_min( zmin )
	2167	! IF( zmin < 0._wp) THEN
	2168	! IF(lwp) CALL ctl_warn('vvl_adv: CFL issue here')
	2169	! IF(lwp) WRITE(numout,*) zmin
	2170	! ENDIF
	2171	! ENDIF
[10116]	2172
	2173	! 3. antidiffusive flux : high order minus low order
	2174	! --------------------------------------------------
	2175	! antidiffusive flux on i and j
	2176	DO jk = 1, jpkm1
	2177	DO jj = 1, jpjm1
	2178	DO ji = 1, fs_jpim1 ! vector opt.
	2179	zwx(ji,jj,jk) = (e2u(ji,jj) * uin(ji,jj,jk) * e3u_n(ji,jj,jk) &
	2180	& - zwx(ji,jj,jk)) * umask(ji,jj,jk)
	2181	zwy(ji,jj,jk) = (e1v(ji,jj) * vin(ji,jj,jk) * e3v_n(ji,jj,jk) &
	2182	& - zwy(ji,jj,jk)) * vmask(ji,jj,jk)
	2183	!
	2184	! Update advective fluxes
	2185	un_td(ji,jj,jk) = un_td(ji,jj,jk) - zwx(ji,jj,jk)
	2186	vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - zwy(ji,jj,jk)
	2187	END DO
	2188	END DO
	2189	END DO
	2190
[10128]	2191	CALL lbc_lnk_multi( zwx, 'U', -1., zwy, 'V', -1. ) !* local domain boundaries
[10116]	2192
	2193	! 4. monotonicity algorithm
	2194	! -------------------------
	2195	CALL nonosc_2d( e3t_b(:,:,:), zwx, zwy, zwi, z2dtt )
	2196
	2197	! 5. final trend with corrected fluxes
	2198	! ------------------------------------
	2199	!
	2200	! Update advective fluxes
	2201	un_td(:,:,:) = (un_td(:,:,:) + zwx(:,:,:))*umask(:,:,:)
	2202	vn_td(:,:,:) = (vn_td(:,:,:) + zwy(:,:,:))*vmask(:,:,:)
	2203	!
	2204	DO jk = 1, jpkm1
	2205	DO jj = 2, jpjm1
	2206	DO ji = fs_2, fs_jpim1 ! vector opt.
	2207	!
	2208	zbtr = r1_e1e2t(ji,jj)
	2209	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	2210	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) )
	2211	! add them to the general tracer trends
	2212	pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
[10164]	2213	zwi(ji,jj,jk) = (e3t_b(ji,jj,jk) + z2dtt * pta(ji,jj,jk)* bdytmask(ji,jj) ) * tmask(ji,jj,jk)
[10116]	2214	END DO
	2215	END DO
	2216	END DO
[10164]	2217
	2218	IF ( ln_vvl_dbg ) THEN
	2219	zmin = MINVAL( zwi(:,:,:), mask = tmask(:,:,:) == 1.e0 )
	2220	IF( lk_mpp ) CALL mpp_min( zmin )
	2221	IF( zmin < 0._wp) THEN
	2222	IF(lwp) CALL ctl_warn('vvl_adv: CFL issue here')
	2223	IF(lwp) WRITE(numout,*) zmin
	2224	ENDIF
	2225	ENDIF
[10116]	2226	!
	2227	IF( ln_timing ) CALL timing_stop('dom_vvl_adv_fct')
	2228	!
	2229	END SUBROUTINE dom_vvl_adv_fct
	2230
	2231	SUBROUTINE dom_vvl_ups_cor( kt, pta, uin, vin )
	2232	!!----------------------------------------------------------------------
	2233	!! * ROUTINE dom_vvl_adv_fct *
	2234	!!
	2235	!! ** Purpose : Correct for addionnal barotropic fluxes
	2236	!! in the upstream direction
	2237	!!
	2238	!! ** Method :
	2239	!!
	2240	!! ** Action : - Update diffusive fluxes uin, vin
	2241	!! - Remove divergence from thickness tendency
	2242	!!----------------------------------------------------------------------
	2243	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	2244	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pta ! thickness baroclinic trend
	2245	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: uin, vin ! input fluxes
	2246	INTEGER :: ji, jj, jk ! dummy loop indices
	2247	INTEGER :: ikbu, ikbv, ibot
	2248	REAL(wp) :: zbtr, ztra ! local scalar
	2249	REAL(wp) :: zdi, zdj ! - -
	2250	REAL(wp) :: zfp_hi, zfp_hj ! - -
	2251	REAL(wp) :: zfm_hi, zfm_hj ! - -
	2252	REAL(wp) :: zfp_ui, zfp_vj ! - -
	2253	REAL(wp) :: zfm_ui, zfm_vj ! - -
	2254	REAL(wp), DIMENSION(jpi,jpj) :: zbu, zbv, zhu_b, zhv_b
	2255	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy
	2256	!!----------------------------------------------------------------------
	2257	!
	2258	IF( ln_timing ) CALL timing_start('dom_vvl_ups_cor')
	2259	!
	2260	! Compute barotropic flux difference:
	2261	zbu(:,:) = 0.e0
	2262	zbv(:,:) = 0.e0
	2263	DO jj = 1, jpj
	2264	DO ji = 1, jpi ! vector opt.
	2265	DO jk = 1, jpkm1
	2266	zbu(ji,jj) = zbu(ji,jj) - uin(ji,jj,jk) * umask(ji,jj,jk)
	2267	zbv(ji,jj) = zbv(ji,jj) - vin(ji,jj,jk) * vmask(ji,jj,jk)
	2268	END DO
	2269	END DO
	2270	ENDDO
[10128]	2271
[10116]	2272	! Compute upstream depths:
	2273	zhu_b(:,:) = 0.e0
	2274	zhv_b(:,:) = 0.e0
	2275
	2276	IF ( ll_shorizd ) THEN
	2277	! Correct bottom value
	2278	! considering "shelf horizon depth"
	2279	DO jj = 1, jpjm1
	2280	DO ji = 1, jpim1 ! vector opt.
	2281	zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))
	2282	zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
	2283	DO jk=1, jpkm1
	2284	zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji ,jj ,jk), 0._wp)
	2285	zfp_hi = MIN(zfp_hi, e3t_b(ji ,jj ,jk))
	2286	zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) )
	2287	!
	2288	zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj ,jk), 0._wp)
	2289	zfm_hi = MIN(zfm_hi, e3t_b(ji+1,jj ,jk))
	2290	zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) )
	2291	!
	2292	zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji ,jj ,jk), 0._wp)
	2293	zfp_hj = MIN(zfp_hj, e3t_b(ji ,jj ,jk))
	2294	zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) )
	2295	!
	2296	zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji ,jj+1,jk), 0._wp)
	2297	zfm_hj = MIN(zfm_hj, e3t_b(ji ,jj+1,jk))
	2298	zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) )
	2299	!
	2300	zhu_b(ji,jj) = zhu_b(ji,jj) + ( zdi * zfp_hi &
	2301	& + (1._wp-zdi) * zfm_hi &
	2302	& ) * umask(ji,jj,jk)
	2303	zhv_b(ji,jj) = zhv_b(ji,jj) + ( zdj * zfp_hj &
	2304	& + (1._wp-zdj) * zfm_hj &
	2305	& ) * vmask(ji,jj,jk)
	2306	END DO
	2307	END DO
	2308	END DO
	2309	ELSE
	2310	DO jj = 1, jpjm1
	2311	DO ji = 1, jpim1 ! vector opt.
	2312	zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))
	2313	zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
	2314	DO jk = 1, jpkm1
	2315	zfp_hi = e3t_b(ji ,jj ,jk)
	2316	zfm_hi = e3t_b(ji+1,jj ,jk)
	2317	zfp_hj = e3t_b(ji ,jj ,jk)
	2318	zfm_hj = e3t_b(ji ,jj+1,jk)
	2319	!
	2320	zhu_b(ji,jj) = zhu_b(ji,jj) + ( zdi * zfp_hi &
	2321	& + (1._wp-zdi) * zfm_hi &
	2322	& ) * umask(ji,jj,jk)
	2323	!
	2324	zhv_b(ji,jj) = zhv_b(ji,jj) + ( zdj * zfp_hj &
	2325	& + (1._wp-zdj) * zfm_hj &
	2326	& ) * vmask(ji,jj,jk)
	2327	END DO
	2328	END DO
	2329	END DO
	2330	ENDIF
	2331
[10128]	2332	CALL lbc_lnk_multi( zhu_b(:,:), 'U', 1., zhv_b(:,:), 'V', 1. ) !* local domain boundaries
	2333
[10116]	2334	! Corrective barotropic velocity (times hor. scale factor)
	2335	zbu(:,:) = zbu(:,:)/ (zhu_b(:,:)*umask(:,:,1)+1._wp-umask(:,:,1))
	2336	zbv(:,:) = zbv(:,:)/ (zhv_b(:,:)*vmask(:,:,1)+1._wp-vmask(:,:,1))
	2337
	2338	! Set corrective fluxes in upstream direction:
	2339	!
	2340	zwx(:,:,:) = 0.e0
	2341	zwy(:,:,:) = 0.e0
[10128]	2342
[10116]	2343	IF ( ll_shorizd ) THEN
	2344	DO jj = 1, jpjm1
	2345	DO ji = 1, fs_jpim1 ! vector opt.
	2346	! upstream scheme
	2347	zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) )
	2348	zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) )
	2349	zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) )
	2350	zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) )
	2351	DO jk = 1, jpkm1
	2352	zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji ,jj ,jk), 0._wp)
	2353	zfp_hi = MIN(e3t_b(ji ,jj ,jk), zfp_hi)
	2354	zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) )
	2355	!
	2356	zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj ,jk), 0._wp)
	2357	zfm_hi = MIN(e3t_b(ji+1,jj ,jk), zfm_hi)
	2358	zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) )
	2359	!
	2360	zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji ,jj ,jk), 0._wp)
	2361	zfp_hj = MIN(e3t_b(ji ,jj ,jk), zfp_hj)
	2362	zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) )
	2363	!
	2364	zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji ,jj+1,jk), 0._wp)
	2365	zfm_hj = MIN(e3t_b(ji ,jj+1,jk), zfm_hj)
	2366	zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) )
	2367	!
	2368	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
	2369	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
	2370	END DO
	2371	END DO
	2372	END DO
	2373	ELSE
	2374	DO jj = 1, jpjm1
	2375	DO ji = 1, fs_jpim1 ! vector opt.
	2376	! upstream scheme
	2377	zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) )
	2378	zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) )
	2379	zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) )
	2380	zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) )
	2381	DO jk = 1, jpkm1
	2382	zfp_hi = e3t_b(ji ,jj ,jk)
	2383	zfm_hi = e3t_b(ji+1,jj ,jk)
	2384	zfp_hj = e3t_b(ji ,jj ,jk)
	2385	zfm_hj = e3t_b(ji ,jj+1,jk)
	2386	!
	2387	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
	2388	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
	2389	END DO
	2390	END DO
	2391	END DO
	2392	ENDIF
	2393
[10128]	2394	CALL lbc_lnk_multi( zwx, 'U', -1., zwy, 'V', -1. ) !* local domain boundaries
	2395
[10116]	2396	uin(:,:,:) = uin(:,:,:) + zwx(:,:,:)
	2397	vin(:,:,:) = vin(:,:,:) + zwy(:,:,:)
	2398	!
	2399	! Update trend with corrective fluxes:
	2400	DO jk = 1, jpkm1
	2401	DO jj = 2, jpjm1
	2402	DO ji = fs_2, fs_jpim1 ! vector opt.
	2403	!
	2404	zbtr = r1_e1e2t(ji,jj)
	2405	! total advective trends
	2406	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	2407	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) )
	2408	! add them to the general tracer trends
	2409	pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
	2410	END DO
	2411	END DO
	2412	END DO
	2413	!
	2414	IF( ln_timing ) CALL timing_stop('dom_vvl_ups_cor')
	2415	!
	2416	END SUBROUTINE dom_vvl_ups_cor
	2417
	2418	SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, p2dt )
	2419	!!---------------------------------------------------------------------
	2420	!! * ROUTINE nonosc_2d *
	2421	!!
	2422	!! ** Purpose : compute monotonic thickness fluxes from the upstream
	2423	!! scheme and the before field by a nonoscillatory algorithm
	2424	!!
	2425	!! ** Method : ... ???
	2426	!! warning : pbef and paft must be masked, but the boundaries
	2427	!! conditions on the fluxes are not necessary zalezak (1979)
	2428	!! drange (1995) multi-dimensional forward-in-time and upstream-
	2429	!! in-space based differencing for fluid
	2430	!!----------------------------------------------------------------------
	2431	!
	2432	!!----------------------------------------------------------------------
	2433	REAL(wp) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step
	2434	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field
	2435	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb ! monotonic fluxes in the 3 directions
	2436	!
	2437	INTEGER :: ji, jj, jk ! dummy loop indices
	2438	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars
	2439	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
	2440	REAL(wp) :: zupip1, zupim1, zupjp1, zupjm1, zupb, zupa
	2441	REAL(wp) :: zdoip1, zdoim1, zdojp1, zdojm1, zdob, zdoa
	2442	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo
	2443	!!----------------------------------------------------------------------
	2444	!
	2445	IF( ln_timing ) CALL timing_start('nonosc2')
	2446	!
	2447	zbig = 1.e+40_wp
	2448	zrtrn = 1.e-15_wp
	2449	zbetup(:,:,jpk) = 0._wp ; zbetdo(:,:,jpk) = 0._wp
	2450
	2451
	2452	! Search local extrema
	2453	! --------------------
	2454	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
	2455	zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ), &
	2456	& paft * tmask - zbig * ( 1.e0 - tmask ) )
	2457	zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ), &
	2458	& paft * tmask + zbig * ( 1.e0 - tmask ) )
	2459
	2460	DO jk = 1, jpkm1
	2461	z2dtt = p2dt
	2462	DO jj = 2, jpjm1
	2463	DO ji = fs_2, fs_jpim1 ! vector opt.
	2464
	2465	! search maximum in neighbourhood
	2466	zup = MAX( zbup(ji ,jj ,jk ), &
	2467	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
	2468	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ))
	2469
	2470	! search minimum in neighbourhood
	2471	zdo = MIN( zbdo(ji ,jj ,jk ), &
	2472	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
	2473	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ))
	2474
	2475	! positive part of the flux
	2476	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
	2477	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) )
	2478
	2479	! negative part of the flux
	2480	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
	2481	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) )
	2482
	2483	! up & down beta terms
	2484	zbt = e1t(ji,jj) * e2t(ji,jj) / z2dtt
	2485	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
	2486	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
	2487	END DO
	2488	END DO
	2489	END DO
[10128]	2490	CALL lbc_lnk_multi( zbetup, 'T', 1. , zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
[10116]	2491
	2492	! 3. monotonic flux in the i & j direction (paa & pbb)
	2493	! ----------------------------------------
	2494	DO jk = 1, jpkm1
	2495	DO jj = 2, jpjm1
	2496	DO ji = fs_2, fs_jpim1 ! vector opt.
	2497	zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
	2498	zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
	2499	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
	2500	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
	2501
	2502	zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
	2503	zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
	2504	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
	2505	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
	2506	END DO
	2507	END DO
	2508	END DO
[10128]	2509	CALL lbc_lnk_multi( paa, 'U', -1., pbb, 'V', -1. ) !* local domain boundaries
[10116]	2510	!
	2511	IF( ln_timing ) CALL timing_stop('nonosc2')
	2512	!
	2513	END SUBROUTINE nonosc_2d
	2514
[592]	2515	!!======================================================================
	2516	END MODULE domvvl

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

Download in other formats: