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

Last change on this file since 9529 was 9529, checked in by jchanut, 7 years ago
ztilde stability update
Property svn:keywords set to `Id`
File size: 155.7 KB

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

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

Download in other formats: