Context Navigation

domvvl.F90 @ 9551

Last change on this file since 9551 was 9551, checked in by jchanut, 6 years ago
Remove tracer horizontal diffusion in test case + correct interfaces interpolation after understanding it is related to ORCA2 edited scale factors
Property svn:keywords set to `Id`
File size: 156.6 KB

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

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

Download in other formats: