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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

domvvl.F90 in NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM/domvvl.F90 @ 10164

Last change on this file since 10164 was 10164, checked in by jchanut, 6 years ago
ztilde update, #2126: Add higher order filtering + alternative scale factor decomposition
Property svn:keywords set to `Id`
File size: 120.8 KB

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

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

Download in other formats: