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 branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 9545

Last change on this file since 9545 was 9545, checked in by jchanut, 7 years ago
Add diagnostics in diawri - revert to old interface interpolation - change default regriding parameters - add linear ramp at startup
Property svn:keywords set to `Id`
File size: 158.0 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:
9	!! vvl option includes z_star and z_tilde coordinates
10	!! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability
11	!! ! 2018-01 (J. Chanut) improve ztilde robustness
12	!!----------------------------------------------------------------------
13	!! 'key_vvl' variable volume
14	!!----------------------------------------------------------------------
15	!!----------------------------------------------------------------------
16	!! dom_vvl_init : define initial vertical scale factors, depths and column thickness
17	!! dom_vvl_sf_nxt : Compute next vertical scale factors
18	!! dom_vvl_sf_swp : Swap vertical scale factors and update the vertical grid
19	!! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
20	!! dom_vvl_rst : read/write restart file
21	!! dom_vvl_ctl : Check the vvl options
22	!! dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors
23	!! : to account for manual changes to e[1,2][u,v] in some Straits
24	!!----------------------------------------------------------------------
25	!! * Modules used
26	USE oce ! ocean dynamics and tracers
27	USE dom_oce ! ocean space and time domain
28	USE sbc_oce ! ocean surface boundary condition
29	USE in_out_manager ! I/O manager
30	USE iom ! I/O manager library
31	USE restart ! ocean restart
32	USE lib_mpp ! distributed memory computing library
33	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
34	USE wrk_nemo ! Memory allocation
35	USE timing ! Timing
36	USE bdy_oce ! ocean open boundary conditions
37	USE sbcrnf ! river runoff
38
39	IMPLICIT NONE
40	PRIVATE
41
42	!! * Routine accessibility
43	PUBLIC dom_vvl_init ! called by domain.F90
44	PUBLIC dom_vvl_sf_nxt ! called by step.F90
45	PUBLIC dom_vvl_sf_swp ! called by step.F90
46	PUBLIC dom_vvl_interpol ! called by dynnxt.F90
47	PRIVATE dom_vvl_orca_fix ! called by dom_vvl_interpol
48
49	!!* Namelist nam_vvl
50	LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate
51	LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate
52	LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate
53	LOGICAL :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate
54	LOGICAL :: ln_vvl_zstar_at_eqtor = .FALSE. ! revert to zstar at equator
55	LOGICAL :: ln_vvl_zstar_on_shelf =.FALSE. ! revert to zstar on shelves
56	LOGICAL :: ln_vvl_kepe =.FALSE. ! kinetic/potential energy transfer
57	LOGICAL :: ln_vvl_adv_fct =.FALSE. ! Centred thickness advection
58	LOGICAL :: ln_vvl_adv_cn2 =.TRUE. ! FCT thickness advection
59	LOGICAL :: ln_vvl_dbg = .FALSE. ! debug control prints
60	LOGICAL :: ln_vvl_ramp = .FALSE. ! Ramp on interfaces displacement
61	LOGICAL :: ln_vvl_lap ! Laplacian thickness diffusion
62	LOGICAL :: ln_vvl_blp ! Bilaplacian thickness diffusion
63	LOGICAL :: ln_vvl_regrid ! ensure layer separation
64	LOGICAL :: ll_shorizd=.FALSE. ! Use "shelf horizon depths"
65	! ! conservation: not used yet
66	REAL(wp) :: rn_ahe3_lap ! thickness diffusion coefficient (Laplacian)
67	REAL(wp) :: rn_ahe3_blp ! thickness diffusion coefficient (Bilaplacian)
68	REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days]
69	REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days]
70	REAL(wp) :: rn_day_ramp ! Duration of linear ramp [days]
71	REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation
72	REAL(wp) :: hsmall=0.01 ! small thickness
73
74	!! * Module variables
75	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport
76	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_lf, vn_lf, hdivn_lf ! 1st order filtered arrays
77	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors
78	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors
79	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! restoring period for scale factors
80	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! restoring period for low freq. divergence
81	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: tildemask ! mask tilde tendency
82	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: hsm, dsm !
83	INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: i_int_bot
84
85	!! * Substitutions
86	# include "domzgr_substitute.h90"
87	# include "vectopt_loop_substitute.h90"
88	!!----------------------------------------------------------------------
89	!! NEMO/OPA 3.3 , NEMO-Consortium (2010)
90	!! $Id$
91	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
92	!!----------------------------------------------------------------------
93
94	CONTAINS
95
96	INTEGER FUNCTION dom_vvl_alloc()
97	!!----------------------------------------------------------------------
98	!! * FUNCTION dom_vvl_alloc *
99	!!----------------------------------------------------------------------
100	IF( ln_vvl_zstar ) dom_vvl_alloc = 0
101	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
102	ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , &
103	& dtilde_e3t_a(jpi,jpj,jpk) , un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , &
104	& i_int_bot(jpi,jpj), tildemask(jpi,jpj), hsm(jpi,jpj), dsm(jpi,jpj), &
105	& 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	un_td = 0.0_wp
109	vn_td = 0.0_wp
110	ENDIF
111	IF( ln_vvl_ztilde ) THEN
112	ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj), &
113	& hdivn_lf(jpi,jpj,jpk), &
114	& un_lf(jpi,jpj,jpk), &
115	& vn_lf(jpi,jpj,jpk), STAT= dom_vvl_alloc )
116	IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc )
117	IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
118	ENDIF
119
120	END FUNCTION dom_vvl_alloc
121
122
123	SUBROUTINE dom_vvl_init
124	!!----------------------------------------------------------------------
125	!! * ROUTINE dom_vvl_init *
126	!!
127	!! ** Purpose : Initialization of all scale factors, depths
128	!! and water column heights
129	!!
130	!! ** Method : - use restart file and/or initialize
131	!! - interpolate scale factors
132	!!
133	!! ** Action : - fse3t_(n/b) and tilde_e3t_(n/b)
134	!! - Regrid: fse3(u/v)_n
135	!! fse3(u/v)_b
136	!! fse3w_n
137	!! fse3(u/v)w_b
138	!! fse3(u/v)w_n
139	!! fsdept_n, fsdepw_n and fsde3w_n
140	!! - h(t/u/v)_0
141	!! - frq_rst_e3t and frq_rst_hdv
142	!!
143	!! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling.
144	!!----------------------------------------------------------------------
145	USE phycst, ONLY : rpi, rsmall, rad
146	!! * Local declarations
147	INTEGER :: ji, jj, jk
148	INTEGER :: ii0, ii1, ij0, ij1
149	REAL(wp):: zcoef, zwgt, ztmp, zhmin, zhmax
150	!!----------------------------------------------------------------------
151	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init')
152
153	IF(lwp) WRITE(numout,*)
154	IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
155	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
156
157	! choose vertical coordinate (z_star, z_tilde or layer)
158	! ==========================
159	CALL dom_vvl_ctl
160
161	! Allocate module arrays
162	! ======================
163	IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
164
165	! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk))
166	! ============================================================================
167	CALL dom_vvl_rst( nit000, 'READ' )
168	fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
169
170	! Reconstruction of all vertical scale factors at now and before time steps
171	! =============================================================================
172	! Horizontal scale factor interpolations
173	! --------------------------------------
174	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
175	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
176	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
177	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
178	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n(:,:,:), 'F' )
179	! Vertical scale factor interpolations
180	! ------------------------------------
181	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )
182	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
183	CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
184	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' )
185	CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
186	CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
187	! t- and w- points depth
188	! ----------------------
189	! set the isf depth as it is in the initial step
190	fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
191	fsdepw_n(:,:,1) = 0.0_wp
192	fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
193	fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1)
194	fsdepw_b(:,:,1) = 0.0_wp
195
196	DO jk = 2, jpk
197	DO jj = 1,jpj
198	DO ji = 1,jpi
199	! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
200	! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
201	! 0.5 where jk = mikt
202	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
203	fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
204	fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) &
205	& + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))
206	fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
207	fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1)
208	fsdept_b(ji,jj,jk) = zcoef * ( fsdepw_b(ji,jj,jk ) + 0.5 * fse3w_b(ji,jj,jk)) &
209	& + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk))
210	END DO
211	END DO
212	END DO
213
214	! Before depth and Inverse of the local depth of the water column at u- and v- points
215	! -----------------------------------------------------------------------------------
216	hu_b(:,:) = 0.
217	hv_b(:,:) = 0.
218	DO jk = 1, jpkm1
219	hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
220	hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
221	END DO
222	hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
223	hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
224
225	! Restoring frequencies for z_tilde coordinate
226	! ============================================
227	tildemask(:,:) = 1._wp
228
229	IF( ln_vvl_ztilde ) THEN
230	! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings
231	frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp )
232	frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
233	!
234	IF( ln_vvl_ztilde_as_zstar ) THEN
235	! Ignore namelist settings and use these next two to emulate z-star using z-tilde
236	frq_rst_e3t(:,:) = 0.0_wp
237	frq_rst_hdv(:,:) = 1.0_wp / rdt
238	tildemask(:,:) = 0._wp
239	ENDIF
240
241	IF ( ln_vvl_zstar_at_eqtor ) THEN
242	DO jj = 1, jpj
243	DO ji = 1, jpi
244	IF( ABS(gphit(ji,jj)) >= 6.) THEN
245	! values outside the equatorial band and transition zone (ztilde)
246	frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp )
247	! frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
248
249	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
250	! values inside the equatorial band (ztilde as zstar)
251	frq_rst_e3t(ji,jj) = 0.0_wp
252	! frq_rst_hdv(ji,jj) = 1.0_wp / rdt
253	tildemask(ji,jj) = 0._wp
254	ELSE
255	! values in the transition band (linearly vary from ztilde to ztilde as zstar values)
256	frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp &
257	& * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
258	& * 180._wp / 3.5_wp ) )
259	! frq_rst_hdv(ji,jj) = (1.0_wp / rdt) &
260	! & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp &
261	! & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
262	! & * 180._wp / 3.5_wp ) )
263	tildemask(ji,jj) = 0.5_wp * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
264	& * 180._wp / 3.5_wp ) )
265	ENDIF
266	END DO
267	END DO
268	IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN
269	ii0 = 103 ; ii1 = 111 ! Suppress ztilde in the Foxe Basin for ORCA2
270	ij0 = 128 ; ij1 = 135 ;
271	frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp
272	tildemask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp
273	! frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rdt
274	ENDIF
275	ENDIF
276	!
277	IF ( ln_vvl_zstar_on_shelf ) THEN
278	zhmin = 50._wp
279	zhmax = 100._wp
280	DO jj = 1, jpj
281	DO ji = 1, jpi
282	zwgt = 1._wp
283	IF(( ht_0(ji,jj)>zhmin).AND.(ht_0(ji,jj) <=zhmax)) THEN
284	zwgt = (ht_0(ji,jj)-zhmin)/(zhmax-zhmin)
285	ELSEIF ( ht_0(ji,jj)<=zhmin) THEN
286	zwgt = 0._wp
287	ENDIF
288	frq_rst_e3t(ji,jj) = MIN(frq_rst_e3t(ji,jj), frq_rst_e3t(ji,jj)*zwgt)
289	tildemask(ji,jj) = MIN(tildemask(ji,jj), zwgt)
290	END DO
291	END DO
292	ENDIF
293	!
294	ztmp = MAXVAL( frq_rst_hdv(:,:) )
295	IF( lk_mpp ) CALL mpp_max( ztmp ) ! max over the global domain
296	!
297	IF ( (ztmp*rdt) > 1._wp) CALL ctl_stop( 'dom_vvl_init: rn_lf_cuttoff is too small' )
298	!
299	ENDIF
300
301	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_init')
302
303	END SUBROUTINE dom_vvl_init
304
305
306	SUBROUTINE dom_vvl_sf_nxt( kt, kcall )
307	!!----------------------------------------------------------------------
308	!! * ROUTINE dom_vvl_sf_nxt *
309	!!
310	!! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt,
311	!! tranxt and dynspg routines
312	!!
313	!! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness.
314	!! - z_tilde_case: after scale factor increment =
315	!! high frequency part of horizontal divergence
316	!! + retsoring towards the background grid
317	!! + thickness difusion
318	!! Then repartition of ssh INCREMENT proportionnaly
319	!! to the "baroclinic" level thickness.
320	!!
321	!! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case
322	!! - tilde_e3t_a: after increment of vertical scale factor
323	!! in z_tilde case
324	!! - fse3(t/u/v)_a
325	!!
326	!! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling.
327	!!----------------------------------------------------------------------
328	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t, ztu, ztv
329	REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv
330	!! * Arguments
331	INTEGER, INTENT( in ) :: kt ! time step
332	INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence
333	!! * Local declarations
334	INTEGER :: ji, jj, jk ! dummy loop indices
335	INTEGER :: ib, ib_bdy, ip, jp ! " " "
336	INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers
337	INTEGER :: ncall
338	REAL(wp) :: z2dt ! temporary scalars
339	REAL(wp) :: z_tmin, z_tmax ! temporary scalars
340	REAL(wp) :: zalpha, zwgt ! temporary scalars
341	REAL(wp) :: zdu, zdv, zramp
342	LOGICAL :: ll_do_bclinic ! temporary logical
343	REAL(wp) :: ztra, zbtr, ztout, ztin, zfac, zmsku, zmskv
344	!!----------------------------------------------------------------------
345	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt')
346	CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
347	CALL wrk_alloc( jpi, jpj, jpk, ze3t, ztu, ztv )
348
349	IF(kt == nit000) THEN
350	IF(lwp) WRITE(numout,*)
351	IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
352	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
353	ENDIF
354
355	ll_do_bclinic = .TRUE.
356	ncall=1
357	IF( PRESENT(kcall) ) THEN
358	! comment line below if tilda coordinate has to be computed at each call
359	IF (kcall == 2 .AND. (ln_vvl_ztilde.OR.ln_vvl_layer) ) ll_do_bclinic = .FALSE.
360	IF (kcall == 2) ncall=2
361	ENDIF
362
363	IF( neuler == 0 .AND. kt == nit000 ) THEN
364	z2dt = rdt
365	ELSE
366	z2dt = 2.0_wp * rdt
367	ENDIF
368
369	! ******************************* !
370	! After scale factors at t-points !
371	! ******************************* !
372
373	! ! --------------------------------------------- !
374	! z_star coordinate and barotropic z-tilde part !
375	! ! --------------------------------------------- !
376
377	z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - ssmask(:,:) )
378	DO jk = 1, jpkm1
379	! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
380	fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
381	END DO
382
383	IF((ln_vvl_ztilde .OR. ln_vvl_layer).AND.ll_do_bclinic ) THEN ! z_tilde or layer coordinate !
384
385	tilde_e3t_a(:,:,:) = 0.0_wp ! tilde_e3t_a used to store tendency terms
386	un_td(:,:,:) = 0.0_wp ! Transport corrections
387	vn_td(:,:,:) = 0.0_wp
388
389	zhdiv(:,:) = 0.
390	DO jk = 1, jpkm1
391	zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
392	END DO
393	zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) )
394
395	ze3t(:,:,:) = 0._wp
396	IF( ln_rnf ) THEN
397	CALL sbc_rnf_div( ze3t ) ! runoffs
398	DO jk=1,jpkm1
399	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ze3t(:,:,jk) * fse3t_n(:,:,jk)
400	END DO
401	ENDIF
402
403	! Thickness advection:
404	! --------------------
405	! Set advection velocities and source term
406	IF ( ln_vvl_ztilde ) THEN
407	IF ( ncall==1 ) THEN
408	zalpha = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
409	DO jk = 1, jpkm1
410	ztu(:,:,jk) = un(:,:,jk) * fse3u_n(:,:,jk) * e2u(:,:)
411	ztv(:,:,jk) = vn(:,:,jk) * fse3v_n(:,:,jk) * e1v(:,:)
412	ze3t(:,:,jk) = -tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * zhdiv(:,:)
413	END DO
414	!
415	un_lf(:,:,:) = un_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztu(:,:,:)
416	vn_lf(:,:,:) = vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztv(:,:,:)
417	hdivn_lf(:,:,:) = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:)
418	! un_lf(:,:,:) = un_lf(:,:,:) * (1._wp - zalpha) + zalpha * un_lf2(:,:,:)
419	! vn_lf(:,:,:) = vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * vn_lf2(:,:,:)
420	! hdivn_lf(:,:,:) = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * hdivn_lf2(:,:,:)
421	ENDIF
422	!
423	DO jk = 1, jpkm1
424	ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk)/fse3u_n(:,:,jk)r1_e2u(:,:))umask(:,:,jk)
425	ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk)/fse3v_n(:,:,jk)r1_e1v(:,:))vmask(:,:,jk)
426	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
427	END DO
428
429	!
430	ELSEIF ( ln_vvl_layer ) THEN
431	!
432	DO jk = 1, jpkm1
433	ztu(:,:,jk) = un(:,:,jk)
434	ztv(:,:,jk) = vn(:,:,jk)
435	END DO
436	!
437	ENDIF
438	!
439	! Block fluxes through small layers:
440	!! DO jk=1,jpkm1
441	!! DO ji = 1, jpi
442	!! DO jj= 1, jpj
443	!! zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, fse3u_n(ji,jj,jk) - hsmall) )
444	!! un_td(ji,jj,jk) = un_td(ji,jj,jk) - (1. - zmsku) * un(ji,jj,jk) * fse3u_n(ji,jj,jk) * e2u(ji,jj)
445	!! ztu(ji,jj,jk) = zmsku * ztu(ji,jj,jk)
446	!! IF ( ln_vvl_ztilde ) un_lf(ji,jj,jk) = zmsku * un_lf(ji,jj,jk)
447	!! !
448	!! zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, fse3v_n(ji,jj,jk) - hsmall) )
449	!! vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - (1. - zmskv) * vn(ji,jj,jk) * fse3v_n(ji,jj,jk) * e1v(ji,jj)
450	!! ztv(ji,jj,jk) = zmskv * ztv(ji,jj,jk)
451	!! IF ( ln_vvl_ztilde ) vn_lf(ji,jj,jk) = zmskv * vn_lf(ji,jj,jk)
452	!! END DO
453	!! END DO
454	!! END DO
455	!
456	! Do advection
457	IF (ln_vvl_adv_fct) THEN
458	CALL dom_vvl_adv_fct( kt, tilde_e3t_a, ztu, ztv )
459	!
460	ELSEIF (ln_vvl_adv_cn2) THEN
461	DO jk = 1, jpkm1
462	DO jj = 2, jpjm1
463	DO ji = fs_2, fs_jpim1 ! vector opt.
464	tilde_e3t_a(ji,jj,jk) = &
465	-( e2u(ji,jj)fse3u(ji,jj,jk) ztu(ji,jj,jk) - e2u(ji-1,jj )fse3u(ji-1,jj ,jk) ztu(ji-1,jj ,jk) &
466	+ e1v(ji,jj)fse3v(ji,jj,jk) ztv(ji,jj,jk) - e1v(ji ,jj-1)fse3v(ji ,jj-1,jk) ztv(ji ,jj-1,jk) ) &
467	/ ( e1t(ji,jj) * e2t(ji,jj) )
468	END DO
469	END DO
470	END DO
471	ENDIF
472	!
473	! Thickness anomaly diffusion:
474	! -----------------------------
475	ztu(:,:,:) = 0.0_wp
476	ztv(:,:,:) = 0.0_wp
477
478	ze3t(:,:,1) = 0.e0
479	DO jj = 1, jpj
480	DO ji = 1, jpi
481	DO jk = 2, jpk
482	ze3t(ji,jj,jk) = ze3t(ji,jj,jk-1) + tilde_e3t_b(ji,jj,jk-1) * tmask(ji,jj,jk-1)
483	END DO
484	END DO
485	END DO
486
487	IF ( ln_vvl_blp ) THEN ! Bilaplacian
488	DO jk = 1, jpkm1
489	DO jj = 1, jpjm1 ! First derivative (gradient)
490	DO ji = 1, fs_jpim1 ! vector opt.
491	! ztu(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) &
492	! & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) )
493	! ztv(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) &
494	! & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) )
495	ztu(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) &
496	& * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj ,jk) )
497	ztv(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) &
498	& * ( ze3t(ji,jj,jk) - ze3t(ji ,jj+1,jk) )
499	END DO
500	END DO
501
502	DO jj = 2, jpjm1 ! Second derivative (divergence) time the eddy diffusivity coefficient
503	DO ji = fs_2, fs_jpim1 ! vector opt.
504	zht(ji,jj) = rn_ahe3_blp * r1_e12t(ji,jj) * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) &
505	& + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )
506	END DO
507	END DO
508
509	#if defined key_bdy
510	DO ib_bdy=1, nb_bdy
511	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
512	ji = idx_bdy(ib_bdy)%nbi(ib,1)
513	jj = idx_bdy(ib_bdy)%nbj(ib,1)
514	zht(ji,jj) = 0._wp
515	END DO
516	END DO
517	#endif
518	CALL lbc_lnk( zht, 'T', 1. ) ! Lateral boundary conditions (unchanged sgn)
519
520	DO jj = 1, jpjm1 ! third derivative (gradient)
521	DO ji = 1, fs_jpim1 ! vector opt.
522	ztu(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) * ( zht(ji+1,jj ) - zht(ji,jj) )
523	ztv(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) * ( zht(ji ,jj+1) - zht(ji,jj) )
524	END DO
525	END DO
526	END DO
527	ENDIF
528
529	IF ( ln_vvl_lap ) THEN ! Laplacian
530	DO jk = 1, jpkm1 ! First derivative (gradient)
531	DO jj = 1, jpjm1
532	DO ji = 1, fs_jpim1 ! vector opt.
533	! zdu = rn_ahe3_lap * umask(ji,jj,jk) * re2u_e1u(ji,jj) &
534	! & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) )
535	! zdv = rn_ahe3_lap * vmask(ji,jj,jk) * re1v_e2v(ji,jj) &
536	! & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) )
537	zdu = rn_ahe3_lap * umask(ji,jj,jk) * re2u_e1u(ji,jj) &
538	& * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj ,jk) )
539	zdv = rn_ahe3_lap * vmask(ji,jj,jk) * re1v_e2v(ji,jj) &
540	& * ( ze3t(ji,jj,jk) - ze3t(ji ,jj+1,jk) )
541	ztu(ji,jj,jk) = ztu(ji,jj,jk) + zdu
542	ztv(ji,jj,jk) = ztv(ji,jj,jk) + zdv
543	END DO
544	END DO
545	END DO
546	ENDIF
547
548	! divergence of diffusive fluxes
549	DO jk = 1, jpkm1
550	DO jj = 2, jpjm1
551	DO ji = fs_2, fs_jpim1 ! vector opt.
552	! tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( ztu(ji-1,jj ,jk) - ztu(ji,jj,jk) &
553	! & + ztv(ji ,jj-1,jk) - ztv(ji,jj,jk) &
554	! & ) * r1_e12t(ji,jj)
555	un_td(ji,jj,jk) = un_td(ji,jj,jk) + ztu(ji,jj,jk+1) - ztu(ji,jj,jk )
556	vn_td(ji,jj,jk) = vn_td(ji,jj,jk) + ztv(ji,jj,jk+1) - ztv(ji,jj,jk )
557	tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( ztu(ji-1,jj ,jk+1) - ztu(ji,jj,jk+1) &
558	& +ztv(ji ,jj-1,jk+1) - ztv(ji,jj,jk+1) &
559	& -ztu(ji-1,jj ,jk ) + ztu(ji,jj,jk ) &
560	& -ztv(ji ,jj-1,jk ) + ztv(ji,jj,jk ) &
561	& ) * r1_e12t(ji,jj)
562	END DO
563	END DO
564	END DO
565
566
567	! un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:)
568	! vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:)
569
570	CALL lbc_lnk( un_td , 'U' , -1.)
571	CALL lbc_lnk( vn_td , 'V' , -1.)
572	!
573	CALL dom_vvl_ups_cor( kt, tilde_e3t_a, un_td, vn_td )
574
575	! IF ( ln_vvl_ztilde ) THEN
576	! ztu(:,:,:) = -un_lf(:,:,:)
577	! ztv(:,:,:) = -vn_lf(:,:,:)
578	! CALL dom_vvl_ups_cor( kt, tilde_e3t_a, ztu, ztv )
579	! ENDIF
580	!
581	! Remove "external thickness" tendency:
582	DO jk = 1, jpkm1
583	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * zhdiv(:,:)
584	END DO
585
586	! Leapfrog time stepping
587	! ~~~~~~~~~~~~~~~~~~~~~~
588	zramp = 1._wp
589	IF ((.NOT.ln_rstart).AND.ln_vvl_ramp) zramp = MIN(MAX( ((kt-nit000)rdt)/(rn_day_ramprday),0._wp),1._wp)
590
591	DO jk=1,jpkm1
592	tilde_e3t_a(:,:,jk) = tilde_e3t_b(:,:,jk) + z2dt * tmask(:,:,jk) * tilde_e3t_a(:,:,jk) &
593	& * tildemask(:,:) * zramp
594	END DO
595
596	! Ensure layer separation:
597	! ~~~~~~~~~~~~~~~~~~~~~~~~
598	IF ( ln_vvl_regrid ) CALL dom_vvl_regrid( kt )
599
600	! Boundary conditions:
601	! ~~~~~~~~~~~~~~~~~~~~
602	#if defined key_bdy
603	DO ib_bdy=1, nb_bdy
604	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
605	!! DO ib = 1, idx_bdy(ib_bdy)%nblen(1)
606	ji = idx_bdy(ib_bdy)%nbi(ib,1)
607	jj = idx_bdy(ib_bdy)%nbj(ib,1)
608	zwgt = idx_bdy(ib_bdy)%nbw(ib,1)
609	ip = bdytmask(ji+1,jj ) - bdytmask(ji-1,jj )
610	jp = bdytmask(ji ,jj+1) - bdytmask(ji ,jj-1)
611	DO jk = 1, jpkm1
612	tilde_e3t_a(ji,jj,jk) = 0.e0
613	!! tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) * (1._wp - zwgt)
614	!! tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji+ip,jj+jp,jk) * tmask(ji+ip,jj+jp,jk)
615	END DO
616	END DO
617	END DO
618	#endif
619	CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )
620
621	! Maximum deformation control
622	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
623	ze3t(:,:,jpk) = 0.0_wp
624	DO jk = 1, jpkm1
625	ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
626	END DO
627	z_tmax = MAXVAL( ze3t(:,:,:) )
628	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
629	z_tmin = MINVAL( ze3t(:,:,:) )
630	IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain
631	! - ML - test: for the moment, stop simulation for too large e3_t variations
632	IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN
633	IF( lk_mpp ) THEN
634	CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
635	CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
636	ELSE
637	ijk_max = MAXLOC( ze3t(:,:,:) )
638	ijk_max(1) = ijk_max(1) + nimpp - 1
639	ijk_max(2) = ijk_max(2) + njmpp - 1
640	ijk_min = MINLOC( ze3t(:,:,:) )
641	ijk_min(1) = ijk_min(1) + nimpp - 1
642	ijk_min(2) = ijk_min(2) + njmpp - 1
643	ENDIF
644	IF (lwp) THEN
645	WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
646	WRITE(numout, *) 'at i, j, k=', ijk_max
647	WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
648	WRITE(numout, *) 'at i, j, k=', ijk_min
649	CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
650	ENDIF
651	ENDIF
652	ENDIF
653
654	! IF( ln_vvl_ztilde ) THEN
655	! IF ( ncall==1 ) THEN
656	! zalpha = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
657	! DO jk = 1, jpkm1
658	! ztu(:,:,jk) = un(:,:,jk) * fse3u_n(:,:,jk) * e2u(:,:) + un_td(:,:,jk)
659	! ztv(:,:,jk) = vn(:,:,jk) * fse3v_n(:,:,jk) * e1v(:,:) + vn_td(:,:,jk)
660	! ze3t(:,:,jk) = -fse3t_n(:,:,jk) * zhdiv(:,:)
661	! END DO
662	! !
663	! un_lf(:,:,:) = un_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztu(:,:,:)
664	! vn_lf(:,:,:) = vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztv(:,:,:)
665	! hdivn_lf(:,:,:) = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:)
666	!! un_lf(:,:,:) = un_lf(:,:,:) * (1._wp - zalpha) + zalpha * un_lf2(:,:,:)
667	!! vn_lf(:,:,:) = vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * vn_lf2(:,:,:)
668	!! hdivn_lf(:,:,:) = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * hdivn_lf2(:,:,:)
669	! ENDIF
670	! ENDIF
671
672	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate !
673	! ! ---baroclinic part--------- !
674	DO jk = 1, jpkm1
675	fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk))
676	END DO
677	ENDIF
678
679	IF((ln_vvl_ztilde .OR. ln_vvl_layer).AND.(.NOT.ll_do_bclinic) ) THEN
680	zhdiv(:,:) = 0.
681	DO jk = 1, jpkm1
682	zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * (hdivn(:,:,jk) - hdivb(:,:,jk))
683	END DO
684	zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
685
686	IF( neuler == 0 .AND. kt == nit000 ) THEN
687	z2dt = rdt
688	ELSE
689	z2dt = 2.0_wp * rdt
690	ENDIF
691
692	DO jk = 1, jpkm1
693	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - z2dt * fse3t_n(:,:,jk) * &
694	& (hdivn(:,:,jk) - hdivb(:,:,jk) - zhdiv(:,:))
695	END DO
696	DO jk = 1, jpkm1
697	fse3t_a(:,:,jk) = fse3t_a(:,:,jk) - z2dt * fse3t_n(:,:,jk) * &
698	& (hdivn(:,:,jk) - hdivb(:,:,jk) - zhdiv(:,:))
699	END DO
700	ENDIF
701
702
703	IF( (ln_vvl_dbg .AND. ( ncall==2 )).AND.(ln_vvl_ztilde .OR. ln_vvl_layer) ) THEN ! - ML - test: control prints for debuging
704	!
705	zht(:,:) = 0.0_wp
706	DO jk = 1, jpkm1
707	zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
708	END DO
709	IF( lwp ) WRITE(numout, *) 'kt =', kt
710	IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
711	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ), mask = tmask(:,:,1) == 1.e0 )
712	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
713	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
714	END IF
715	!
716	z_tmin = MINVAL( fse3t_n(:,:,:), mask = tmask(:,:,:) == 1.e0 )
717	IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain
718	IF( lwp ) WRITE(numout, *) kt,' MINVAL(fse3t_n) =', z_tmin
719	IF ( z_tmin .LE. 0._wp ) THEN
720	IF( lk_mpp ) THEN
721	CALL mpp_minloc(fse3t_n(:,:,:), tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
722	ELSE
723	ijk_min = MINLOC( fse3t_n(:,:,:), mask = tmask(:,:,:) == 1.e0 )
724	ijk_min(1) = ijk_min(1) + nimpp - 1
725	ijk_min(2) = ijk_min(2) + njmpp - 1
726	ENDIF
727	IF (lwp) THEN
728	WRITE(numout, *) 'Negative scale factor, fse3t_n =', z_tmin
729	WRITE(numout, *) 'at i, j, k=', ijk_min
730	CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
731	ENDIF
732	ENDIF
733	!
734	z_tmin = MINVAL( fse3u_n(:,:,:), mask = umask(:,:,:) == 1.e0 )
735	IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain
736	IF( lwp ) WRITE(numout, *) kt,' MINVAL(fse3u_n) =', z_tmin
737	IF ( z_tmin .LE. 0._wp ) THEN
738	IF( lk_mpp ) THEN
739	CALL mpp_minloc(fse3u_n(:,:,:), umask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
740	ELSE
741	ijk_min = MINLOC( fse3u_n(:,:,:), mask = umask(:,:,:) == 1.e0 )
742	ijk_min(1) = ijk_min(1) + nimpp - 1
743	ijk_min(2) = ijk_min(2) + njmpp - 1
744	ENDIF
745	IF (lwp) THEN
746	WRITE(numout, *) 'Negative scale factor, fse3u_n =', z_tmin
747	WRITE(numout, *) 'at i, j, k=', ijk_min
748	CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
749	ENDIF
750	ENDIF
751	!
752	z_tmin = MINVAL( fse3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0 )
753	IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain
754	IF( lwp ) WRITE(numout, *) kt,' MINVAL(fse3v_n) =', z_tmin
755	IF ( z_tmin .LE. 0._wp ) THEN
756	IF( lk_mpp ) THEN
757	CALL mpp_minloc(fse3v_n(:,:,:), vmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
758	ELSE
759	ijk_min = MINLOC( fse3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0 )
760	ijk_min(1) = ijk_min(1) + nimpp - 1
761	ijk_min(2) = ijk_min(2) + njmpp - 1
762	ENDIF
763	IF (lwp) THEN
764	WRITE(numout, *) 'Negative scale factor, fse3v_n =', z_tmin
765	WRITE(numout, *) 'at i, j, k=', ijk_min
766	CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
767	ENDIF
768	ENDIF
769	!
770	z_tmin = MINVAL( fse3f_n(:,:,:), mask = fmask(:,:,:) == 1.e0 )
771	IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain
772	IF( lwp ) WRITE(numout, *) kt,' MINVAL(fse3f_n) =', z_tmin
773	IF ( z_tmin .LE. 0._wp ) THEN
774	IF( lk_mpp ) THEN
775	CALL mpp_minloc(fse3f_n(:,:,:), fmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
776	ELSE
777	ijk_min = MINLOC( fse3f_n(:,:,:), mask = fmask(:,:,:) == 1.e0 )
778	ijk_min(1) = ijk_min(1) + nimpp - 1
779	ijk_min(2) = ijk_min(2) + njmpp - 1
780	ENDIF
781	IF (lwp) THEN
782	WRITE(numout, *) 'Negative scale factor, fse3f_n =', z_tmin
783	WRITE(numout, *) 'at i, j, k=', ijk_min
784	CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor')
785	ENDIF
786	ENDIF
787	!
788	zht(:,:) = 0.0_wp
789	DO jk = 1, jpkm1
790	zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
791	END DO
792	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
793	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
794	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn -SUM(fse3t_n))) =', z_tmax
795	!
796	zht(:,:) = 0.0_wp
797	DO jk = 1, jpkm1
798	zht(:,:) = zht(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
799	END DO
800	zwu(:,:) = 0._wp
801	DO jj = 1, jpjm1
802	DO ji = 1, fs_jpim1 ! vector opt.
803	zwu(ji,jj) = 0.5_wp * umask(ji,jj,1) * r1_e12u(ji,jj) &
804	& * ( e12t(ji,jj) * sshn(ji,jj) + e12t(ji+1,jj) * sshn(ji+1,jj) )
805	END DO
806	END DO
807	CALL lbc_lnk( zwu(:,:), 'U', 1._wp )
808	z_tmax = MAXVAL( umask(:,:,1) * umask_i(:,:) * ABS( hu_0(:,:) + zwu(:,:) - zht(:,:) ) )
809	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
810	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hu_0+sshu_n-SUM(fse3u_n))) =', z_tmax
811	!
812	zht(:,:) = 0.0_wp
813	DO jk = 1, jpkm1
814	zht(:,:) = zht(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
815	END DO
816	zwv(:,:) = 0._wp
817	DO jj = 1, jpjm1
818	DO ji = 1, fs_jpim1 ! vector opt.
819	zwv(ji,jj) = 0.5_wp * vmask(ji,jj,1) * r1_e12v(ji,jj) &
820	& * ( e12t(ji,jj) * sshn(ji,jj) + e12t(ji,jj+1) * sshn(ji,jj+1) )
821	END DO
822	END DO
823	CALL lbc_lnk( zwv(:,:), 'V', 1._wp )
824	z_tmax = MAXVAL( vmask(:,:,1) * vmask_i(:,:) * ABS( hv_0(:,:) + zwv(:,:) - zht(:,:) ) )
825	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
826	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hv_0+sshv_n-SUM(fse3v_n))) =', z_tmax
827	!
828	zht(:,:) = 0.0_wp
829	DO jk = 1, jpkm1
830	DO jj = 1, jpjm1
831	zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk)*umask(:,jj+1,jk)
832	END DO
833	END DO
834	zwu(:,:) = 0._wp
835	DO jj = 1, jpjm1
836	DO ji = 1, fs_jpim1 ! vector opt.
837	zwu(ji,jj) = 0.25_wp * umask(ji,jj,1) * umask(ji,jj+1,1) * r1_e12f(ji,jj) &
838	& * ( e12t(ji ,jj) * sshn(ji ,jj) + e12t(ji ,jj+1) * sshn(ji ,jj+1) &
839	& + e12t(ji+1,jj) * sshn(ji+1,jj) + e12t(ji+1,jj+1) * sshn(ji+1,jj+1) )
840	END DO
841	END DO
842	CALL lbc_lnk( zht(:,:), 'F', 1._wp )
843	CALL lbc_lnk( zwu(:,:), 'F', 1._wp )
844	z_tmax = MAXVAL( fmask(:,:,1) * fmask_i(:,:) * ABS( hf_0(:,:) + zwu(:,:) - zht(:,:) ) )
845	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
846	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hf_0+sshf_n-SUM(fse3f_n))) =', z_tmax
847	!
848	END IF
849
850	! *********************************** !
851	! After scale factors at u- v- points !
852	! *********************************** !
853
854	CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )
855	CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
856
857	! *********************************** !
858	! After depths at u- v points !
859	! *********************************** !
860
861	hu_a(:,:) = 0._wp ! Ocean depth at U-points
862	hv_a(:,:) = 0._wp ! Ocean depth at V-points
863	DO jk = 1, jpkm1
864	hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
865	hv_a(:,:) = hv_a(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
866	END DO
867	! ! Inverse of the local depth
868	hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
869	hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
870
871	CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
872	CALL wrk_dealloc( jpi, jpj, jpk, ze3t, ztu, ztv )
873
874	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_nxt')
875
876	END SUBROUTINE dom_vvl_sf_nxt
877
878
879	SUBROUTINE dom_vvl_sf_swp( kt )
880	!!----------------------------------------------------------------------
881	!! * ROUTINE dom_vvl_sf_swp *
882	!!
883	!! ** Purpose : compute time filter and swap of scale factors
884	!! compute all depths and related variables for next time step
885	!! write outputs and restart file
886	!!
887	!! ** Method : - swap of e3t with trick for volume/tracer conservation
888	!! - reconstruct scale factor at other grid points (interpolate)
889	!! - recompute depths and water height fields
890	!!
891	!! ** Action : - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step
892	!! - Recompute:
893	!! fse3(u/v)_b
894	!! fse3w_n
895	!! fse3(u/v)w_b
896	!! fse3(u/v)w_n
897	!! fsdept_n, fsdepw_n and fsde3w_n
898	!! h(u/v) and h(u/v)r
899	!!
900	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
901	!! Leclair, M., and G. Madec, 2011, Ocean Modelling.
902	!!----------------------------------------------------------------------
903	!! * Arguments
904	INTEGER, INTENT( in ) :: kt ! time step
905	!! * Local declarations
906	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwrk
907	INTEGER :: jk ! dummy loop indices
908	!!----------------------------------------------------------------------
909
910	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp')
911	!
912	IF( kt == nit000 ) THEN
913	IF(lwp) WRITE(numout,*)
914	IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
915	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step'
916	ENDIF
917	!
918	! Time filter and swap of scale factors
919	! =====================================
920	! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
921	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
922	IF( neuler == 0 .AND. kt == nit000 ) THEN
923	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
924	ELSE
925	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &
926	& + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
927	ENDIF
928	tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
929	ENDIF
930
931	fsdept_b(:,:,:) = fsdept_n(:,:,:)
932	fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
933
934	fse3t_n(:,:,:) = fse3t_a(:,:,:)
935	fse3u_n(:,:,:) = fse3u_a(:,:,:)
936	fse3v_n(:,:,:) = fse3v_a(:,:,:)
937
938	! Compute all missing vertical scale factor and depths
939	! ====================================================
940	! Horizontal scale factor interpolations
941	! --------------------------------------
942	! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
943	! - JC - hu_b, hv_b, hur_b, hvr_b also
944	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n (:,:,:), 'F' )
945	! Vertical scale factor interpolations
946	! ------------------------------------
947	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )
948	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
949	CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
950	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' )
951	CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
952	CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
953	! t- and w- points depth
954	! ----------------------
955	fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
956	fsdepw_n(:,:,1) = 0.0_wp
957	fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
958	DO jk = 2, jpk
959	fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
960	fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
961	fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:)
962	END DO
963	! Local depth and Inverse of the local depth of the water column at u- and v- points
964	! ----------------------------------------------------------------------------------
965	hu (:,:) = hu_a (:,:)
966	hv (:,:) = hv_a (:,:)
967
968	! Inverse of the local depth
969	hur(:,:) = hur_a(:,:)
970	hvr(:,:) = hvr_a(:,:)
971
972	! Local depth of the water column at t- points
973	! --------------------------------------------
974	ht(:,:) = 0.
975	DO jk = 1, jpkm1
976	ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
977	END DO
978
979	! Write additional diagnostics
980	! ============================
981	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) CALL dom_vvl_dia( kt)
982
983	! write restart file
984	! ==================
985	IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
986	!
987	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp')
988
989	END SUBROUTINE dom_vvl_sf_swp
990
991	SUBROUTINE dom_vvl_dia( kt )
992	!!----------------------------------------------------------------------
993	!! * ROUTINE dom_vvl_dia *
994	!!
995	!! ** Purpose : Output some diagnostics in ztilde/zlayer cases
996	!!
997	!!----------------------------------------------------------------------
998	!! * Arguments
999	INTEGER, INTENT( in ) :: kt ! time step
1000	!! * Local declarations
1001	INTEGER :: ji,jj,jk ! dummy loop indices
1002	REAL(wp) :: zufim1, zufi, zvfjm1, zvfj, ztmp1, z2dt
1003	REAL(wp), DIMENSION(4) :: zr1
1004	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwdw, zout
1005	!!----------------------------------------------------------------------
1006	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_dia')
1007	!
1008	CALL wrk_alloc( jpi, jpj, jpk, zwdw, zout )
1009	!
1010	! Compute internal interfaces depths:
1011	!------------------------------------
1012	IF ( iom_use("dh_tilde").OR.iom_use("depw_tilde").OR.iom_use("stiff_tilde")) THEN
1013	zwdw(:,:,1) = 0.e0
1014	DO jj = 1, jpj
1015	DO ji = 1, jpi
1016	DO jk = 2, jpkm1
1017	zwdw(ji,jj,jk) = zwdw(ji,jj,jk-1) + &
1018	& (tilde_e3t_n(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
1019	END DO
1020	END DO
1021	END DO
1022	ENDIF
1023	!
1024	! Output interface depth anomaly:
1025	! -------------------------------
1026	IF ( iom_use("depw_tilde") ) CALL iom_put( "depw_tilde", (zwdw(:,:,:)-gdepw_0(:,:,:))*tmask(:,:,:) )
1027	!
1028	! Output grid stiffness (T-points):
1029	! ---------------------------------
1030	IF ( iom_use("stiff_tilde" ) ) THEN
1031	zr1(:) = 0.e0
1032	zout(:,:,jpk) = 0.e0
1033	DO ji = 2, jpim1
1034	DO jj = 2, jpjm1
1035	DO jk = 1, jpkm1
1036	zr1(1) = umask(ji-1,jj ,jk) *abs( (zwdw(ji ,jj ,jk )-zwdw(ji-1,jj ,jk ) &
1037	& +zwdw(ji ,jj ,jk+1)-zwdw(ji-1,jj ,jk+1)) &
1038	& /(zwdw(ji ,jj ,jk )+zwdw(ji-1,jj ,jk ) &
1039	& -zwdw(ji ,jj ,jk+1)-zwdw(ji-1,jj ,jk+1) + rsmall) )
1040	zr1(2) = umask(ji ,jj ,jk) *abs( (zwdw(ji+1,jj ,jk )-zwdw(ji ,jj ,jk ) &
1041	& +zwdw(ji+1,jj ,jk+1)-zwdw(ji ,jj ,jk+1)) &
1042	& /(zwdw(ji+1,jj ,jk )+zwdw(ji ,jj ,jk ) &
1043	& -zwdw(ji+1,jj ,jk+1)-zwdw(ji ,jj ,jk+1) + rsmall) )
1044	zr1(3) = vmask(ji ,jj ,jk) *abs( (zwdw(ji ,jj+1,jk )-zwdw(ji ,jj ,jk ) &
1045	& +zwdw(ji ,jj+1,jk+1)-zwdw(ji ,jj ,jk+1)) &
1046	& /(zwdw(ji ,jj+1,jk )+zwdw(ji ,jj ,jk ) &
1047	& -zwdw(ji ,jj+1,jk+1)-zwdw(ji ,jj ,jk+1) + rsmall) )
1048	zr1(4) = vmask(ji ,jj-1,jk) *abs( (zwdw(ji ,jj ,jk )-zwdw(ji ,jj-1,jk ) &
1049	& +zwdw(ji ,jj ,jk+1)-zwdw(ji ,jj-1,jk+1)) &
1050	& /(zwdw(ji ,jj ,jk )+zwdw(ji ,jj-1,jk ) &
1051	& -zwdw(ji, jj ,jk+1)-zwdw(ji ,jj-1,jk+1) + rsmall) )
1052	zout(ji,jj,jk) = MAXVAL(zr1(1:4))
1053	END DO
1054	END DO
1055	END DO
1056
1057	CALL lbc_lnk( zout, 'T', 1. )
1058	CALL iom_put( "stiff_tilde", zout(:,:,:) )
1059	END IF
1060	! Output Horizontal Laplacian of interfaces depths (W-points):
1061	! ------------------------------------------------------------
1062	IF ( iom_use("dh_tilde") ) THEN
1063	!
1064	zout(:,:,1 )=0._wp
1065	zout(:,:,jpk)=0._wp
1066	DO jk = 2, jpkm1
1067	DO jj = 1, jpjm1
1068	DO ji = 1, fs_jpim1 ! vector opt.
1069	ua(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) &
1070	& * ( zwdw(ji,jj,jk) - zwdw(ji+1,jj ,jk) )
1071	va(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) &
1072	& * ( zwdw(ji,jj,jk) - zwdw(ji ,jj+1,jk) )
1073	END DO
1074	END DO
1075	END DO
1076
1077	DO jk = 2, jpkm1
1078	DO jj = 2, jpjm1
1079	DO ji = fs_2, fs_jpim1 ! vector opt.
1080	ztmp1 = ( (ua(ji-1,jj ,jk) - ua(ji,jj,jk)) &
1081	& + (va(ji ,jj-1,jk) - va(ji,jj,jk)) ) * SQRT(r1_e12t(ji,jj))
1082	zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk)
1083	END DO
1084	END DO
1085	END DO
1086	! Mask open boundaries:
1087	#if defined key_bdy
1088	IF (lk_bdy) THEN
1089	DO jk = 1, jpkm1
1090	zout(:,:,jk) = zout(:,:,jk) * bdytmask(:,:)
1091	END DO
1092	ENDIF
1093	#endif
1094	CALL lbc_lnk( zout(:,:,:), 'T', 1. )
1095	!
1096	CALL iom_put( "dh_tilde", zout(:,:,:) )
1097	ENDIF
1098	!
1099	! Output vertical Laplacian of interfaces depths (W-points):
1100	! ----------------------------------------------------------
1101	IF ( iom_use("dz_tilde" ) ) THEN
1102	zout(:,:,1 ) = 0._wp
1103	zout(:,:,jpk) = 0._wp
1104	DO jk=2,jpkm1
1105	zout(:,:,jk) = 2._wp*ABS(tilde_e3t_n(:,:,jk)+e3t_0(:,:,jk)-tilde_e3t_n(:,:,jk-1)-e3t_0(:,:,jk-1)) &
1106	& /(tilde_e3t_n(:,:,jk)+e3t_0(:,:,jk)+tilde_e3t_n(:,:,jk-1)+e3t_0(:,:,jk-1)) &
1107	& * tmask(:,:,jk)
1108	END DO
1109	CALL iom_put( "dz_tilde", zout(:,:,:) )
1110	END IF
1111	!
1112	!
1113	! Output low pass U-velocity:
1114	! ---------------------------
1115	IF ( iom_use("un_lf_tilde" ).AND.ln_vvl_ztilde ) THEN
1116	zout(:,:,jpk) = 0.e0
1117	DO jk=1,jpkm1
1118	zout(:,:,jk) = un_lf(:,:,jk)/fse3u_n(:,:,jk)*r1_e2u(:,:)
1119	END DO
1120	CALL iom_put( "un_lf_tilde", zout(:,:,:) )
1121	END IF
1122	!
1123	! Output low pass V-velocity:
1124	! ---------------------------
1125	IF ( iom_use("vn_lf_tilde" ).AND.ln_vvl_ztilde ) THEN
1126	zout(:,:,jpk) = 0.e0
1127	DO jk=1,jpkm1
1128	zout(:,:,jk) = vn_lf(:,:,jk)/fse3v_n(:,:,jk)*r1_e1v(:,:)
1129	END DO
1130	CALL iom_put( "vn_lf_tilde", zout(:,:,:) )
1131	END IF
1132	!
1133	CALL wrk_dealloc( jpi, jpj, jpk, zwdw, zout )
1134	!
1135	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_dia')
1136	!
1137	END SUBROUTINE dom_vvl_dia
1138
1139	SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
1140	!!---------------------------------------------------------------------
1141	!! * ROUTINE dom_vvl__interpol *
1142	!!
1143	!! ** Purpose : interpolate scale factors from one grid point to another
1144	!!
1145	!! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0)
1146	!! - horizontal interpolation: grid cell surface averaging
1147	!! - vertical interpolation: simple averaging
1148	!!----------------------------------------------------------------------
1149	!! * Arguments
1150	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated
1151	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3
1152	CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors
1153	! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'
1154	!! * Local declarations
1155	INTEGER :: nmet ! horizontal interpolation method
1156	INTEGER :: ji, jj, jk ! dummy loop indices
1157	INTEGER :: jkbot ! bottom level
1158	LOGICAL :: l_is_orca ! local logical
1159	REAL(wp) :: zmin, zdo, zup, ztap
1160	REAL(wp), POINTER, DIMENSION(:,:) :: zs ! Surface interface depth anomaly
1161	REAL(wp), POINTER, DIMENSION(:,:,:) :: zw ! Interface depth anomaly
1162	!!----------------------------------------------------------------------
1163	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol')
1164	!
1165	l_is_orca = .FALSE.
1166	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE. ! ORCA R2 configuration - will need to correct some locations
1167
1168
1169	nmet=1
1170	! nmet=0 ! Original method (Surely wrong)
1171	! nmet=1 ! Interface interpolation
1172	! nmet=2 ! Internal interfaces interpolation only, spread barotropic increment
1173	ztap = 0.1_wp ! Minimum fraction of T-point thickness at cell interfaces
1174
1175	IF ( (nmet==1).OR.(nmet==2) ) THEN
1176	SELECT CASE ( pout )
1177	!
1178	CASE( 'U', 'V', 'F' )
1179	! Compute interface depth anomaly at T-points
1180	CALL wrk_alloc( jpi, jpj, jpk, zw )
1181	CALL wrk_alloc( jpi, jpj, zs )
1182	!
1183	zw(:,:,:) = 0._wp
1184	!
1185	DO jj = 1, jpj
1186	DO ji = 1, jpi
1187	jkbot = mbkt(ji,jj)
1188	DO jk=jkbot,1,-1
1189	zw(ji,jj,jk) = zw(ji,jj,jk+1) - pe3_in(ji,jj,jk) + e3t_0(ji,jj,jk)
1190	END DO
1191	END DO
1192	END DO
1193	!
1194	!
1195	IF (nmet==2) THEN ! Consider "internal" interfaces only
1196	zs(:,:) = - zw(:,:,1) ! Save surface anomaly (ssh)
1197	!
1198	zw(:,:,:) = 0._wp
1199	DO jj = 1, jpj
1200	DO ji = 1, jpi
1201	jkbot = mbkt(ji,jj)
1202	DO jk=jkbot,1,-1
1203	zw(ji,jj,jk) = zw(ji,jj,jk+1) - ( pe3_in(ji,jj,jk) &
1204	& * ht_0(ji,jj) / (ht_0(ji,jj) + zs(ji,jj) + 1._wp - tmask(ji,jj,1)) &
1205	& - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
1206	END DO
1207	END DO
1208	END DO
1209	ENDIF
1210	!
1211	END SELECT
1212	END IF
1213
1214	pe3_out(:,:,:) = 0._wp
1215
1216	SELECT CASE ( pout )
1217	! ! ------------------------------------- !
1218	CASE( 'U' ) ! interpolation from T-point to U-point !
1219	! ! ------------------------------------- !
1220	! horizontal surface weighted interpolation
1221	IF (nmet==0) THEN
1222	DO jk = 1, jpk
1223	DO jj = 1, jpjm1
1224	DO ji = 1, fs_jpim1 ! vector opt.
1225	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) &
1226	& * ( e12t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1227	& + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
1228	END DO
1229	END DO
1230	END DO
1231	ENDIF
1232	!
1233	IF ( (nmet==1).OR.(nmet==2) ) THEN
1234	DO jj = 1, jpjm1
1235	DO ji = 1, fs_jpim1 ! vector opt.
1236	! Correction at last level:
1237	jkbot = mbku(ji,jj)
1238	pe3_out(ji,jj,jkbot) = -0.5_wp * umask(ji,jj,jkbot) * r1_e12u(ji,jj) &
1239	& * ( e12t(ji ,jj) * zw(ji ,jj,jkbot) &
1240	& + e12t(ji+1,jj) * zw(ji+1,jj,jkbot) )
1241	!
1242	! If there is a step, taper bottom interface:
1243	IF (hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ) THEN
1244	IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN
1245	! zmin = ztap * pe3_in(ji+1,jj,jkbot)
1246	zmin = ztap * (-zw(ji+1,jj,jkbot)+e3t_0(ji+1,jj,jkbot))
1247	ELSE
1248	! zmin = ztap * pe3_in(ji ,jj,jkbot)
1249	zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot))
1250	ENDIF
1251	zmin = -e3u_0(ji,jj,jkbot) + zmin
1252	pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin)
1253	ENDIF
1254	!
1255	zdo = -pe3_out(ji,jj,jkbot)
1256	DO jk=jkbot-1,1,-1
1257	zup = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji ,jj) &
1258	& ( e12t(ji ,jj) zw(ji ,jj,jk) &
1259	& +e12t(ji+1,jj) * zw(ji+1,jj,jk) )
1260	pe3_out(ji,jj,jk) = zdo - zup
1261	zdo = zdo - pe3_out(ji,jj,jk)
1262	END DO
1263	END DO
1264	END DO
1265	END IF
1266	!
1267	IF (nmet==2) THEN ! Spread sea level anomaly
1268	DO jj = 1, jpjm1
1269	DO ji = 1, fs_jpim1 ! vector opt.
1270	DO jk=1,jpk
1271	pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) &
1272	& + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) ) &
1273	& / ( hu_0(ji,jj) + 1._wp - umask_i(ji,jj) ) &
1274	& * 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) &
1275	& * ( e12t(ji,jj)zs(ji,jj) + e12t(ji+1,jj)zs(ji+1,jj) )
1276	END DO
1277	END DO
1278	END DO
1279	!
1280	ENDIF
1281	!
1282	IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1283	! boundary conditions
1284	CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
1285	pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
1286	!
1287	IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs )
1288	IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw )
1289	!
1290	! ! ------------------------------------- !
1291	CASE( 'V' ) ! interpolation from T-point to V-point !
1292	! ! ------------------------------------- !
1293	! horizontal surface weighted interpolation
1294	IF (nmet==0) THEN
1295	DO jk = 1, jpk
1296	DO jj = 1, jpjm1
1297	DO ji = 1, fs_jpim1 ! vector opt.
1298	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) &
1299	& * ( e12t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1300	& + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
1301	END DO
1302	END DO
1303	END DO
1304	ENDIF
1305	!
1306	IF ( (nmet==1).OR.(nmet==2) ) THEN
1307	DO jj = 1, jpjm1
1308	DO ji = 1, fs_jpim1 ! vector opt.
1309	! Correction at last level:
1310	jkbot = mbkv(ji,jj)
1311	pe3_out(ji,jj,jkbot) = -0.5_wp * vmask(ji,jj,jkbot) * r1_e12v(ji,jj) &
1312	& * ( e12t(ji,jj ) * zw(ji,jj ,jkbot) &
1313	& + e12t(ji,jj+1) * zw(ji,jj+1,jkbot) )
1314	!
1315	! If there is a step, taper bottom interface:
1316	IF (hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ) THEN
1317	IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN
1318	! zmin = ztap * pe3_in(ji,jj+1,jkbot)
1319	zmin = ztap * (-zw(ji,jj+1,jkbot)+e3t_0(ji,jj+1,jkbot))
1320	ELSE
1321	! zmin = ztap * pe3_in(ji,jj ,jkbot)
1322	zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot))
1323	ENDIF
1324	zmin = -e3v_0(ji,jj,jkbot) + zmin
1325	pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin)
1326	ENDIF
1327	!
1328	zdo = -pe3_out(ji,jj,jkbot)
1329	DO jk=jkbot-1,1,-1
1330	zup = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj ) &
1331	& * ( e12t(ji,jj ) * zw(ji,jj ,jk) &
1332	& +e12t(ji,jj+1) * zw(ji,jj+1,jk) )
1333	!
1334	pe3_out(ji,jj,jk) = zdo - zup
1335	zdo = zdo - pe3_out(ji,jj,jk)
1336	END DO
1337	END DO
1338	END DO
1339	END IF
1340	!
1341	IF (nmet==2) THEN ! Spread sea level anomaly
1342	!
1343	DO jj = 1, jpjm1
1344	DO ji = 1, fs_jpim1 ! vector opt.
1345	DO jk=1,jpk
1346	pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) &
1347	& + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) ) &
1348	& / ( hv_0(ji,jj) + 1._wp - vmask_i(ji,jj) ) &
1349	& * 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) &
1350	& * ( e12t(ji,jj)zs(ji,jj) + e12t(ji,jj+1)zs(ji,jj+1) )
1351	END DO
1352	END DO
1353	END DO
1354	!
1355	ENDIF
1356	!
1357	IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1358	! boundary conditions
1359	CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
1360	pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
1361	!
1362	IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs )
1363	IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw )
1364	!
1365	! ! ------------------------------------- !
1366	CASE( 'F' ) ! interpolation from T-point to F-point !
1367	! ! ------------------------------------- !
1368	! horizontal surface weighted interpolation
1369	IF (nmet==0) THEN
1370	DO jk=1,jpk
1371	DO jj = 1, jpjm1
1372	DO ji = 1, fs_jpim1 ! vector opt.
1373	pe3_out(ji,jj,jk) = 0.25_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) &
1374	& * ( e12t(ji ,jj ) * ( pe3_in(ji ,jj ,jk)-e3t_0(ji ,jj ,jk) ) &
1375	& + e12t(ji ,jj+1) * ( pe3_in(ji ,jj+1,jk)-e3t_0(ji ,jj+1,jk) ) &
1376	& + e12t(ji+1,jj ) * ( pe3_in(ji+1,jj ,jk)-e3t_0(ji+1,jj ,jk) ) &
1377	& + e12t(ji+1,jj+1) * ( pe3_in(ji+1,jj+1,jk)-e3t_0(ji+1,jj+1,jk) ))
1378	END DO
1379	END DO
1380	END DO
1381	ENDIF
1382	!
1383	IF ( (nmet==1).OR.(nmet==2) ) THEN
1384	DO jj = 1, jpjm1
1385	DO ji = 1, fs_jpim1 ! vector opt.
1386	! bottom correction:
1387	jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1))
1388	pe3_out(ji,jj,jkbot) = -0.25_wp * umask(ji,jj,jkbot) * umask(ji,jj+1,jkbot) * r1_e12f(ji,jj) &
1389	& * ( e12t(ji ,jj ) * zw(ji ,jj ,jkbot) &
1390	& + e12t(ji+1,jj ) * zw(ji+1,jj ,jkbot) &
1391	& + e12t(ji ,jj+1) * zw(ji ,jj+1,jkbot) &
1392	& + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jkbot) )
1393	!
1394	! If there is a step, taper bottom interface:
1395	IF (hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj ) + hu_0(ji,jj+1) ) ) THEN
1396	IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN
1397	IF ( ht_0(ji+1,jj+1) < ht_0(ji ,jj+1) ) THEN
1398	! zmin = ztap * pe3_in(ji+1,jj+1,jkbot)
1399	zmin = ztap * (-zw(ji+1,jj+1,jkbot)+e3t_0(ji+1,jj+1,jkbot))
1400	ELSE
1401	! zmin = ztap * pe3_in(ji ,jj+1,jkbot)
1402	zmin = ztap * (-zw(ji,jj+1,jkbot)+e3t_0(ji,jj+1,jkbot))
1403	ENDIF
1404	ELSE
1405	IF ( ht_0(ji+1,jj ) < ht_0(ji ,jj ) ) THEN
1406	! zmin = ztap * pe3_in(ji+1,jj ,jkbot)
1407	zmin = ztap * (-zw(ji+1,jj,jkbot)+e3t_0(ji+1,jj,jkbot))
1408	ELSE
1409	! zmin = ztap * pe3_in(ji ,jj ,jkbot)
1410	zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot))
1411	ENDIF
1412	ENDIF
1413	zmin = -e3f_0(ji,jj,jkbot) + zmin
1414	pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin)
1415	ENDIF
1416	!
1417	zdo = -pe3_out(ji,jj,jkbot)
1418	DO jk=jkbot-1,1,-1
1419	zup = 0.25_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) &
1420	& * ( e12t(ji ,jj ) * zw(ji ,jj ,jk) &
1421	& + e12t(ji+1,jj ) * zw(ji+1,jj ,jk) &
1422	& + e12t(ji ,jj+1) * zw(ji ,jj+1,jk) &
1423	& + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jk) )
1424	pe3_out(ji,jj,jk) = zdo - zup
1425	!
1426	zdo = zdo - pe3_out(ji,jj,jk)
1427	END DO
1428	!
1429	END DO
1430	END DO
1431	ENDIF
1432	!
1433	IF (nmet==2) THEN ! Spread sea level anomaly
1434	!
1435	DO jj = 1, jpjm1
1436	DO ji = 1, fs_jpim1 ! vector opt.
1437	DO jk=1,jpk
1438	pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) &
1439	& + ( pe3_out(ji,jj,jk) + e3f_0(ji,jj,jk) ) &
1440	& / ( hf_0(ji,jj) + 1._wp - umask_i(ji,jj)*umask_i(ji,jj+1) ) &
1441	& * 0.25_wp * umask(ji,jj,jk)umask(ji,jj+1,jk)r1_e12f(ji,jj) &
1442	& * ( e12t(ji ,jj)zs(ji ,jj) + e12t(ji ,jj+1)zs(ji ,jj+1) &
1443	& +e12t(ji+1,jj)zs(ji+1,jj) + e12t(ji+1,jj+1)zs(ji+1,jj+1) )
1444	END DO
1445	END DO
1446	END DO
1447	END IF
1448	!
1449	IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1450	! boundary conditions
1451	CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
1452	pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
1453	!
1454	IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs )
1455	IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw )
1456	!
1457	! ! ------------------------------------- !
1458	CASE( 'W' ) ! interpolation from T-point to W-point !
1459	! ! ------------------------------------- !
1460	! vertical simple interpolation
1461	pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
1462	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
1463	DO jk = 2, jpk
1464	pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) &
1465	& + 0.5_wp * tmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) )
1466	END DO
1467	! ! -------------------------------------- !
1468	CASE( 'UW' ) ! interpolation from U-point to UW-point !
1469	! ! -------------------------------------- !
1470	! vertical simple interpolation
1471	pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
1472	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
1473	DO jk = 2, jpk
1474	pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) &
1475	& + 0.5_wp * umask(:,:,jk) * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) )
1476	END DO
1477	! ! -------------------------------------- !
1478	CASE( 'VW' ) ! interpolation from V-point to VW-point !
1479	! ! -------------------------------------- !
1480	! vertical simple interpolation
1481	pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
1482	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
1483	DO jk = 2, jpk
1484	pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) &
1485	& + 0.5_wp * vmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) )
1486	END DO
1487	END SELECT
1488	!
1489
1490	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol')
1491
1492	END SUBROUTINE dom_vvl_interpol
1493
1494
1495	SUBROUTINE dom_vvl_rst( kt, cdrw )
1496	!!---------------------------------------------------------------------
1497	!! * ROUTINE dom_vvl_rst *
1498	!!
1499	!! ** Purpose : Read or write VVL file in restart file
1500	!!
1501	!! ** Method : use of IOM library
1502	!! if the restart does not contain vertical scale factors,
1503	!! they are set to the _0 values
1504	!! if the restart does not contain vertical scale factors increments (z_tilde),
1505	!! they are set to 0.
1506	!!----------------------------------------------------------------------
1507	!! * Arguments
1508	INTEGER , INTENT(in) :: kt ! ocean time-step
1509	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1510	!! * Local declarations
1511	INTEGER :: jk
1512	INTEGER, DIMENSION(7) :: id ! local integers
1513	REAL(wp), POINTER, DIMENSION(:,:) :: zhdiv
1514	!!----------------------------------------------------------------------
1515	!
1516	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_rst')
1517	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
1518	! ! ===============
1519	IF( ln_rstart ) THEN !* Read the restart file
1520	CALL rst_read_open ! open the restart file if necessary
1521	CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn )
1522	!
1523	id(1) = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
1524	id(2) = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
1525	id(3) = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
1526	id(4) = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
1527	id(5) = iom_varid( numror, 'un_lf' , ldstop = .FALSE. )
1528	id(6) = iom_varid( numror, 'vn_lf' , ldstop = .FALSE. )
1529	id(7) = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. )
1530	! ! --------- !
1531	! ! all cases !
1532	! ! --------- !
1533	IF( MIN( id(1), id(2) ) > 0 ) THEN ! all required arrays exist
1534	CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
1535	CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
1536	! needed to restart if land processor not computed
1537	IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'
1538	WHERE ( tmask(:,:,:) == 0.0_wp )
1539	fse3t_n(:,:,:) = e3t_0(:,:,:)
1540	fse3t_b(:,:,:) = e3t_0(:,:,:)
1541	END WHERE
1542	IF( neuler == 0 ) THEN
1543	fse3t_b(:,:,:) = fse3t_n(:,:,:)
1544	ENDIF
1545	ELSE IF( id(1) > 0 ) THEN
1546	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
1547	IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
1548	IF(lwp) write(numout,*) 'neuler is forced to 0'
1549	CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
1550	fse3t_b(:,:,:) = fse3t_n(:,:,:)
1551	neuler = 0
1552	ELSE IF( id(2) > 0 ) THEN
1553	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
1554	IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
1555	IF(lwp) write(numout,*) 'neuler is forced to 0'
1556	CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
1557	fse3t_b(:,:,:) = fse3t_n(:,:,:)
1558	neuler = 0
1559	ELSE
1560	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'
1561	IF(lwp) write(numout,*) 'Compute scale factor from sshn'
1562	IF(lwp) write(numout,*) 'neuler is forced to 0'
1563	DO jk=1,jpk
1564	fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
1565	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
1566	& + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
1567	END DO
1568	fse3t_b(:,:,:) = fse3t_n(:,:,:)
1569	neuler = 0
1570	ENDIF
1571	! ! ----------- !
1572	IF( ln_vvl_zstar ) THEN ! z_star case !
1573	! ! ----------- !
1574	IF( MIN( id(3), id(4) ) > 0 ) THEN
1575	CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
1576	ENDIF
1577	! ! ----------------------- !
1578	ELSE ! z_tilde and layer cases !
1579	! ! ----------------------- !
1580	IF( MIN( id(3), id(4) ) > 0 ) THEN ! all required arrays exist
1581	CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
1582	CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
1583	ELSE ! one at least array is missing
1584	tilde_e3t_b(:,:,:) = 0.0_wp
1585	tilde_e3t_n(:,:,:) = 0.0_wp
1586	!
1587	neuler = 0
1588	ENDIF ! ------------ !
1589	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
1590	! ! ------------ !
1591	IF( MINVAL(id(5:7) ) > 0 ) THEN ! all required arrays exist
1592	CALL iom_get( numror, jpdom_autoglo, 'un_lf' , un_lf(:,:,:) )
1593	CALL iom_get( numror, jpdom_autoglo, 'vn_lf' , vn_lf(:,:,:) )
1594	CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:) )
1595	ELSE ! one at least array is missing
1596	un_lf(:,:,:) = 0.0_wp
1597	vn_lf(:,:,:) = 0.0_wp
1598	hdivn_lf(:,:,:) = 0.0_wp
1599	neuler = 0
1600	ENDIF
1601	ENDIF
1602	!
1603	ENDIF
1604	!
1605	ELSE !* Initialize at "rest"
1606	fse3t_b(:,:,:) = e3t_0(:,:,:)
1607	fse3t_n(:,:,:) = e3t_0(:,:,:)
1608	sshn(:,:) = 0.0_wp
1609	IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
1610	tilde_e3t_b(:,:,:) = 0.0_wp
1611	tilde_e3t_n(:,:,:) = 0.0_wp
1612	END IF
1613	IF( ln_vvl_ztilde ) THEN
1614	un_lf(:,:,:) = 0.0_wp
1615	vn_lf(:,:,:) = 0.0_wp
1616	hdivn_lf(:,:,:) = 0.0_wp
1617	ENDIF
1618	ENDIF
1619
1620	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
1621	! ! ===================
1622	IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
1623	! ! --------- !
1624	! ! all cases !
1625	! ! --------- !
1626	CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
1627	CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
1628	! ! ----------------------- !
1629	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
1630	! ! ----------------------- !
1631	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
1632	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
1633	END IF
1634	! ! -------------!
1635	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
1636	! ! ------------ !
1637	CALL iom_rstput( kt, nitrst, numrow, 'un_lf' , un_lf(:,:,:) )
1638	CALL iom_rstput( kt, nitrst, numrow, 'vn_lf' , vn_lf(:,:,:) )
1639	ENDIF
1640
1641	ENDIF
1642	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_rst')
1643
1644	END SUBROUTINE dom_vvl_rst
1645
1646	SUBROUTINE dom_vvl_ctl
1647	!!---------------------------------------------------------------------
1648	!! * ROUTINE dom_vvl_ctl *
1649	!!
1650	!! ** Purpose : Control the consistency between namelist options
1651	!! for vertical coordinate
1652	!!----------------------------------------------------------------------
1653	INTEGER :: ioptio
1654	INTEGER :: ios
1655
1656	NAMELIST/nam_vvl/ ln_vvl_zstar , ln_vvl_ztilde , &
1657	& ln_vvl_layer , ln_vvl_ztilde_as_zstar , &
1658	& ln_vvl_zstar_at_eqtor , ln_vvl_zstar_on_shelf , &
1659	& ln_vvl_adv_cn2 , ln_vvl_adv_fct , &
1660	& ln_vvl_lap , ln_vvl_blp , &
1661	& rn_ahe3_lap , rn_ahe3_blp , &
1662	& rn_rst_e3t , rn_lf_cutoff , &
1663	& ln_vvl_regrid , rn_zdef_max , &
1664	& ln_vvl_ramp , rn_day_ramp , &
1665	& ln_vvl_dbg ! not yet implemented: ln_vvl_kepe
1666	!!----------------------------------------------------------------------
1667
1668	REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist :
1669	READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
1670	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
1671
1672	REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run
1673	READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
1674	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
1675	IF(lwm) WRITE ( numond, nam_vvl )
1676
1677	IF(lwp) THEN ! Namelist print
1678	WRITE(numout,*)
1679	WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
1680	WRITE(numout,*) '~~~~~~~~~~~'
1681	WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate'
1682	WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar
1683	WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde
1684	WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer
1685	WRITE(numout,*) ' ztilde as zstar ln_vvl_ztilde_as_zstar = ', ln_vvl_ztilde_as_zstar
1686	! WRITE(numout,*) ' Namelist nam_vvl : chose kinetic-to-potential energy conservation'
1687	! WRITE(numout,*) ' ln_vvl_kepe = ', ln_vvl_kepe
1688	WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor
1689	WRITE(numout,*) ' ztilde on shelves ln_vvl_zstar_on_shelf = ', ln_vvl_zstar_on_shelf
1690	WRITE(numout,*) ' Namelist nam_vvl : thickness advection scheme'
1691	WRITE(numout,*) ' 2nd order ln_vvl_adv_cn2 = ', ln_vvl_adv_cn2
1692	WRITE(numout,*) ' 2nd order FCT ln_vvl_adv_fct = ', ln_vvl_adv_fct
1693	WRITE(numout,*) ' Namelist nam_vvl : thickness diffusion scheme'
1694	WRITE(numout,*) ' Laplacian ln_vvl_lap = ', ln_vvl_lap
1695	WRITE(numout,*) ' Bilaplacian ln_vvl_blp = ', ln_vvl_blp
1696	WRITE(numout,*) ' Laplacian coefficient rn_ahe3_lap = ', rn_ahe3_lap
1697	WRITE(numout,*) ' Bilaplacian coefficient rn_ahe3_blp = ', rn_ahe3_blp
1698	WRITE(numout,*) ' Namelist nam_vvl : layers regriding'
1699	WRITE(numout,*) ' ln_vvl_regrid = ', ln_vvl_regrid
1700	WRITE(numout,*) ' Namelist nam_vvl : linear ramp at startup'
1701	WRITE(numout,*) ' ln_vvl_ramp = ', ln_vvl_ramp
1702	WRITE(numout,*) ' rn_day_ramp = ', rn_day_ramp
1703	WRITE(numout,*) ' Namelist nam_vvl : maximum e3t deformation fractional change'
1704	WRITE(numout,*) ' rn_zdef_max = ', rn_zdef_max
1705	IF( ln_vvl_ztilde_as_zstar ) THEN
1706	WRITE(numout,*) ' ztilde running in zstar emulation mode; '
1707	WRITE(numout,*) ' ignoring namelist timescale parameters and using:'
1708	WRITE(numout,*) ' hard-wired : z-tilde to zstar restoration timescale (days)'
1709	WRITE(numout,*) ' rn_rst_e3t = 0.0'
1710	WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1711	WRITE(numout,*) ' rn_lf_cutoff = 1.0/rdt'
1712	ELSE
1713	WRITE(numout,*) ' Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
1714	WRITE(numout,*) ' rn_rst_e3t = ', rn_rst_e3t
1715	WRITE(numout,*) ' Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
1716	WRITE(numout,*) ' rn_lf_cutoff = ', rn_lf_cutoff
1717	ENDIF
1718	WRITE(numout,*) ' Namelist nam_vvl : debug prints'
1719	WRITE(numout,*) ' ln_vvl_dbg = ', ln_vvl_dbg
1720	ENDIF
1721
1722	ioptio = 0 ! Parameter control
1723	IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
1724	IF( ln_vvl_zstar ) ioptio = ioptio + 1
1725	IF( ln_vvl_ztilde ) ioptio = ioptio + 1
1726	IF( ln_vvl_layer ) ioptio = ioptio + 1
1727
1728	IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1729
1730	IF ( ln_vvl_ztilde.OR.ln_vvl_layer ) THEN
1731	ioptio = 0 ! Choose one advection scheme at most
1732	IF( ln_vvl_adv_cn2 ) ioptio = ioptio + 1
1733	IF( ln_vvl_adv_fct ) ioptio = ioptio + 1
1734	IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE thickness advection scheme in namelist nam_vvl' )
1735	ENDIF
1736
1737	IF(lwp) THEN ! Print the choice
1738	WRITE(numout,*)
1739	IF( ln_vvl_zstar ) WRITE(numout,*) ' zstar vertical coordinate is used'
1740	IF( ln_vvl_ztilde ) WRITE(numout,*) ' ztilde vertical coordinate is used'
1741	IF( ln_vvl_layer ) WRITE(numout,*) ' layer vertical coordinate is used'
1742	IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) ' to emulate a zstar coordinate'
1743	! - ML - Option not developed yet
1744	! IF( ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option used'
1745	! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option not used'
1746	ENDIF
1747
1748	! Use of "shelf horizon depths" should be allowed with s-z coordinates, but we restrict it to zco and zps
1749	! for the time being
1750	IF ( ln_sco ) THEN
1751	ll_shorizd=.FALSE.
1752	ELSE
1753	ll_shorizd=.TRUE.
1754	ENDIF
1755
1756	#if defined key_agrif
1757	IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )
1758	#endif
1759
1760	END SUBROUTINE dom_vvl_ctl
1761
1762	SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1763	!!---------------------------------------------------------------------
1764	!! * ROUTINE dom_vvl_orca_fix *
1765	!!
1766	!! ** Purpose : Correct surface weighted, horizontally interpolated,
1767	!! scale factors at locations that have been individually
1768	!! modified in domhgr. Such modifications break the
1769	!! relationship between e12t and e1u*e2u etc.
1770	!! Recompute some scale factors ignoring the modified metric.
1771	!!----------------------------------------------------------------------
1772	!! * Arguments
1773	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated
1774	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3
1775	CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors
1776	! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'
1777	!! * Local declarations
1778	INTEGER :: ji, jj, jk ! dummy loop indices
1779	INTEGER :: ij0, ij1, ii0, ii1 ! dummy loop indices
1780	INTEGER :: isrow ! index for ORCA1 starting row
1781	!! acc
1782	!! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
1783	!! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
1784	!!
1785	! ! =====================
1786	IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN ! ORCA R2 configuration
1787	! ! =====================
1788	!! acc
1789	IF( nn_cla == 0 ) THEN
1790	!
1791	ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u was modified)
1792	ij0 = 102 ; ij1 = 102
1793	DO jk = 1, jpkm1
1794	DO jj = mj0(ij0), mj1(ij1)
1795	DO ji = mi0(ii0), mi1(ii1)
1796	SELECT CASE ( pout )
1797	CASE( 'U' )
1798	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1799	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1800	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1801	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1802	CASE( 'F' )
1803	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1804	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1805	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1806	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1807	END SELECT
1808	END DO
1809	END DO
1810	END DO
1811	!
1812	ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u and e1v were modified)
1813	ij0 = 88 ; ij1 = 88
1814	DO jk = 1, jpkm1
1815	DO jj = mj0(ij0), mj1(ij1)
1816	DO ji = mi0(ii0), mi1(ii1)
1817	SELECT CASE ( pout )
1818	CASE( 'U' )
1819	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1820	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1821	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1822	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1823	CASE( 'V' )
1824	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1825	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1826	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1827	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1828	CASE( 'F' )
1829	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1830	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1831	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1832	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1833	END SELECT
1834	END DO
1835	END DO
1836	END DO
1837	ENDIF
1838
1839	ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u was modified)
1840	ij0 = 116 ; ij1 = 116
1841	DO jk = 1, jpkm1
1842	DO jj = mj0(ij0), mj1(ij1)
1843	DO ji = mi0(ii0), mi1(ii1)
1844	SELECT CASE ( pout )
1845	CASE( 'U' )
1846	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1847	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1848	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1849	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1850	CASE( 'F' )
1851	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1852	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1853	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1854	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1855	END SELECT
1856	END DO
1857	END DO
1858	END DO
1859	ENDIF
1860	!
1861	! ! =====================
1862	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
1863	! ! =====================
1864	! This dirty section will be suppressed by simplification process:
1865	! all this will come back in input files
1866	! Currently these hard-wired indices relate to configuration with
1867	! extend grid (jpjglo=332)
1868	! which had a grid-size of 362x292.
1869	isrow = 332 - jpjglo
1870	!
1871	ii0 = 282 ; ii1 = 283 ! Gibraltar Strait (e2u was modified)
1872	ij0 = 241 - isrow ; ij1 = 241 - isrow
1873	DO jk = 1, jpkm1
1874	DO jj = mj0(ij0), mj1(ij1)
1875	DO ji = mi0(ii0), mi1(ii1)
1876	SELECT CASE ( pout )
1877	CASE( 'U' )
1878	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1879	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1880	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1881	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1882	CASE( 'F' )
1883	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1884	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1885	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1886	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1887	END SELECT
1888	END DO
1889	END DO
1890	END DO
1891	!
1892	ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u was modified)
1893	ij0 = 248 - isrow ; ij1 = 248 - isrow
1894	DO jk = 1, jpkm1
1895	DO jj = mj0(ij0), mj1(ij1)
1896	DO ji = mi0(ii0), mi1(ii1)
1897	SELECT CASE ( pout )
1898	CASE( 'U' )
1899	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1900	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1901	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1902	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1903	CASE( 'F' )
1904	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1905	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1906	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1907	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1908	END SELECT
1909	END DO
1910	END DO
1911	END DO
1912	!
1913	ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v was modified)
1914	ij0 = 164 - isrow ; ij1 = 165 - isrow
1915	DO jk = 1, jpkm1
1916	DO jj = mj0(ij0), mj1(ij1)
1917	DO ji = mi0(ii0), mi1(ii1)
1918	SELECT CASE ( pout )
1919	CASE( 'V' )
1920	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1921	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1922	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1923	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1924	END SELECT
1925	END DO
1926	END DO
1927	END DO
1928	!
1929	ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
1930	ij0 = 164 - isrow ; ij1 = 165 - isrow
1931	DO jk = 1, jpkm1
1932	DO jj = mj0(ij0), mj1(ij1)
1933	DO ji = mi0(ii0), mi1(ii1)
1934	SELECT CASE ( pout )
1935	CASE( 'V' )
1936	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1937	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1938	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1939	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1940	END SELECT
1941	END DO
1942	END DO
1943	END DO
1944	!
1945	ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v was modified)
1946	ij0 = 164 - isrow ; ij1 = 165 - isrow
1947	DO jk = 1, jpkm1
1948	DO jj = mj0(ij0), mj1(ij1)
1949	DO ji = mi0(ii0), mi1(ii1)
1950	SELECT CASE ( pout )
1951	CASE( 'V' )
1952	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1953	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1954	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1955	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1956	END SELECT
1957	END DO
1958	END DO
1959	END DO
1960	!
1961	ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v was modified)
1962	ij0 = 164 - isrow ; ij1 = 165 - isrow
1963	DO jk = 1, jpkm1
1964	DO jj = mj0(ij0), mj1(ij1)
1965	DO ji = mi0(ii0), mi1(ii1)
1966	SELECT CASE ( pout )
1967	CASE( 'V' )
1968	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1969	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1970	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1971	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1972	END SELECT
1973	END DO
1974	END DO
1975	END DO
1976	!
1977	ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v was modified)
1978	ij0 = 181 - isrow ; ij1 = 182 - isrow
1979	DO jk = 1, jpkm1
1980	DO jj = mj0(ij0), mj1(ij1)
1981	DO ji = mi0(ii0), mi1(ii1)
1982	SELECT CASE ( pout )
1983	CASE( 'V' )
1984	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1985	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1986	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1987	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1988	END SELECT
1989	END DO
1990	END DO
1991	END DO
1992	!
1993	ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v was modified)
1994	ij0 = 181 - isrow ; ij1 = 182 - isrow
1995	DO jk = 1, jpkm1
1996	DO jj = mj0(ij0), mj1(ij1)
1997	DO ji = mi0(ii0), mi1(ii1)
1998	SELECT CASE ( pout )
1999	CASE( 'V' )
2000	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
2001	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
2002	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
2003	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
2004	END SELECT
2005	END DO
2006	END DO
2007	END DO
2008	ENDIF
2009	! ! =====================
2010	IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration
2011	! ! =====================
2012	!
2013	ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u was modified)
2014	ij0 = 327 ; ij1 = 327
2015	DO jk = 1, jpkm1
2016	DO jj = mj0(ij0), mj1(ij1)
2017	DO ji = mi0(ii0), mi1(ii1)
2018	SELECT CASE ( pout )
2019	CASE( 'U' )
2020	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
2021	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
2022	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
2023	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
2024	CASE( 'F' )
2025	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
2026	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
2027	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
2028	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
2029	END SELECT
2030	END DO
2031	END DO
2032	END DO
2033	!
2034	ii0 = 627 ; ii1 = 628 ! Bosphorus Strait (e2u was modified)
2035	ij0 = 343 ; ij1 = 343
2036	DO jk = 1, jpkm1
2037	DO jj = mj0(ij0), mj1(ij1)
2038	DO ji = mi0(ii0), mi1(ii1)
2039	SELECT CASE ( pout )
2040	CASE( 'U' )
2041	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
2042	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
2043	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
2044	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
2045	CASE( 'F' )
2046	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
2047	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
2048	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
2049	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
2050	END SELECT
2051	END DO
2052	END DO
2053	END DO
2054	!
2055	ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u was modified)
2056	ij0 = 232 ; ij1 = 232
2057	DO jk = 1, jpkm1
2058	DO jj = mj0(ij0), mj1(ij1)
2059	DO ji = mi0(ii0), mi1(ii1)
2060	SELECT CASE ( pout )
2061	CASE( 'U' )
2062	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
2063	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
2064	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
2065	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
2066	CASE( 'F' )
2067	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
2068	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
2069	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
2070	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
2071	END SELECT
2072	END DO
2073	END DO
2074	END DO
2075	!
2076	ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u was modified)
2077	ij0 = 232 ; ij1 = 232
2078	DO jk = 1, jpkm1
2079	DO jj = mj0(ij0), mj1(ij1)
2080	DO ji = mi0(ii0), mi1(ii1)
2081	SELECT CASE ( pout )
2082	CASE( 'U' )
2083	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
2084	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
2085	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
2086	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
2087	CASE( 'F' )
2088	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
2089	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
2090	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
2091	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
2092	END SELECT
2093	END DO
2094	END DO
2095	END DO
2096	!
2097	ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u was modified)
2098	ij0 = 270 ; ij1 = 270
2099	DO jk = 1, jpkm1
2100	DO jj = mj0(ij0), mj1(ij1)
2101	DO ji = mi0(ii0), mi1(ii1)
2102	SELECT CASE ( pout )
2103	CASE( 'U' )
2104	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
2105	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
2106	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
2107	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
2108	CASE( 'F' )
2109	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
2110	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
2111	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
2112	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
2113	END SELECT
2114	END DO
2115	END DO
2116	END DO
2117	!
2118	ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v was modified)
2119	ij0 = 232 ; ij1 = 233
2120	DO jk = 1, jpkm1
2121	DO jj = mj0(ij0), mj1(ij1)
2122	DO ji = mi0(ii0), mi1(ii1)
2123	SELECT CASE ( pout )
2124	CASE( 'V' )
2125	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
2126	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
2127	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
2128	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
2129	END SELECT
2130	END DO
2131	END DO
2132	END DO
2133	!
2134	ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v was modified)
2135	ij0 = 276 ; ij1 = 276
2136	DO jk = 1, jpkm1
2137	DO jj = mj0(ij0), mj1(ij1)
2138	DO ji = mi0(ii0), mi1(ii1)
2139	SELECT CASE ( pout )
2140	CASE( 'V' )
2141	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
2142	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
2143	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
2144	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
2145	END SELECT
2146	END DO
2147	END DO
2148	END DO
2149	ENDIF
2150	END SUBROUTINE dom_vvl_orca_fix
2151
2152	SUBROUTINE dom_vvl_zdf( kt, p2dt )
2153	!!----------------------------------------------------------------------
2154	!! * ROUTINE dom_vvl_zdf *
2155	!!
2156	!! ** Purpose : Do vertical thicknesses anomaly diffusion
2157	!!
2158	!! ** Method :
2159	!!
2160	!! ** Action :
2161	!!---------------------------------------------------------------------
2162	USE oce , ONLY: zwd => ua , zws => va ! (ua,va) used as 3D workspace
2163	!
2164	INTEGER , INTENT(in) :: kt ! ocean time-step index
2165	REAL(wp) , INTENT(in) :: p2dt ! time step
2166	!
2167	INTEGER :: ji, jj, jk ! dummy loop indices
2168	REAL(wp) :: zr_tscale
2169	REAL(wp) :: za1, za2, za3, za4, zdiff
2170	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwt
2171	!!---------------------------------------------------------------------
2172	!
2173	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_zdf')
2174	!
2175	zr_tscale = 0.5
2176	!
2177	CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt)
2178	!
2179	!
2180
2181	zwt(:,:,1:jpk) = 1.e-10
2182
2183	zwt(:,:,1) = 0.0_wp
2184	! DO jk = 2, jpkm1
2185	! DO jj = 1, jpj
2186	! DO ji = 1, jpi
2187	! zwt(ji,jj,jk) = zwt(ji,jj,jk-1) + (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
2188	! END DO
2189	! END DO
2190	! END DO
2191	!
2192	!
2193	! Set diffusivity (homogeneous to an inverse time scale)
2194	!
2195	DO jk = 2, jpkm1
2196	DO jj = 2, jpjm1
2197	DO ji = fs_2, fs_jpim1
2198	! Taper a little bit:
2199	za1 = tilde_e3t_n(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)
2200	za2 = tilde_e3t_n(ji,jj,jk )+e3t_0(ji,jj,jk )
2201	za4 = 0.5_wp * (e3t_0(ji,jj,jk-1) + e3t_0(ji,jj,jk ))
2202	za3 = 0.5_wp * (za1 + za2)
2203	zdiff = ABS(za3-za4)/za4
2204	IF (zdiff>=0.8) THEN
2205	zwt(ji,jj,jk) = zr_tscale * MIN(zdiff,1._wp) * za3 / p2dt * tmask(ji,jj,jk)
2206	! zwt(ji,jj,jk) = dsm(ji,jj)/ht_0(ji,jj)(1._wp-tanh((mbkt(ji,jj)+1-jk)0.2))*tmask(ji,jj,jk)
2207
2208	ELSE
2209	zwt(ji,jj,jk) = 0.e0*tmask(ji,jj,jk)
2210	ENDIF
2211	END DO
2212	END DO
2213	END DO
2214	!
2215	!
2216	! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
2217	DO jk = 1, jpkm1
2218	DO jj = 2, jpjm1
2219	DO ji = fs_2, fs_jpim1 ! vector opt.
2220	zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / fse3w(ji,jj,jk )
2221	zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / fse3w(ji,jj,jk+1)
2222	zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
2223	END DO
2224	END DO
2225	END DO
2226	!
2227	!! Matrix inversion from the first level
2228	!!----------------------------------------------------------------------
2229	DO jj = 2, jpjm1
2230	DO ji = fs_2, fs_jpim1
2231	zwt(ji,jj,1) = zwd(ji,jj,1)
2232	END DO
2233	END DO
2234	DO jk = 2, jpkm1
2235	DO jj = 2, jpjm1
2236	DO ji = fs_2, fs_jpim1
2237	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
2238	END DO
2239	END DO
2240	END DO
2241	!
2242	! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
2243	DO jk = 2, jpkm1
2244	DO jj = 2, jpjm1
2245	DO ji = fs_2, fs_jpim1
2246	tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * tilde_e3t_a(ji,jj,jk-1)
2247	END DO
2248	END DO
2249	END DO
2250	! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
2251	DO jj = 2, jpjm1
2252	DO ji = fs_2, fs_jpim1
2253	tilde_e3t_a(ji,jj,jpkm1) = tilde_e3t_a(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
2254	END DO
2255	END DO
2256	DO jk = jpk-2, 1, -1
2257	DO jj = 2, jpjm1
2258	DO ji = fs_2, fs_jpim1
2259	tilde_e3t_a(ji,jj,jk) = ( tilde_e3t_a(ji,jj,jk) - zws(ji,jj,jk) * tilde_e3t_a(ji,jj,jk+1) ) &
2260	& / zwt(ji,jj,jk) * tmask(ji,jj,jk)
2261	END DO
2262	END DO
2263	END DO
2264	!
2265	CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt)
2266	!
2267	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_zdf')
2268	!
2269	END SUBROUTINE dom_vvl_zdf
2270
2271	SUBROUTINE dom_vvl_zdf2( kt, p2dt )
2272	!!----------------------------------------------------------------------
2273	!! * ROUTINE dom_vvl_zdf *
2274	!!
2275	!! ** Purpose : Do vertical interface diffusion
2276	!!
2277	!! ** Method :
2278	!!
2279	!! ** Action :
2280	!!---------------------------------------------------------------------
2281	USE oce , ONLY: zwd => ua , zws => va ! (ua,va) used as 3D workspace
2282	!
2283	INTEGER , INTENT(in) :: kt ! ocean time-step index
2284	REAL(wp) , INTENT(in) :: p2dt ! time step
2285	!
2286	INTEGER :: ji, jj, jk, kbot, kbotm1 ! dummy loop indices
2287	REAL(wp) :: zr_tscale
2288	REAL(wp) :: za1, za2, za3, za4, zdiff
2289	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwt, zw
2290	!!---------------------------------------------------------------------
2291	!
2292	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_zdf')
2293	!
2294	zr_tscale = 0.5
2295	!
2296	CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt, zw)
2297	!
2298	! Compute internal interfaces depths:
2299	zw(:,:,1) = 0._wp
2300	DO jk = 2, jpk
2301	DO jj = 2, jpjm1
2302	DO ji = fs_2, fs_jpim1 ! vector opt.
2303	zw(ji,jj,jk) = zw(ji,jj,jk-1) + (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
2304	END DO
2305	END DO
2306	END DO
2307	!
2308	! Set diffusivities at interfaces
2309	zwt(:,:,:) = 0.00000001_wp * tmask(:,:,:)
2310	zwt(:,:,1) = 0._wp
2311	!
2312	!
2313	! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
2314	DO jk = 2, jpkm1
2315	DO jj = 2, jpjm1
2316	DO ji = fs_2, fs_jpim1 ! vector opt.
2317	zwi(ji,jj,jk) = - 0.5_wp * p2dt * (zwt(ji,jj,jk-1) + zwt(ji,jj,jk )) &
2318	& / fse3w(ji,jj,jk-1) / fse3t(ji,jj,jk-1) &
2319	& * ht(ji,jj)
2320	zws(ji,jj,jk) = - 0.5_wp * p2dt * (zwt(ji,jj,jk ) + zwt(ji,jj,jk+1)) &
2321	& / fse3w(ji,jj,jk ) / fse3t(ji,jj,jk ) &
2322	& * ht(ji,jj)
2323
2324	zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
2325	END DO
2326	END DO
2327	END DO
2328	! Boundary conditions (Neumann)
2329	DO jj = 2, jpjm1
2330	DO ji = fs_2, fs_jpim1
2331	zwi(ji,jj,1) = 0._wp
2332	zws(ji,jj,1) = 0._wp
2333	zwd(ji,jj,1) = 1._wp
2334	zw (ji,jj,1) = 0._wp
2335	!
2336	zwd(ji,jj,2) = zwd(ji,jj,2) + zwi(ji,jj,2)
2337	zwi(ji,jj,2) = 0._wp
2338	! zwi(ji,jj,2) = 0._wp
2339	! zws(ji,jj,2) = 0._wp
2340	! zwd(ji,jj,2) = 1._wp
2341	END DO
2342	END DO
2343	DO jj = 2, jpjm1
2344	DO ji = fs_2, fs_jpim1
2345	kbot = mbkt(ji,jj) + 1
2346	kbotm1 = mbkt(ji,jj)
2347	zwi(ji,jj,kbot ) = 0._wp
2348	zws(ji,jj,kbot ) = 0._wp
2349	zwd(ji,jj,kbot ) = 1._wp
2350	!
2351	zwd(ji,jj,kbotm1) = zwd(ji,jj,kbotm1) + zws(ji,jj,kbotm1)
2352	zws(ji,jj,kbotm1) = 0._wp
2353	! zwi(ji,jj,kbotm1) = 0._wp
2354	! zws(ji,jj,kbotm1) = 0._wp
2355	! zwd(ji,jj,kbotm1) = 1._wp
2356	END DO
2357	END DO
2358	!
2359	DO jk = 2, jpkm1
2360	DO jj = 2, jpjm1
2361	DO ji = fs_2, fs_jpim1
2362	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
2363	END DO
2364	END DO
2365	END DO
2366	!
2367	DO jk = 2, jpk
2368	DO jj = 2, jpjm1
2369	DO ji = fs_2, fs_jpim1 ! vector opt.
2370	zwi(ji,jj,jk) = zw(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) *zwi(ji,jj,jk-1)
2371	END DO
2372	END DO
2373	END DO
2374	!
2375	DO jk = jpk-1, 2, -1
2376	DO jj = 2, jpjm1
2377	DO ji = fs_2, fs_jpim1 ! vector opt.
2378	zw(ji,jj,jk) = ( zwi(ji,jj,jk) - zws(ji,jj,jk) * zw(ji,jj,jk+1) ) / zwd(ji,jj,jk)
2379	END DO
2380	END DO
2381	END DO
2382	!
2383	! Revert to thicknesses anomalies:
2384	DO jk = 1, jpkm1
2385	DO jj = 2, jpjm1
2386	DO ji = fs_2, fs_jpim1 ! vector opt.
2387	tilde_e3t_a(ji,jj,jk) = (zw(ji,jj,jk+1)-zw(ji,jj,jk)-e3t_0(ji,jj,jk))* tmask(ji,jj,jk)
2388	END DO
2389	END DO
2390	END DO
2391	!
2392	CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt, zw)
2393	!
2394	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_zdf')
2395	!
2396	END SUBROUTINE dom_vvl_zdf2
2397
2398	SUBROUTINE dom_vvl_regrid( kt )
2399	!!----------------------------------------------------------------------
2400	!! * ROUTINE dom_vvl_regrid *
2401	!!
2402	!! ** Purpose : Ensure "well-behaved" vertical grid
2403	!!
2404	!! ** Method : More or less adapted from references below.
2405	!!regrid
2406	!! ** Action : Ensure that thickness are above a given value, spaced enough
2407	!! and revert to Eulerian coordinates near the bottom.
2408	!!
2409	!! References : Bleck, R. and S. Benjamin, 1993: Regional Weather Prediction
2410	!! with a Model Combining Terrain-following and Isentropic
2411	!! coordinates. Part I: Model Description. Monthly Weather Rev.,
2412	!! 121, 1770-1785.
2413	!! Toy, M., 2011: Incorporating Condensational Heating into a
2414	!! Nonhydrostatic Atmospheric Model Based on a Hybrid Isentropic-
2415	!! Sigma Vertical Coordinate. Monthly Weather Rev., 139, 2940-2954.
2416	!!----------------------------------------------------------------------
2417	!! * Arguments
2418	INTEGER, INTENT( in ) :: kt ! time step
2419
2420	!! * Local declarations
2421	INTEGER :: ji, jj, jk ! dummy loop indices
2422	LOGICAL :: ll_chk_bot2top, ll_chk_top2bot, ll_lapdiff_cond
2423	LOGICAL :: ll_zdiff_cond, ll_blpdiff_cond
2424	INTEGER :: jkbot
2425	REAL(wp) :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd
2426	REAL(wp) :: zufim1, zufi, zvfjm1, zvfj, dzmin_int, dzmin_surf
2427	REAL(wp) :: zh_new, zh_old, zh_bef, ztmp, ztmp1, z2dt, zh_up, zh_dwn
2428	REAL(wp) :: zeu2, zev2, zfrch_stp, zfrch_rel, zfrac_bot, zscal_bot
2429	REAL(wp) :: zhdiff, zhdiff2, zvdiff, zhlim, zhlim2, zvlim
2430	REAL(wp), POINTER, DIMENSION(:,:) :: zdw
2431	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwdw, zwdw_b
2432	!!----------------------------------------------------------------------
2433
2434	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_regrid')
2435	!
2436	CALL wrk_alloc( jpi, jpj, jpk, zwdw)
2437	!
2438	! Some user defined parameters below:
2439	ll_chk_bot2top = .TRUE.
2440	ll_chk_top2bot = .TRUE.
2441	dzmin_int = 0.1_wp ! Absolute minimum depth in the interior (in meters)
2442	dzmin_surf = 1.0_wp ! Absolute minimum depth at the surface (in meters)
2443	zfrch_stp = 0.5_wp ! Maximum fractionnal thickness change in one time step (<= 1.)
2444	zfrch_rel = 0.4_wp ! Maximum relative thickness change in the vertical (<= 1.)
2445	zfrac_bot = 0.05_wp ! Fraction of bottom level allowed to change
2446	zscal_bot = 2.0_wp ! Depth lengthscale
2447	ll_zdiff_cond = .TRUE. ! Conditionnal vertical diffusion of interfaces
2448	zvdiff = 0.1_wp ! m
2449	zvlim = 0.5_wp ! max d2h/dh
2450	ll_lapdiff_cond = .TRUE. ! Conditionnal Laplacian diffusion of interfaces
2451	zhdiff = 0.2_wp ! ad.
2452	zhlim = 0.03_wp ! ad. max lap(z)*e1
2453	ll_blpdiff_cond = .TRUE. ! Conditionnal Bilaplacian diffusion of interfaces
2454	zhdiff2 = 0.2_wp ! ad.
2455	! zhlim2 = 3.e-11_wp ! max bilap(z)
2456	zhlim2 = 0.03_wp ! ad. max bilap(z)e1*3
2457	! ---------------------------------------------------------------------------------------
2458	!
2459	! Set arrays determining maximum vertical displacement at the bottom:
2460	!--------------------------------------------------------------------
2461	IF ( kt==nit000 ) THEN
2462	DO jj = 2, jpjm1
2463	DO ji = 2, jpim1
2464	jk = MIN(mbkt(ji,jj), mbkt(ji+1,jj), mbkt(ji-1,jj), mbkt(ji,jj+1), mbkt(ji,jj-1))
2465	i_int_bot(ji,jj) = jk
2466	END DO
2467	END DO
2468	dsm(:,:) = REAL( i_int_bot(:,:), wp ) ; CALL lbc_lnk(dsm(:,:),'T',1.)
2469	i_int_bot(:,:) = MAX( INT( dsm(:,:) ), 1 )
2470
2471	CALL wrk_alloc( jpi, jpj, zdw )
2472	DO jj = 2, jpjm1
2473	DO ji = 2, jpim1
2474	zdw(ji,jj) = MAX(ABS(ht_0(ji,jj)-ht_0(ji+1,jj))*umask(ji ,jj,1), &
2475	& ABS(ht_0(ji,jj)-ht_0(ji-1,jj))*umask(ji-1,jj,1), &
2476	& ABS(ht_0(ji,jj)-ht_0(ji,jj+1))*vmask(ji,jj ,1), &
2477	& ABS(ht_0(ji,jj)-ht_0(ji,jj-1))*vmask(ji,jj-1,1) )
2478	zdw(ji,jj) = MAX(zscal_bot * zdw(ji,jj), rsmall )
2479	END DO
2480	END DO
2481	CALL lbc_lnk( zdw(:,:), 'T', 1. )
2482
2483	DO jj = 2, jpjm1
2484	DO ji = 2, jpim1
2485	dsm(ji,jj) = 1._wp/16._wp * ( zdw(ji-1,jj-1) + zdw(ji+1,jj-1) &
2486	& + zdw(ji-1,jj+1) + zdw(ji+1,jj+1) &
2487	& + 2._wp*( zdw(ji ,jj-1) + zdw(ji-1,jj ) &
2488	& + zdw(ji+1,jj ) + zdw(ji ,jj+1) ) &
2489	& + 4._wp* zdw(ji ,jj ) )
2490	END DO
2491	END DO
2492
2493	CALL lbc_lnk( dsm(:,:), 'T', 1. )
2494	CALL wrk_dealloc( jpi, jpj, zdw)
2495
2496	IF (ln_zps) THEN
2497	DO jj = 1, jpj
2498	DO ji = 1, jpi
2499	jk = i_int_bot(ji,jj)
2500	hsm(ji,jj) = zfrac_bot * e3w_1d(jk)
2501	! dsm(ji,jj) = MAX(dsm(ji,jj), 0.1_wp*ht_0(ji,jj))
2502	END DO
2503	END DO
2504	ELSE
2505	DO jj = 1, jpj
2506	DO ji = 1, jpi
2507	jk = i_int_bot(ji,jj)
2508	hsm(ji,jj) = zfrac_bot * e3w_0(ji,jj,jk)
2509	! dsm(ji,jj) = MAX(dsm(ji,jj), 0.1_wp*ht_0(ji,jj))
2510	END DO
2511	END DO
2512	ENDIF
2513	END IF
2514
2515	! Provisionnal interface depths:
2516	!-------------------------------
2517	zwdw(:,:,1) = 0.e0
2518	DO jj = 1, jpj
2519	DO ji = 1, jpi
2520	DO jk = 2, jpk
2521	zwdw(ji,jj,jk) = zwdw(ji,jj,jk-1) + &
2522	& (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
2523	END DO
2524	END DO
2525	END DO
2526	!
2527	! Conditionnal horizontal Laplacian diffusion:
2528	!---------------------------------------------
2529	IF ( ll_lapdiff_cond ) THEN
2530	CALL wrk_alloc( jpi, jpj, jpk, zwdw_b)
2531	!
2532	zwdw_b(:,:,1) = 0._wp
2533	DO jj = 1, jpj
2534	DO ji = 1, jpi
2535	DO jk=2,jpk
2536	zwdw_b(ji,jj,jk) = zwdw_b(ji,jj,jk-1) + &
2537	& (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
2538	END DO
2539	END DO
2540	END DO
2541	!
2542	DO jk = 2, jpkm1
2543	DO jj = 1, jpjm1
2544	DO ji = 1, fs_jpim1 ! vector opt.
2545	ua(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) &
2546	& * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj ,jk) )
2547	va(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) &
2548	& * ( zwdw_b(ji,jj,jk) - zwdw_b(ji ,jj+1,jk) )
2549	END DO
2550	END DO
2551	END DO
2552
2553	DO jk = 2, jpkm1
2554	DO jj = 2, jpjm1
2555	DO ji = fs_2, fs_jpim1 ! vector opt.
2556	ztmp1 = ( (ua(ji-1,jj ,jk) - ua(ji,jj,jk)) &
2557	& + (va(ji ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj)
2558	zh2 = MAX(abs(ztmp1)-zhlim*SQRT(r1_e12t(ji,jj)), 0._wp)
2559	ztmp = SIGN(zh2, ztmp1)
2560	zeu2 = zhdiff * e12t(ji,jj)e12t(ji,jj)/(e1t(ji,jj)e1t(ji,jj) + e2t(ji,jj)*e2t(ji,jj))
2561	zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp * tmask(ji,jj,jk)
2562	END DO
2563	END DO
2564	END DO
2565	!
2566	CALL wrk_dealloc( jpi, jpj, jpk, zwdw_b)
2567	!
2568	ENDIF
2569
2570	! Conditionnal horizontal Bilaplacian diffusion:
2571	!-----------------------------------------------
2572	IF ( ll_blpdiff_cond ) THEN
2573	CALL wrk_alloc( jpi, jpj, jpk, zwdw_b)
2574	!
2575	zwdw_b(:,:,1) = 0._wp
2576	DO jj = 1, jpj
2577	DO ji = 1, jpi
2578	DO jk=2,jpk
2579	zwdw_b(ji,jj,jk) = zwdw_b(ji,jj,jk-1) + &
2580	& (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1)
2581	END DO
2582	END DO
2583	END DO
2584	!
2585	DO jk = 2, jpkm1
2586	DO jj = 1, jpjm1
2587	DO ji = 1, fs_jpim1 ! vector opt.
2588	ua(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) &
2589	& * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj ,jk) )
2590	va(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) &
2591	& * ( zwdw_b(ji,jj,jk) - zwdw_b(ji ,jj+1,jk) )
2592	END DO
2593	END DO
2594	END DO
2595
2596	DO jk = 2, jpkm1
2597	DO jj = 2, jpjm1
2598	DO ji = fs_2, fs_jpim1 ! vector opt.
2599	zwdw_b(ji,jj,jk) = -( (ua(ji-1,jj ,jk) - ua(ji,jj,jk)) &
2600	& + (va(ji ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj)
2601	END DO
2602	END DO
2603	END DO
2604	!
2605	CALL lbc_lnk( zwdw_b(:,:,:), 'T', 1. )
2606	!
2607	DO jk = 2, jpkm1
2608	DO jj = 1, jpjm1
2609	DO ji = 1, fs_jpim1 ! vector opt.
2610	ua(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) &
2611	& * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj ,jk) )
2612	va(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) &
2613	& * ( zwdw_b(ji,jj,jk) - zwdw_b(ji ,jj+1,jk) )
2614	END DO
2615	END DO
2616	END DO
2617	!
2618	DO jk = 2, jpkm1
2619	DO jj = 2, jpjm1
2620	DO ji = fs_2, fs_jpim1 ! vector opt.
2621	ztmp1 = ( (ua(ji-1,jj ,jk) - ua(ji,jj,jk)) &
2622	& + (va(ji ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj)
2623	zh2 = MAX(abs(ztmp1)-zhlim2SQRT(r1_e12t(ji,jj))r1_e12t(ji,jj), 0._wp)
2624	ztmp = SIGN(zh2, ztmp1)
2625	zeu2 = zhdiff2 * e12t(ji,jj)*e12t(ji,jj) / 16._wp
2626	zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp * tmask(ji,jj,jk)
2627	END DO
2628	END DO
2629	END DO
2630	!
2631	CALL wrk_dealloc( jpi, jpj, jpk, zwdw_b)
2632	!
2633	ENDIF
2634
2635	! Conditionnal vertical diffusion:
2636	!---------------------------------
2637	IF ( ll_zdiff_cond ) THEN
2638	DO jk = 2, jpkm1
2639	DO jj = 2, jpjm1
2640	DO ji = fs_2, fs_jpim1 ! vector opt.
2641	ztmp = -( (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))*tmask(ji,jj,jk-1) &
2642	-(tilde_e3t_b(ji,jj,jk )+e3t_0(ji,jj,jk ))*tmask(ji,jj,jk ) )
2643	ztmp1 = 0.5_wp * ( tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) &
2644	& +tilde_e3t_b(ji,jj,jk ) + e3t_0(ji,jj,jk ) )
2645	zh2 = MAX(abs(ztmp)-zvlim*ztmp1, 0._wp)
2646	ztmp = SIGN(zh2, ztmp)
2647	IF ((jk==mbkt(ji,jj)).AND.(ln_zps)) ztmp=0.e0
2648	zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zvdiff * ztmp * tmask(ji,jj,jk)
2649	END DO
2650	END DO
2651	END DO
2652	ENDIF
2653	!
2654	! Check grid from the bottom to the surface
2655	!------------------------------------------
2656	IF ( ll_chk_bot2top ) THEN
2657	DO jj = 2, jpjm1
2658	DO ji = 2, jpim1
2659	jkbot = mbkt(ji,jj)
2660	DO jk = jkbot,2,-1
2661	!
2662	zh_0 = e3t_0(ji,jj,jk)
2663	zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1))
2664	zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
2665	zh_min = MIN(zh_0/3._wp, dzmin_int)
2666	!
2667	! Set maximum and minimum vertical excursions
2668	ztmph = hsm(ji,jj)
2669	ztmpd = dsm(ji,jj)
2670	zh2 = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd)
2671	zdiff = cush_max(gdepw_0(ji,jj,jk)-zwdw(ji,jj,jk), zh2 )
2672	zwdw(ji,jj,jk) = MAX(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) - zdiff)
2673	zdiff = cush_max(zwdw(ji,jj,jk)-gdepw_0(ji,jj,jk), zh2 )
2674	zwdw(ji,jj,jk) = MIN(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) + zdiff)
2675	!
2676	! New layer thickness:
2677	zh_new = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
2678	!
2679	! Ensure minimum layer thickness:
2680	! zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
2681	zh_new = cush(zh_new, zh_min)
2682	!
2683	! Final flux:
2684	zdiff = (zh_new - zh_old) * tmask(ji,jj,jk)
2685	!
2686	! Limit thickness change in 1 time step:
2687	ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
2688	zdiff = SIGN(ztmp, zh_new - zh_old)
2689	zh_new = zdiff + zh_old
2690	!
2691	zwdw(ji,jj,jk) = zwdw(ji,jj,jk+1) - zh_new
2692	END DO
2693	END DO
2694	END DO
2695	END IF
2696	!
2697	! Check grid from the surface to the bottom
2698	!------------------------------------------
2699	IF ( ll_chk_top2bot ) THEN
2700	DO jj = 2, jpjm1
2701	DO ji = 2, jpim1
2702	jkbot = mbkt(ji,jj)
2703	DO jk = 1, jkbot-1
2704	!
2705	zh_0 = e3t_0(ji,jj,jk)
2706	zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk+1) + e3t_0(ji,jj,jk+1))
2707	zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
2708	zh_min = MIN(zh_0/3._wp, dzmin_int)
2709	!
2710	zwdw(ji,jj,jk+1) = MAX(zwdw(ji,jj,jk+1), REAL(jk)*dzmin_surf)
2711	!
2712	! New layer thickness:
2713	zh_new = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk)
2714	!
2715	! Ensure minimum layer thickness:
2716	! zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new)
2717	zh_new = cush(zh_new, zh_min)
2718	!
2719	! Final flux:
2720	zdiff = (zh_new -zh_old) * tmask(ji,jj,jk)
2721	!
2722	! Limit flux:
2723	ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef )
2724	zdiff = SIGN(ztmp, zh_new - zh_old)
2725	zh_new = zdiff + zh_old
2726	!
2727	zwdw(ji,jj,jk+1) = zwdw(ji,jj,jk) + zh_new
2728	END DO
2729	!
2730	END DO
2731	END DO
2732	ENDIF
2733	!
2734	DO jj = 2, jpjm1
2735	DO ji = 2, jpim1
2736	DO jk = 1, jpkm1
2737	tilde_e3t_a(ji,jj,jk) = (zwdw(ji,jj,jk+1)-zwdw(ji,jj,jk)-e3t_0(ji,jj,jk)) * tmask(ji,jj,jk)
2738	END DO
2739	END DO
2740	END DO
2741	!
2742	CALL wrk_dealloc( jpi, jpj, jpk, zwdw )
2743	!
2744	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_regrid')
2745	!
2746	END SUBROUTINE dom_vvl_regrid
2747
2748	FUNCTION cush(hin, hmin) RESULT(hout)
2749	!!----------------------------------------------------------------------
2750	!! * FUNCTION cush *
2751	!!
2752	!! ** Purpose :
2753	!!
2754	!! ** Method :
2755	!!
2756	!!----------------------------------------------------------------------
2757	IMPLICIT NONE
2758	REAL(wp), INTENT(in) :: hin, hmin
2759	REAL(wp) :: hout, zx, zh_cri
2760	!!----------------------------------------------------------------------
2761	zh_cri = 3._wp * hmin
2762	!
2763	IF ( hin<=0._wp ) THEN
2764	hout = hmin
2765	!
2766	ELSEIF ( (hin>0._wp).AND.(hin<=zh_cri) ) THEN
2767	zx = hin/zh_cri
2768	hout = hmin * (1._wp + zx + zx*zx)
2769	!
2770	ELSEIF ( hin>zh_cri ) THEN
2771	hout = hin
2772	!
2773	ENDIF
2774	!
2775	END FUNCTION cush
2776
2777	FUNCTION cush_max(hin, hmax) RESULT(hout)
2778	!!----------------------------------------------------------------------
2779	!! * FUNCTION cush *
2780	!!
2781	!! ** Purpose :
2782	!!
2783	!! ** Method :
2784	!!
2785	!!----------------------------------------------------------------------
2786	IMPLICIT NONE
2787	REAL(wp), INTENT(in) :: hin, hmax
2788	REAL(wp) :: hout, hmin, zx, zh_cri
2789	!!----------------------------------------------------------------------
2790	hmin = 0.1_wp * hmax
2791	zh_cri = 3._wp * hmin
2792	!
2793	IF ( (hin>=(hmax-zh_cri)).AND.(hin<=(hmax-hmin))) THEN
2794	zx = (hmax-hin)/zh_cri
2795	hout = hmax - hmin * (1._wp + zx + zx*zx)
2796	!
2797	ELSEIF ( hin>(hmax-zh_cri) ) THEN
2798	hout = hmax - hmin
2799	!
2800	ELSE
2801	hout = hin
2802	!
2803	ENDIF
2804	!
2805	END FUNCTION cush_max
2806
2807	SUBROUTINE dom_vvl_adv_fct( kt, pta, uin, vin )
2808	!!----------------------------------------------------------------------
2809	!! * ROUTINE dom_vvl_adv_fct *
2810	!!
2811	!! ** Purpose : Do thickness advection
2812	!!
2813	!! ** Method : FCT scheme to ensure positivity
2814	!!
2815	!! ** Action : - Update pta thickness tendency and diffusive fluxes
2816	!! - this is the total trend, hence it does include sea level motions
2817	!!----------------------------------------------------------------------
2818	!
2819	INTEGER, INTENT( in ) :: kt ! ocean time-step index
2820	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pta ! thickness baroclinic trend
2821	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: uin, vin ! input velocities
2822	!
2823	INTEGER :: ji, jj, jk, ib, ib_bdy ! dummy loop indices
2824	INTEGER :: ikbu, ikbv, ibot
2825	REAL(wp) :: z2dtt, zbtr, ztra ! local scalar
2826	REAL(wp) :: zdi, zdj, zmin ! - -
2827	REAL(wp) :: zfp_ui, zfp_vj ! - -
2828	REAL(wp) :: zfm_ui, zfm_vj ! - -
2829	REAL(wp) :: zfp_hi, zfp_hj ! - -
2830	REAL(wp) :: zfm_hi, zfm_hj ! - -
2831	REAL(wp) :: ztout, ztin, zfac ! - -
2832	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy, zwi
2833	!!----------------------------------------------------------------------
2834	!
2835	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_adv_fct')
2836	!
2837	CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy, zwi)
2838	!
2839	! 1. Initializations
2840	! ------------------
2841	!
2842	IF( neuler == 0 .AND. kt == nit000 ) THEN
2843	z2dtt = rdt
2844	ELSE
2845	z2dtt = 2.0_wp * rdt
2846	ENDIF
2847	!
2848	zwi(:,:,:) = 0.e0
2849	zwx(:,:,:) = 0.e0
2850	zwy(:,:,:) = 0.e0
2851	!
2852	!
2853	! 2. upstream advection with initial mass fluxes & intermediate update
2854	! --------------------------------------------------------------------
2855	IF ( ll_shorizd ) THEN
2856	DO jk = 1, jpkm1
2857	DO jj = 1, jpjm1
2858	DO ji = 1, fs_jpim1 ! vector opt.
2859	!
2860	zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji ,jj ,jk), 0._wp)
2861	zfp_hi = MIN(zfp_hi, fse3t_b(ji ,jj ,jk))
2862	zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) )
2863	!
2864	zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj ,jk), 0._wp)
2865	zfm_hi = MIN(zfm_hi, fse3t_b(ji+1,jj ,jk))
2866	zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) )
2867	!
2868	zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji ,jj ,jk), 0._wp)
2869	zfp_hj = MIN(zfp_hj, fse3t_b(ji ,jj ,jk))
2870	zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) )
2871	!
2872	zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji ,jj+1,jk), 0._wp)
2873	zfm_hj = MIN(zfm_hj, fse3t_b(ji ,jj+1,jk))
2874	zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) )
2875	!
2876	zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) )
2877	zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) )
2878	zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) )
2879	zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) )
2880	zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
2881	zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
2882	END DO
2883	END DO
2884	END DO
2885	ELSE
2886	DO jk = 1, jpkm1
2887	DO jj = 1, jpjm1
2888	DO ji = 1, fs_jpim1 ! vector opt.
2889	!
2890	zfp_hi = fse3t_b(ji ,jj ,jk)
2891	zfm_hi = fse3t_b(ji+1,jj ,jk)
2892	zfp_hj = fse3t_b(ji ,jj ,jk)
2893	zfm_hj = fse3t_b(ji ,jj+1,jk)
2894	!
2895	zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) )
2896	zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) )
2897	zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) )
2898	zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) )
2899	zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
2900	zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
2901	END DO
2902	END DO
2903	END DO
2904	ENDIF
2905
2906	! total advective trend
2907	DO jk = 1, jpkm1
2908	DO jj = 2, jpjm1
2909	DO ji = fs_2, fs_jpim1 ! vector opt.
2910	zbtr = r1_e12t(ji,jj)
2911	! total intermediate advective trends
2912	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
2913	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) )
2914	!
2915	! update and guess with monotonic sheme
2916	pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
2917	zwi(ji,jj,jk) = (fse3t_b(ji,jj,jk) + z2dtt * ztra ) * tmask(ji,jj,jk)
2918	END DO
2919	END DO
2920	END DO
2921
2922	CALL lbc_lnk( zwi, 'T', 1. )
2923
2924	#if defined key_bdy
2925	DO ib_bdy=1, nb_bdy
2926	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1)
2927	ji = idx_bdy(ib_bdy)%nbi(ib,1)
2928	jj = idx_bdy(ib_bdy)%nbj(ib,1)
2929	DO jk = 1, jpkm1
2930	zwi(ji,jj,jk) = fse3t_a(ji,jj,jk)
2931	END DO
2932	END DO
2933	END DO
2934	#endif
2935
2936	IF ( ln_vvl_dbg ) THEN
2937	zmin = MINVAL( zwi(:,:,:), mask = tmask(:,:,:) == 1.e0 )
2938	IF( lk_mpp ) CALL mpp_min( zmin )
2939	IF( zmin < 0._wp) THEN
2940	IF(lwp) CALL ctl_warn('vvl_adv: CFL issue here')
2941	IF(lwp) WRITE(numout,*) zmin
2942	ENDIF
2943	ENDIF
2944
2945	! 3. antidiffusive flux : high order minus low order
2946	! --------------------------------------------------
2947	! antidiffusive flux on i and j
2948	DO jk = 1, jpkm1
2949	DO jj = 1, jpjm1
2950	DO ji = 1, fs_jpim1 ! vector opt.
2951	zwx(ji,jj,jk) = (e2u(ji,jj) * uin(ji,jj,jk) * fse3u_n(ji,jj,jk) &
2952	& - zwx(ji,jj,jk)) * umask(ji,jj,jk)
2953	zwy(ji,jj,jk) = (e1v(ji,jj) * vin(ji,jj,jk) * fse3v_n(ji,jj,jk) &
2954	& - zwy(ji,jj,jk)) * vmask(ji,jj,jk)
2955	!
2956	! Update advective fluxes
2957	un_td(ji,jj,jk) = un_td(ji,jj,jk) - zwx(ji,jj,jk)
2958	vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - zwy(ji,jj,jk)
2959	END DO
2960	END DO
2961	END DO
2962
2963	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral boundary conditions
2964
2965	! 4. monotonicity algorithm
2966	! -------------------------
2967	CALL nonosc_2d( fse3t_b(:,:,:), zwx, zwy, zwi, z2dtt )
2968
2969	! 5. final trend with corrected fluxes
2970	! ------------------------------------
2971	!
2972	! Update advective fluxes
2973	un_td(:,:,:) = (un_td(:,:,:) + zwx(:,:,:))*umask(:,:,:)
2974	vn_td(:,:,:) = (vn_td(:,:,:) + zwy(:,:,:))*vmask(:,:,:)
2975	!
2976	DO jk = 1, jpkm1
2977	DO jj = 2, jpjm1
2978	DO ji = fs_2, fs_jpim1 ! vector opt.
2979	!
2980	zbtr = r1_e12t(ji,jj)
2981	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
2982	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) )
2983	! add them to the general tracer trends
2984	pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
2985	END DO
2986	END DO
2987	END DO
2988	!
2989	CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy, zwi)
2990	!
2991	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_adv_fct')
2992	!
2993	END SUBROUTINE dom_vvl_adv_fct
2994
2995	SUBROUTINE dom_vvl_ups_cor( kt, pta, uin, vin )
2996	!!----------------------------------------------------------------------
2997	!! * ROUTINE dom_vvl_adv_fct *
2998	!!
2999	!! ** Purpose : Correct for addionnal barotropic fluxes
3000	!! in the upstream direction
3001	!!
3002	!! ** Method :
3003	!!
3004	!! ** Action : - Update diffusive fluxes uin, vin
3005	!! - Remove divergence from thickness tendency
3006	!!----------------------------------------------------------------------
3007	INTEGER, INTENT( in ) :: kt ! ocean time-step index
3008	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pta ! thickness baroclinic trend
3009	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: uin, vin ! input fluxes
3010	INTEGER :: ji, jj, jk ! dummy loop indices
3011	INTEGER :: ikbu, ikbv, ibot
3012	REAL(wp) :: zbtr, ztra ! local scalar
3013	REAL(wp) :: zdi, zdj ! - -
3014	REAL(wp) :: zfp_hi, zfp_hj ! - -
3015	REAL(wp) :: zfm_hi, zfm_hj ! - -
3016	REAL(wp) :: zfp_ui, zfp_vj ! - -
3017	REAL(wp) :: zfm_ui, zfm_vj ! - -
3018	REAL(wp), POINTER, DIMENSION(:,:) :: zbu, zbv, zhu_b, zhv_b
3019	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy
3020	!!----------------------------------------------------------------------
3021	!
3022	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_ups_cor')
3023	!
3024	CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv)
3025	CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy)
3026	!
3027	! Compute barotropic flux difference:
3028	zbu(:,:) = 0.e0
3029	zbv(:,:) = 0.e0
3030	DO jj = 1, jpj
3031	DO ji = 1, jpi ! vector opt.
3032	DO jk = 1, jpkm1
3033	zbu(ji,jj) = zbu(ji,jj) - uin(ji,jj,jk) * umask(ji,jj,jk)
3034	zbv(ji,jj) = zbv(ji,jj) - vin(ji,jj,jk) * vmask(ji,jj,jk)
3035	END DO
3036	END DO
3037	ENDDO
3038
3039	! Compute upstream depths:
3040	zhu_b(:,:) = 0.e0
3041	zhv_b(:,:) = 0.e0
3042
3043	IF ( ll_shorizd ) THEN
3044	! Correct bottom value
3045	! considering "shelf horizon depth"
3046	DO jj = 1, jpjm1
3047	DO ji = 1, jpim1 ! vector opt.
3048	zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))
3049	zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
3050	DO jk=1, jpkm1
3051	zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji ,jj ,jk), 0._wp)
3052	zfp_hi = MIN(zfp_hi, fse3t_b(ji ,jj ,jk))
3053	zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) )
3054	!
3055	zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj ,jk), 0._wp)
3056	zfm_hi = MIN(zfm_hi, fse3t_b(ji+1,jj ,jk))
3057	zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) )
3058	!
3059	zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji ,jj ,jk), 0._wp)
3060	zfp_hj = MIN(zfp_hj, fse3t_b(ji ,jj ,jk))
3061	zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) )
3062	!
3063	zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji ,jj+1,jk), 0._wp)
3064	zfm_hj = MIN(zfm_hj, fse3t_b(ji ,jj+1,jk))
3065	zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) )
3066	!
3067	zhu_b(ji,jj) = zhu_b(ji,jj) + ( zdi * zfp_hi &
3068	& + (1._wp-zdi) * zfm_hi &
3069	& ) * umask(ji,jj,jk)
3070	zhv_b(ji,jj) = zhv_b(ji,jj) + ( zdj * zfp_hj &
3071	& + (1._wp-zdj) * zfm_hj &
3072	& ) * vmask(ji,jj,jk)
3073	END DO
3074	END DO
3075	END DO
3076	ELSE
3077	DO jj = 1, jpjm1
3078	DO ji = 1, jpim1 ! vector opt.
3079	zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj))
3080	zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj))
3081	DO jk = 1, jpkm1
3082	zfp_hi = fse3t_b(ji ,jj ,jk)
3083	zfm_hi = fse3t_b(ji+1,jj ,jk)
3084	zfp_hj = fse3t_b(ji ,jj ,jk)
3085	zfm_hj = fse3t_b(ji ,jj+1,jk)
3086	!
3087	zhu_b(ji,jj) = zhu_b(ji,jj) + ( zdi * zfp_hi &
3088	& + (1._wp-zdi) * zfm_hi &
3089	& ) * umask(ji,jj,jk)
3090	!
3091	zhv_b(ji,jj) = zhv_b(ji,jj) + ( zdj * zfp_hj &
3092	& + (1._wp-zdj) * zfm_hj &
3093	& ) * vmask(ji,jj,jk)
3094	END DO
3095	END DO
3096	END DO
3097	ENDIF
3098
3099	! Corrective barotropic velocity (times hor. scale factor)
3100	zbu(:,:) = zbu(:,:)/ (zhu_b(:,:)*umask(:,:,1)+1._wp-umask(:,:,1))
3101	zbv(:,:) = zbv(:,:)/ (zhv_b(:,:)*vmask(:,:,1)+1._wp-vmask(:,:,1))
3102
3103	CALL lbc_lnk( zbu(:,:), 'U', -1. )
3104	CALL lbc_lnk( zbv(:,:), 'V', -1. )
3105
3106	! Set corrective fluxes in upstream direction:
3107	!
3108	zwx(:,:,:) = 0.e0
3109	zwy(:,:,:) = 0.e0
3110	IF ( ll_shorizd ) THEN
3111	DO jj = 1, jpjm1
3112	DO ji = 1, fs_jpim1 ! vector opt.
3113	! upstream scheme
3114	zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) )
3115	zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) )
3116	zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) )
3117	zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) )
3118	DO jk = 1, jpkm1
3119	zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji ,jj ,jk), 0._wp)
3120	zfp_hi = MIN(fse3t_b(ji ,jj ,jk), zfp_hi)
3121	zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) )
3122	!
3123	zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj ,jk), 0._wp)
3124	zfm_hi = MIN(fse3t_b(ji+1,jj ,jk), zfm_hi)
3125	zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) )
3126	!
3127	zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji ,jj ,jk), 0._wp)
3128	zfp_hj = MIN(fse3t_b(ji ,jj ,jk), zfp_hj)
3129	zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) )
3130	!
3131	zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji ,jj+1,jk), 0._wp)
3132	zfm_hj = MIN(fse3t_b(ji ,jj+1,jk), zfm_hj)
3133	zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) )
3134	!
3135	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
3136	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
3137	END DO
3138	END DO
3139	END DO
3140	ELSE
3141	DO jj = 1, jpjm1
3142	DO ji = 1, fs_jpim1 ! vector opt.
3143	! upstream scheme
3144	zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) )
3145	zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) )
3146	zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) )
3147	zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) )
3148	DO jk = 1, jpkm1
3149	zfp_hi = fse3t_b(ji ,jj ,jk)
3150	zfm_hi = fse3t_b(ji+1,jj ,jk)
3151	zfp_hj = fse3t_b(ji ,jj ,jk)
3152	zfm_hj = fse3t_b(ji ,jj+1,jk)
3153	!
3154	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk)
3155	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk)
3156	END DO
3157	END DO
3158	END DO
3159	ENDIF
3160	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral boundary conditions
3161
3162	uin(:,:,:) = uin(:,:,:) + zwx(:,:,:)
3163	vin(:,:,:) = vin(:,:,:) + zwy(:,:,:)
3164	!
3165	! Update trend with corrective fluxes:
3166	DO jk = 1, jpkm1
3167	DO jj = 2, jpjm1
3168	DO ji = fs_2, fs_jpim1 ! vector opt.
3169	!
3170	zbtr = r1_e12t(ji,jj)
3171	! total advective trends
3172	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
3173	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) )
3174	! add them to the general tracer trends
3175	pta(ji,jj,jk) = pta(ji,jj,jk) + ztra
3176	END DO
3177	END DO
3178	END DO
3179	!
3180	CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv)
3181	CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy)
3182	!
3183	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_ups_cor')
3184	!
3185	END SUBROUTINE dom_vvl_ups_cor
3186
3187	SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, p2dt )
3188	!!---------------------------------------------------------------------
3189	!! * ROUTINE nonosc_2d *
3190	!!
3191	!! ** Purpose : compute monotonic thickness fluxes from the upstream
3192	!! scheme and the before field by a nonoscillatory algorithm
3193	!!
3194	!! ** Method : ... ???
3195	!! warning : pbef and paft must be masked, but the boundaries
3196	!! conditions on the fluxes are not necessary zalezak (1979)
3197	!! drange (1995) multi-dimensional forward-in-time and upstream-
3198	!! in-space based differencing for fluid
3199	!!----------------------------------------------------------------------
3200	!
3201	!!----------------------------------------------------------------------
3202	REAL(wp) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step
3203	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field
3204	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb ! monotonic fluxes in the 3 directions
3205	!
3206	INTEGER :: ji, jj, jk ! dummy loop indices
3207	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars
3208	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
3209	REAL(wp) :: zupip1, zupim1, zupjp1, zupjm1, zupb, zupa
3210	REAL(wp) :: zdoip1, zdoim1, zdojp1, zdojm1, zdob, zdoa
3211	REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo
3212	!!----------------------------------------------------------------------
3213	!
3214	IF( nn_timing == 1 ) CALL timing_start('nonosc2')
3215	!
3216	CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
3217	!
3218
3219	zbig = 1.e+40_wp
3220	zrtrn = 1.e-15_wp
3221	zbetup(:,:,jpk) = 0._wp ; zbetdo(:,:,jpk) = 0._wp
3222
3223
3224	! Search local extrema
3225	! --------------------
3226	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
3227	zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ), &
3228	& paft * tmask - zbig * ( 1.e0 - tmask ) )
3229	zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ), &
3230	& paft * tmask + zbig * ( 1.e0 - tmask ) )
3231
3232	DO jk = 1, jpkm1
3233	z2dtt = p2dt
3234	DO jj = 2, jpjm1
3235	DO ji = fs_2, fs_jpim1 ! vector opt.
3236
3237	! search maximum in neighbourhood
3238	zup = MAX( zbup(ji ,jj ,jk ), &
3239	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
3240	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ))
3241
3242	! search minimum in neighbourhood
3243	zdo = MIN( zbdo(ji ,jj ,jk ), &
3244	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
3245	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ))
3246
3247	! positive part of the flux
3248	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
3249	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) )
3250
3251	! negative part of the flux
3252	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
3253	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) )
3254
3255	! up & down beta terms
3256	zbt = e1t(ji,jj) * e2t(ji,jj) / z2dtt
3257	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
3258	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
3259	END DO
3260	END DO
3261	END DO
3262
3263	CALL lbc_lnk( zbetup, 'T', 1. ) ; CALL lbc_lnk( zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
3264
3265	! 3. monotonic flux in the i & j direction (paa & pbb)
3266	! ----------------------------------------
3267	DO jk = 1, jpkm1
3268	DO jj = 2, jpjm1
3269	DO ji = fs_2, fs_jpim1 ! vector opt.
3270	zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
3271	zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
3272	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
3273	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
3274
3275	zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
3276	zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
3277	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
3278	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( <