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

Line
1	MODULE domvvl
2	!!======================================================================
3	!! * MODULE domvvl *
4	!! Ocean :
5	!!======================================================================
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
8	!! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates
9	!! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability
10	!! 4.0 ! 2018-09 (J. Chanut) improve z_tilde robustness
11	!!----------------------------------------------------------------------
12
13	!!----------------------------------------------------------------------
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	!!----------------------------------------------------------------------
21	USE oce ! ocean dynamics and tracers
22	USE phycst ! physical constant
23	USE dom_oce ! ocean space and time domain
24	USE sbc_oce ! ocean surface boundary condition
25	USE wet_dry ! wetting and drying
26	USE usrdef_istate ! user defined initial state (wad only)
27	USE restart ! ocean restart
28	!
29	USE in_out_manager ! I/O manager
30	USE iom ! I/O manager library
31	USE lib_mpp ! distributed memory computing library
32	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
33	USE timing ! Timing
34	USE bdy_oce ! ocean open boundary conditions
35	USE sbcrnf ! river runoff
36	USE dynspg_ts, ONLY: un_adv, vn_adv
37
38	IMPLICIT NONE
39	PRIVATE
40
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
45
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
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
62	! ! conservation: not used yet
63	INTEGER :: nn_filt_order=1
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]
70
71	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport
72	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: un_lf, vn_lf, hdivn_lf ! low frequency fluxes and divergence
73	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors
74	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a ! baroclinic scale factors
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
80
81	!! * Substitutions
82	# include "vectopt_loop_substitute.h90"
83	!!----------------------------------------------------------------------
84	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
85	!! $Id$
86	!! Software governed by the CeCILL license (see ./LICENSE)
87	!!----------------------------------------------------------------------
88	CONTAINS
89
90	INTEGER FUNCTION dom_vvl_alloc()
91	!!----------------------------------------------------------------------
92	!! * FUNCTION dom_vvl_alloc *
93	!!----------------------------------------------------------------------
94	IF( ln_vvl_zstar ) dom_vvl_alloc = 0
95	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
96	ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , &
97	& un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , &
98	& tildemask(jpi,jpj) , hsm(jpi,jpj) , dsm(jpi,jpj) , i_int_bot(jpi,jpj), STAT = dom_vvl_alloc )
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')
101	un_td = 0._wp
102	vn_td = 0._wp
103	ENDIF
104	IF( ln_vvl_ztilde ) THEN
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 )
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
110	!
111	END FUNCTION dom_vvl_alloc
112
113
114	SUBROUTINE dom_vvl_init
115	!!----------------------------------------------------------------------
116	!! * ROUTINE dom_vvl_init *
117	!!
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	!!
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
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.
135	!!----------------------------------------------------------------------
136	INTEGER :: ji, jj, jk
137	INTEGER :: ii0, ii1, ij0, ij1
138	REAL(wp):: zcoef, zwgt, ztmp, zhmin, zhmax
139	!!----------------------------------------------------------------------
140	!
141	IF(lwp) WRITE(numout,*)
142	IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
143	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
144	!
145	CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer)
146	!
147	! ! Allocate module arrays
148	IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
149	!
150	! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdivn_lf
151	CALL dom_vvl_rst( nit000, 'READ' )
152	e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all
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' )
160	CALL dom_vvl_interpol( e3t_n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F
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' )
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(:,:,:)
173	!
174	! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep)
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
180	DO jk = 2, jpk ! vertical sum
181	DO jj = 1,jpj
182	DO ji = 1,jpi
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))
195	END DO
196	END DO
197	END DO
198	!
199	! !== thickness of the water column !! (ocean portion only)
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)
205	DO jk = 2, jpkm1
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)
211	END DO
212	!
213	! !== inverse of water column thickness ==! (u- and v- points)
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
219	! !== z_tilde coordinate case ==! (Restoring frequencies)
220	tildemask(:,:) = 1._wp
221
222	IF( ln_vvl_ztilde ) THEN
223	!!gm : idea: add here a READ in a file of custumized restoring frequency
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 )
227	!
228	IF( ln_vvl_ztilde_as_zstar ) THEN
229	! Ignore namelist settings and use these next two to emulate z-star using z-tilde
230	frq_rst_e3t(:,:) = 0.0_wp
231	frq_rst_hdv(:,:) = 1.0_wp / rdt
232	rn_lf_cutoff = 2.0_wp * rpi * rdt / 86400._wp
233	tildemask(:,:) = 0._wp
234	ENDIF
235
236	IF ( ln_vvl_zstar_at_eqtor ) THEN
237	DO jj = 1, jpj
238	DO ji = 1, jpi
239	!!gm case \|gphi\| >= 6 degrees is useless initialized just above by default
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 )
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
246	! values inside the equatorial band (ztilde as zstar)
247	frq_rst_e3t(ji,jj) = 0.0_wp
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)
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 ) )
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 ) )
261	ENDIF
262	END DO
263	END DO
264	ENDIF
265	!
266	IF ( ln_vvl_zstar_on_shelf ) THEN
267	zhmin = 50._wp
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	!
288	ENDIF
289
290	IF( ln_vvl_layer ) THEN
291	IF ( ln_vvl_zstar_on_shelf ) THEN
292	zhmin = 50._wp
293	zhmax = 100._wp
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
324	!
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	! ! ------------ !
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')
341	ENDIF
342	!
343	ENDIF
344	!
345	END SUBROUTINE dom_vvl_init
346
347
348	SUBROUTINE dom_vvl_sf_nxt( kt, kcall )
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	!!
363	!! ** Action : - hdivn_lf : restoring towards full baroclinic divergence in z_tilde case
364	!! - tilde_e3t_a: after increment of vertical scale factor
365	!! in z_tilde case
366	!! - e3(t/u/v)_a
367	!!
368	!! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling.
369	!!----------------------------------------------------------------------
370	INTEGER, INTENT( in ) :: kt ! time step
371	INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence
372	!
373	LOGICAL :: ll_do_bclinic ! local logical
374	INTEGER :: ji, jj, jk, jo ! dummy loop indices
375	INTEGER :: ib, ib_bdy, ip, jp ! " " "
376	INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers
377	INTEGER :: ncall
378	REAL(wp) :: z2dt , z_tmin, z_tmax! local scalars
379	REAL(wp) :: zalpha, zwgt ! temporary scalars
380	REAL(wp) :: zdu, zdv, zramp, zmet
381	REAL(wp) :: ztra, zbtr, ztout, ztin, zfac, zmsku, zmskv
382	REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv
383	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, ztu, ztv
384	!!----------------------------------------------------------------------
385	!
386	IF( ln_linssh ) RETURN ! No calculation in linear free surface
387	!
388	IF( ln_timing ) CALL timing_start('dom_vvl_sf_nxt')
389	!
390	IF( kt == nit000 ) THEN
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
396	zmet = 1._wp
397
398	ll_do_bclinic = .TRUE.
399	ncall = 1
400
401	IF( PRESENT(kcall) ) THEN
402	! comment line below if tilda coordinate has to be computed at each call
403	IF (kcall == 2 .AND. (ln_vvl_ztilde.OR.ln_vvl_layer) ) ll_do_bclinic = .FALSE.
404	ncall = kcall
405	ENDIF
406
407	IF( neuler == 0 .AND. kt == nit000 ) THEN
408	z2dt = rdt
409	ELSE
410	z2dt = 2.0_wp * rdt
411	ENDIF
412
413	! ******************************* !
414	! After scale factors at t-points !
415	! ******************************* !
416	! ! --------------------------------------------- !
417	! ! z_star coordinate and barotropic z-tilde part !
418	! ! --------------------------------------------- !
419	!
420	z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
421	DO jk = 1, jpkm1
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)
424	END DO
425	!
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
432	IF( (ln_vvl_ztilde.OR.ln_vvl_layer) .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !
433	! ! ------baroclinic part------ !
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
437
438	zhdiv(:,:) = 0.
439	DO jk = 1, jpkm1
440	zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
441	END DO
442	zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) )
443
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 )
458	DO jk = 1, jpkm1
459	ztu(:,:,jk) = un(:,:,jk) * e3u_n(:,:,jk) * e2u(:,:)
460	ztv(:,:,jk) = vn(:,:,jk) * e3v_n(:,:,jk) * e1v(:,:)
461	ze3t(:,:,jk) = -tilde_e3t_a(:,:,jk) - (e3t_n(:,:,jk)-zmettilde_e3t_n(:,:,jk)) zhdiv(:,:)
462	END DO
463	!
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
473	ENDIF
474	!
475	DO jk = 1, jpkm1
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)
479	END DO
480	!
481	ELSEIF ( ln_vvl_layer ) THEN
482	!
483	DO jk = 1, jpkm1
484	ztu(:,:,jk) = un(:,:,jk)
485	ztv(:,:,jk) = vn(:,:,jk)
486	END DO
487	!
488	ENDIF
489	!
490	! Block fluxes through small layers:
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
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
512	DO jk = 1, jpkm1
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
521	END DO
522	ENDIF
523	!
524	! Thickness anomaly diffusion:
525	! ----------------------------
526	ztu(:,:,:) = 0.0_wp
527	ztv(:,:,:) = 0.0_wp
528
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
535	END DO
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
576	ENDIF
577
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
589	END DO
590	END DO
591	ENDIF
592
593	! divergence of diffusive fluxes
594	DO jk = 1, jpkm1
595	DO jj = 2, jpjm1
596	DO ji = fs_2, fs_jpim1 ! vector opt.
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 ) &
603	& ) * r1_e1e2t(ji,jj)
604	END DO
605	END DO
606	END DO
607
608	CALL lbc_lnk_multi( un_td, 'U', -1., vn_td, 'V', -1. ) !* local domain boundaries
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
620	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + (e3t_n(:,:,jk)-zmettilde_e3t_n(:,:,jk)) zhdiv(:,:) * tmask(:,:,jk)
621	END DO
622
623	! Leapfrog time stepping
624	! ~~~~~~~~~~~~~~~~~~~~~~
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)
627
628	DO jk=1,jpkm1
629	tilde_e3t_a(:,:,jk) = tilde_e3t_b(:,:,jk) + z2dt * tmask(:,:,jk) * tilde_e3t_a(:,:,jk) &
630	& * tildemask(:,:) * zramp
631	END DO
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
655	ENDIF
656
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 )
663	!
664	! divergence of diffusive fluxes
665	DO jk = 1, jpkm1
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
673	END DO
674	CALL lbc_lnk( ze3t(:,:,:), 'T', 1. )
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
679	ENDIF
680
681	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate !
682	! ! ---baroclinic part--------- !
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	!
711	DO jk = 1, jpkm1
712	e3t_a(:,:,jk) = e3t_a(:,:,jk) + tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
713	END DO
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
721	ENDIF
722
723	IF( ln_vvl_dbg.AND.(ln_vvl_ztilde .OR. ln_vvl_layer) ) THEN ! - ML - test: control prints for debuging
724	!
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
732	IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
733	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ), mask = tmask(:,:,1) == 1.e0 )
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	!
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
750	ji = ijk_min(1) ; jj = ijk_min(2) ; jk = ijk_min(3)
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	!
811	zht(:,:) = 0.0_wp
812	DO jk = 1, jpkm1
813	zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
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
817	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn -SUM(e3t_n))) =', z_tmax
818	!
819	zht(:,:) = 0.0_wp
820	DO jk = 1, jpkm1
821	zht(:,:) = zht(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
822	END DO
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(:,:) ) )
832	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
833	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hu_0+sshu_n-SUM(e3u_n))) =', z_tmax
834	!
835	zht(:,:) = 0.0_wp
836	DO jk = 1, jpkm1
837	zht(:,:) = zht(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
838	END DO
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(:,:) ) )
848	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
849	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hv_0+sshv_n-SUM(e3v_n))) =', z_tmax
850	!
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(:,:) ) )
868	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
869	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hf_0+sshf_n-SUM(e3f_n))) =', z_tmax
870	!
871	END IF
872
873	! *********************************** !
874	! After scale factors at u- v- points !
875	! *********************************** !
876
877	CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
878	CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
879
880	! *********************************** !
881	! After depths at u- v points !
882	! *********************************** !
883
884	hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
885	hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
886	DO jk = 2, jpkm1
887	hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
888	hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
889	END DO
890	! ! Inverse of the local depth
891	!!gm BUG ? don't understand the use of umask_i here .....
892	r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
893	r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
894	!
895	IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt')
896	!
897	END SUBROUTINE dom_vvl_sf_nxt
898
899
900	SUBROUTINE dom_vvl_sf_swp( kt )
901	!!----------------------------------------------------------------------
902	!! * ROUTINE dom_vvl_sf_swp *
903	!!
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
907	!!
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	!!
912	!! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
913	!! - Recompute:
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
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.
923	!!----------------------------------------------------------------------
924	INTEGER, INTENT( in ) :: kt ! time step
925	!
926	INTEGER :: ji, jj, jk ! dummy loop indices
927	REAL(wp) :: zcoef ! local scalar
928	!!----------------------------------------------------------------------
929	!
930	IF( ln_linssh ) RETURN ! No calculation in linear free surface
931	!
932	IF( ln_timing ) CALL timing_start('dom_vvl_sf_swp')
933	!
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'
938	ENDIF
939	!
940	! Time filter and swap of scale factors
941	! =====================================
942	! - ML - e3(t/u/v)_b are allready computed in dynnxt.
943	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
944	IF( neuler == 0 .AND. kt == nit000 ) THEN
945	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
946	ELSE
947	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &
948	& + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
949	ENDIF
950	tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
951	ENDIF
952	gdept_b(:,:,:) = gdept_n(:,:,:)
953	gdepw_b(:,:,:) = gdepw_n(:,:,:)
954
955	e3t_n(:,:,:) = e3t_a(:,:,:)
956	e3u_n(:,:,:) = e3u_a(:,:,:)
957	e3v_n(:,:,:) = e3v_a(:,:,:)
958
959	! Compute all missing vertical scale factor and depths
960	! ====================================================
961	! Horizontal scale factor interpolations
962	! --------------------------------------
963	! - ML - e3u_b and e3v_b are allready computed in dynnxt
964	! - JC - hu_b, hv_b, hur_b, hvr_b also
965
966	CALL dom_vvl_interpol( e3t_n(:,:,:), e3f_n(:,:,:), 'F' )
967
968	! Vertical scale factor interpolations
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' )
975
976	! t- and w- points depth (set the isf depth as it is in the initial step)
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(:,:)
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))
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)
990	END DO
991	END DO
992	END DO
993
994	! Local depth and Inverse of the local depth of the water
995	! -------------------------------------------------------
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)
1000	DO jk = 2, jpkm1
1001	ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
1002	END DO
1003
1004	! Output some diagnostics:
1005	! ------------------------
1006	IF (ln_vvl_ztilde .OR. ln_vvl_layer) CALL dom_vvl_dia( kt )
1007
1008	! write restart file
1009	! ==================
1010	IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
1011	!
1012	IF( ln_timing ) CALL timing_stop('dom_vvl_sf_swp')
1013	!
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	!!----------------------------------------------------------------------
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	!
1032	INTEGER :: ji, jj, jk, jkbot ! dummy loop indices
1033	INTEGER :: nmet ! horizontal interpolation method
1034	REAL(wp) :: zlnwd ! =1./0. when ln_wd_il = T/F
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
1040	!!----------------------------------------------------------------------
1041	!
1042	! nmet = 0 ! Original method (Surely wrong)
1043	nmet = 2 ! Interface interpolation
1044	! nmet = 2 ! Internal interfaces interpolation only, spread barotropic increment
1045	! Note that we kept surface weighted interpolation for barotropic increment to be compliant
1046	! with what is done in surface pressure module.
1047	!
1048	IF(ln_wd_il) THEN
1049	zlnwd = 1.0_wp
1050	ELSE
1051	zlnwd = 0.0_wp
1052	END IF
1053	!
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)
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
1087	zw(:,:,:) = (zw(:,:,:) - gdepw_0(:,:,:))*tmask(:,:,:)
1088	!
1089	END SELECT
1090	END IF
1091	!
1092	pe3_out(:,:,:) = 0.0_wp
1093	!
1094	SELECT CASE ( pout ) !== type of interpolation ==!
1095	!
1096	CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean
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
1108	DO jj = 1, jpjm1
1109	DO ji = 1, fs_jpim1 ! vector opt.
1110	! Correction at last level:
1111	jkbot = mbku(ji,jj)
1112	zdo = 0._wp
1113	DO jk=jkbot,1,-1
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)
1116	!
1117	! If there is a step, taper bottom interface:
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 )
1128	zdo = zup
1129	END DO
1130	END DO
1131	END DO
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) ) &
1140	& / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) ) &
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	!
1150	CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
1151	pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
1152	!
1153	CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean
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
1165	DO jj = 1, jpjm1
1166	DO ji = 1, fs_jpim1 ! vector opt.
1167	! Correction at last level:
1168	jkbot = mbkv(ji,jj)
1169	zdo = 0._wp
1170	DO jk=jkbot,1,-1
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)
1173	!
1174	! If there is a step, taper bottom interface:
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 )
1185	zdo = zup
1186	END DO
1187	END DO
1188	END DO
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) ) &
1197	& / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) ) &
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	!
1207	CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
1208	pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
1209	!
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
1225	DO jj = 1, jpjm1
1226	DO ji = 1, fs_jpim1 ! vector opt.
1227	! bottom correction:
1228	jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1))
1229	zdo = 0._wp
1230	DO jk=jkbot,1,-1
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)
1235	!
1236	! If there is a step, taper bottom interface:
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)
1254	!
1255	pe3_out(ji,jj,jk) = ( zdo - zup ) &
1256	& ( umask(ji,jj,jk) umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd )
1257	zdo = zup
1258	END DO
1259	END DO
1260	END DO
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	!
1280	CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
1281	pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
1282	!
1283	CASE( 'W' ) !* from T- to W-point : vertical simple mean
1284	!
1285	pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
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
1288	DO jk = 2, jpk
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 ) )
1293	END DO
1294	!
1295	CASE( 'UW' ) !* from U- to UW-point : vertical simple mean
1296	!
1297	pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
1298	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
1299	!!gm BUG? use here wumask in case of ISF ? to be checked
1300	DO jk = 2, jpk
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 ) )
1305	END DO
1306	!
1307	CASE( 'VW' ) !* from V- to VW-point : vertical simple mean
1308	!
1309	pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
1310	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
1311	!!gm BUG? use here wvmask in case of ISF ? to be checked
1312	DO jk = 2, jpk
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 ) )
1317	END DO
1318	END SELECT
1319	!
1320	END SUBROUTINE dom_vvl_interpol
1321
1322
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
1337	!
1338	INTEGER :: ji, jj, jk
1339	INTEGER :: id1, id2, id3, id4, id5, id6, id7 ! local integers
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
1346	CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn, ldxios = lrxios )
1347	!
1348	id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
1349	id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
1350	id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
1351	id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
1352	id5 = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. )
1353	id6 = iom_varid( numror, 'un_lf', ldstop = .FALSE. )
1354	id7 = iom_varid( numror, 'vn_lf', ldstop = .FALSE. )
1355	! ! --------- !
1356	! ! all cases !
1357	! ! --------- !
1358	IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist
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 )
1361	! needed to restart if land processor not computed
1362	IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
1363	WHERE ( tmask(:,:,:) == 0.0_wp )
1364	e3t_n(:,:,:) = e3t_0(:,:,:)
1365	e3t_b(:,:,:) = e3t_0(:,:,:)
1366	END WHERE
1367	IF( neuler == 0 ) THEN
1368	e3t_b(:,:,:) = e3t_n(:,:,:)
1369	ENDIF
1370	ELSE IF( id1 > 0 ) THEN
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.'
1373	IF(lwp) write(numout,*) 'neuler is forced to 0'
1374	CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
1375	e3t_n(:,:,:) = e3t_b(:,:,:)
1376	neuler = 0
1377	ELSE IF( id2 > 0 ) THEN
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.'
1380	IF(lwp) write(numout,*) 'neuler is forced to 0'
1381	CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
1382	e3t_b(:,:,:) = e3t_n(:,:,:)
1383	neuler = 0
1384	ELSE
1385	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
1386	IF(lwp) write(numout,*) 'Compute scale factor from sshn'
1387	IF(lwp) write(numout,*) 'neuler is forced to 0'
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))
1392	END DO
1393	e3t_b(:,:,:) = e3t_n(:,:,:)
1394	neuler = 0
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
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 )
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	! ! ------------ !
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 )
1419	ELSE ! array is missing
1420	hdivn_lf(:,:,:,:) = 0.0_wp
1421	un_lf(:,:,:,:) = 0.0_wp
1422	vn_lf(:,:,:,:) = 0.0_wp
1423	ENDIF
1424	ENDIF
1425	ENDIF
1426	!
1427	ELSE !* Initialize at "rest"
1428	!
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
1457	DO jk = 1, jpk
1458	e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
1459	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
1460	& + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
1461	END DO
1462	e3t_b(:,:,:) = e3t_n(:,:,:)
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
1471	!
1472	ELSE
1473	!
1474	! Just to read set ssh in fact, called latter once vertical grid
1475	! is set up:
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(:,:,:)
1486	!
1487	END IF ! end of ll_wd edits
1488
1489	IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
1490	tilde_e3t_b(:,:,:) = 0._wp
1491	tilde_e3t_n(:,:,:) = 0._wp
1492	IF( ln_vvl_ztilde ) THEN
1493	hdivn_lf(:,:,:,:) = 0._wp
1494	un_lf(:,:,:,:) = 0._wp
1495	vn_lf(:,:,:,:) = 0._wp
1496	ENDIF
1497	END IF
1498	ENDIF
1499	!
1500	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
1501	! ! ===================
1502	IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
1503	IF( lwxios ) CALL iom_swap( cwxios_context )
1504	! ! --------- !
1505	! ! all cases !
1506	! ! --------- !
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 )
1509	! ! ----------------------- !
1510	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
1511	! ! ----------------------- !
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)
1514	END IF
1515	! ! -------------!
1516	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
1517	! ! ------------ !
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)
1521	ENDIF
1522	!
1523	IF( lwxios ) CALL iom_swap( cxios_context )
1524	ENDIF
1525	!
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	!!----------------------------------------------------------------------
1536	INTEGER :: ioptio, ios
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	!!----------------------------------------------------------------------
1549	!
1550	REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist :
1551	READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
1552	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
1553	REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run
1554	READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
1555	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
1556	IF(lwm) WRITE ( numond, nam_vvl )
1557	!
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,*) '~~~~~~~~~~~'
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
1569	WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor
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
1584	IF( ln_vvl_ztilde_as_zstar ) THEN
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'
1591	ELSE
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
1596	ENDIF
1597	WRITE(numout,*) ' Namelist nam_vvl : debug prints'
1598	WRITE(numout,*) ' ln_vvl_dbg = ', ln_vvl_dbg
1599	ENDIF
1600	!
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	!
1608	ioptio = 0 ! Parameter control
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	!
1614	IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1615	IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1616	!
1617	IF(lwp) THEN ! Print the choice
1618	WRITE(numout,*)
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'
1623	ENDIF
1624	!
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	!
1633	#if defined key_agrif
1634	IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) ) CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1635	#endif
1636	!
1637	END SUBROUTINE dom_vvl_ctl
1638
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.
1681	dzmin_int = 1.0_wp ! Absolute minimum depth in the interior (in meters)
1682	dzmin_surf = 1.0_wp ! Absolute minimum depth at the surface (in meters)
1683	zfrch_stp = 5._wp ! Maximum fractionnal thickness change in one time step (<= 1.)
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
1691	zhdiff = 0.01_wp ! ad.
1692	zhlim = 0.03_wp ! ad. max lap(z)*e1
1693	ll_blpdiff_cond = .TRUE. ! Conditionnal Bilaplacian diffusion of interfaces
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)
1739	dsm(ji,jj) = MAX(dsm(ji,jj), 0.05_wp*ht_0(ji,jj))
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)
1747	dsm(ji,jj) = MAX(dsm(ji,jj), 0.05_wp*ht_0(ji,jj))
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)
1895	! zh_min = MAX(zh_min, zh_min-e3t_a(ji,jj,jk)+e3t_0(ji,jj,jk))
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)
1913	! zh_new = cush(zh_new, zh_min)
1914	zh_new = MAX(zh_new, zh_min)
1915	!
1916	! Final flux:
1917	zdiff = (zh_new - zh_old) * tmask(ji,jj,jk)
1918	!
1919	! Limit thickness change in 1 time step:
1920	! ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
1921	! zdiff = SIGN(ztmp, zh_new - zh_old)
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)
1942	! zh_min = MAX(zh_min, zh_min-e3t_a(ji,jj,jk)+e3t_0(ji,jj,jk))
1943	!
1944	! zwdw(ji,jj,jk+1) = MAX(zwdw(ji,jj,jk+1), REAL(jk)*dzmin_surf)
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)
1951	! zh_new = cush(zh_new, zh_min)
1952	zh_new = MAX(zh_new, zh_min)
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
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
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
2196	CALL lbc_lnk_multi( zwx, 'U', -1., zwy, 'V', -1. ) !* local domain boundaries
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
2218	zwi(ji,jj,jk) = (e3t_b(ji,jj,jk) + z2dtt * pta(ji,jj,jk)* bdytmask(ji,jj) ) * tmask(ji,jj,jk)
2219	END DO
2220	END DO
2221	END DO
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
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
2276
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
2337	CALL lbc_lnk_multi( zhu_b(:,:), 'U', 1., zhv_b(:,:), 'V', 1. ) !* local domain boundaries
2338
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
2347
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
2399	CALL lbc_lnk_multi( zwx, 'U', -1., zwy, 'V', -1. ) !* local domain boundaries
2400
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
2495	CALL lbc_lnk_multi( zbetup, 'T', 1. , zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
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
2514	CALL lbc_lnk_multi( paa, 'U', -1., pbb, 'V', -1. ) !* local domain boundaries
2515	!
2516	IF( ln_timing ) CALL timing_stop('nonosc2')
2517	!
2518	END SUBROUTINE nonosc_2d
2519
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
2696	!!======================================================================
2697	END MODULE domvvl

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

Download in other formats: