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

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

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

Download in other formats: