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

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

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

Download in other formats: