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

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

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

Download in other formats: