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

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

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

Download in other formats: