Context Navigation

domvvl.F90 @ 15610

Last change on this file since 15610 was 15610, checked in by jchanut, 3 years ago
#1791: add namelist parameter nn_vvl_interp to control how scale factors are set at U-V-F points. nn_vvl_interp=0 is the old interpolation scheme that do not account for steps. nn_vvl_interp=1, corrects the bottom thicknesses, but does not ensure that they get too small for stability. nn_vvl_interp=2 is a "qco like" formulation, for which scale factors anomalies are set proportional to the scale factors at rest. Set nn_vvl_interp=2 as the default.
Property svn:keywords set to `Id`
File size: 64.1 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	!!----------------------------------------------------------------------
11
12	!!----------------------------------------------------------------------
13	!! dom_vvl_init : define initial vertical scale factors, depths and column thickness
14	!! dom_vvl_zgr : most part of dom_vvl_init
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
35	IMPLICIT NONE
36	PRIVATE
37
38	PUBLIC dom_vvl_init ! called by domain.F90
39	PUBLIC dom_vvl_zgr ! called by iceistate.F90
40	PUBLIC dom_vvl_sf_nxt ! called by step.F90
41	PUBLIC dom_vvl_sf_swp ! called by step.F90
42	PUBLIC dom_vvl_interpol ! called by dynnxt.F90
43
44	! !!* Namelist nam_vvl
45	LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate
46	LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate
47	LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate
48	LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate
49	LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate
50	LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer
51	!
52	INTEGER :: nn_vvl_interp ! scale factors anomaly interpolation method at U-V-F points
53	! =0 linear with no bottom correction over steps (old)
54	! =1 linear with bottom correction over steps
55	! =2 "qco like", i.e. proportional to thicknesses at rest
56	!
57	! ! conservation: not used yet
58	REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient
59	REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days]
60	REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days]
61	REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation
62	LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints
63
64	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport
65	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence
66	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors
67	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors
68	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors
69	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence
70
71	!! * Substitutions
72	# include "vectopt_loop_substitute.h90"
73	!!----------------------------------------------------------------------
74	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
75	!! $Id$
76	!! Software governed by the CeCILL license (see ./LICENSE)
77	!!----------------------------------------------------------------------
78	CONTAINS
79
80	INTEGER FUNCTION dom_vvl_alloc()
81	!!----------------------------------------------------------------------
82	!! * FUNCTION dom_vvl_alloc *
83	!!----------------------------------------------------------------------
84	IF( ln_vvl_zstar ) dom_vvl_alloc = 0
85	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
86	ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , &
87	& dtilde_e3t_a(jpi,jpj,jpk) , un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , &
88	& STAT = dom_vvl_alloc )
89	CALL mpp_sum ( 'domvvl', dom_vvl_alloc )
90	IF( dom_vvl_alloc /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_alloc: failed to allocate arrays' )
91	un_td = 0._wp
92	vn_td = 0._wp
93	ENDIF
94	IF( ln_vvl_ztilde ) THEN
95	ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
96	CALL mpp_sum ( 'domvvl', dom_vvl_alloc )
97	IF( dom_vvl_alloc /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_alloc: failed to allocate arrays' )
98	ENDIF
99	!
100	END FUNCTION dom_vvl_alloc
101
102
103	SUBROUTINE dom_vvl_init
104	!!----------------------------------------------------------------------
105	!! * ROUTINE dom_vvl_init *
106	!!
107	!! ** Purpose : Initialization of all scale factors, depths
108	!! and water column heights
109	!!
110	!! ** Method : - use restart file and/or initialize
111	!! - interpolate scale factors
112	!!
113	!! ** Action : - e3t_(n/b) and tilde_e3t_(n/b)
114	!! - Regrid: e3(u/v)_n
115	!! e3(u/v)_b
116	!! e3w_n
117	!! e3(u/v)w_b
118	!! e3(u/v)w_n
119	!! gdept_n, gdepw_n and gde3w_n
120	!! - h(t/u/v)_0
121	!! - frq_rst_e3t and frq_rst_hdv
122	!!
123	!! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling.
124	!!----------------------------------------------------------------------
125	!
126	IF(lwp) WRITE(numout,*)
127	IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
128	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
129	!
130	CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer)
131	!
132	! ! Allocate module arrays
133	IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
134	!
135	! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
136	CALL dom_vvl_rst( nit000, 'READ' )
137	e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all
138	!
139	! !== Set of all other vertical scale factors ==! (now and before)
140	CALL dom_vvl_zgr
141	!
142	IF(lwxios) THEN
143	! define variables in restart file when writing with XIOS
144	CALL iom_set_rstw_var_active('e3t_b')
145	CALL iom_set_rstw_var_active('e3t_n')
146	! ! ----------------------- !
147	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
148	! ! ----------------------- !
149	CALL iom_set_rstw_var_active('tilde_e3t_b')
150	CALL iom_set_rstw_var_active('tilde_e3t_n')
151	END IF
152	! ! -------------!
153	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
154	! ! ------------ !
155	CALL iom_set_rstw_var_active('hdiv_lf')
156	ENDIF
157	!
158	ENDIF
159	!
160	END SUBROUTINE dom_vvl_init
161
162	SUBROUTINE dom_vvl_zgr
163	!!----------------------------------------------------------------------
164	!! * ROUTINE dom_vvl_zgr *
165	!!
166	!! ** Purpose : Interpolation of all scale factors,
167	!! depths and water column heights
168	!!
169	!! ** Method : - interpolate scale factors
170	!!
171	!! ** Action : - e3t_(n/b) and tilde_e3t_(n/b)
172	!! - Regrid: e3(u/v)_n
173	!! e3(u/v)_b
174	!! e3w_n
175	!! e3(u/v)w_b
176	!! e3(u/v)w_n
177	!! gdept_n, gdepw_n and gde3w_n
178	!! - h(t/u/v)_0
179	!! - frq_rst_e3t and frq_rst_hdv
180	!!
181	!! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling.
182	!!----------------------------------------------------------------------
183	INTEGER :: ji, jj, jk
184	INTEGER :: ii0, ii1, ij0, ij1
185	REAL(wp):: zcoef
186	!!----------------------------------------------------------------------
187	! !== Set of all other vertical scale factors ==! (now and before)
188	! ! Horizontal interpolation of e3t
189	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U
190	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
191	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V
192	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
193	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F
194	! ! Vertical interpolation of e3t,u,v
195	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W
196	CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W' )
197	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW
198	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
199	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW
200	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
201
202	! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
203	e3t_a(:,:,:) = e3t_n(:,:,:)
204	e3u_a(:,:,:) = e3u_n(:,:,:)
205	e3v_a(:,:,:) = e3v_n(:,:,:)
206	!
207	! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep)
208	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration)
209	gdepw_n(:,:,1) = 0.0_wp
210	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg
211	gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
212	gdepw_b(:,:,1) = 0.0_wp
213	DO jk = 2, jpk ! vertical sum
214	DO jj = 1, jpj
215	DO ji = 1, jpi
216	! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
217	! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
218	! ! 0.5 where jk = mikt
219	!!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ??
220	zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
221	gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
222	gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) &
223	& + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk))
224	gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
225	gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
226	gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) &
227	& + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk))
228	END DO
229	END DO
230	END DO
231	!
232	! !== thickness of the water column !! (ocean portion only)
233	ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) ....
234	hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
235	hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
236	hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
237	hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
238	DO jk = 2, jpkm1
239	ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
240	hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
241	hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
242	hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
243	hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
244	END DO
245	!
246	! !== inverse of water column thickness ==! (u- and v- points)
247	r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF
248	r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
249	r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
250	r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
251
252	! !== z_tilde coordinate case ==! (Restoring frequencies)
253	IF( ln_vvl_ztilde ) THEN
254	!!gm : idea: add here a READ in a file of custumized restoring frequency
255	! ! Values in days provided via the namelist
256	! ! use rsmall to avoid possible division by zero errors with faulty settings
257	frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp )
258	frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
259	!
260	IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile
261	frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings
262	frq_rst_hdv(:,:) = 1._wp / rdt
263	ENDIF
264	IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator
265	DO jj = 1, jpj
266	DO ji = 1, jpi
267	!!gm case \|gphi\| >= 6 degrees is useless initialized just above by default
268	IF( ABS(gphit(ji,jj)) >= 6.) THEN
269	! values outside the equatorial band and transition zone (ztilde)
270	frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp )
271	frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
272	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star
273	! values inside the equatorial band (ztilde as zstar)
274	frq_rst_e3t(ji,jj) = 0.0_wp
275	frq_rst_hdv(ji,jj) = 1.0_wp / rdt
276	ELSE ! transition band (2.5 to 6 degrees N/S)
277	! ! (linearly transition from z-tilde to z-star)
278	frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp &
279	& * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
280	& * 180._wp / 3.5_wp ) )
281	frq_rst_hdv(ji,jj) = (1.0_wp / rdt) &
282	& + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp &
283	& * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
284	& * 180._wp / 3.5_wp ) )
285	ENDIF
286	END DO
287	END DO
288	IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN
289	IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
290	ii0 = 103 ; ii1 = 111
291	ij0 = 128 ; ij1 = 135 ;
292	frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp
293	frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rdt
294	ENDIF
295	ENDIF
296	ENDIF
297	ENDIF
298	!
299	END SUBROUTINE dom_vvl_zgr
300
301	SUBROUTINE dom_vvl_sf_nxt( kt, kcall )
302	!!----------------------------------------------------------------------
303	!! * ROUTINE dom_vvl_sf_nxt *
304	!!
305	!! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt,
306	!! tranxt and dynspg routines
307	!!
308	!! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness.
309	!! - z_tilde_case: after scale factor increment =
310	!! high frequency part of horizontal divergence
311	!! + retsoring towards the background grid
312	!! + thickness difusion
313	!! Then repartition of ssh INCREMENT proportionnaly
314	!! to the "baroclinic" level thickness.
315	!!
316	!! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case
317	!! - tilde_e3t_a: after increment of vertical scale factor
318	!! in z_tilde case
319	!! - e3(t/u/v)_a
320	!!
321	!! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling.
322	!!----------------------------------------------------------------------
323	INTEGER, INTENT( in ) :: kt ! time step
324	INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence
325	!
326	INTEGER :: ji, jj, jk ! dummy loop indices
327	INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers
328	REAL(wp) :: z2dt, z_tmin, z_tmax ! local scalars
329	LOGICAL :: ll_do_bclinic ! local logical
330	REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv
331	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t
332	!!----------------------------------------------------------------------
333	!
334	IF( ln_linssh ) RETURN ! No calculation in linear free surface
335	!
336	IF( ln_timing ) CALL timing_start('dom_vvl_sf_nxt')
337	!
338	IF( kt == nit000 ) THEN
339	IF(lwp) WRITE(numout,*)
340	IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
341	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
342	ENDIF
343
344	ll_do_bclinic = .TRUE.
345	IF( PRESENT(kcall) ) THEN
346	IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.
347	ENDIF
348
349	! ******************************* !
350	! After acale factors at t-points !
351	! ******************************* !
352	! ! --------------------------------------------- !
353	! ! z_star coordinate and barotropic z-tilde part !
354	! ! --------------------------------------------- !
355	!
356	z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
357	DO jk = 1, jpkm1
358	! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
359	e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
360	END DO
361	!
362	IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !
363	! ! ------baroclinic part------ !
364	! I - initialization
365	! ==================
366
367	! 1 - barotropic divergence
368	! -------------------------
369	zhdiv(:,:) = 0._wp
370	zht(:,:) = 0._wp
371	DO jk = 1, jpkm1
372	zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
373	zht (:,:) = zht (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
374	END DO
375	zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
376
377	! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only)
378	! --------------------------------------------------
379	IF( ln_vvl_ztilde ) THEN
380	IF( kt > nit000 ) THEN
381	DO jk = 1, jpkm1
382	hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) &
383	& * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
384	END DO
385	ENDIF
386	ENDIF
387
388	! II - after z_tilde increments of vertical scale factors
389	! =======================================================
390	tilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms
391
392	! 1 - High frequency divergence term
393	! ----------------------------------
394	IF( ln_vvl_ztilde ) THEN ! z_tilde case
395	DO jk = 1, jpkm1
396	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
397	END DO
398	ELSE ! layer case
399	DO jk = 1, jpkm1
400	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
401	END DO
402	ENDIF
403
404	! 2 - Restoring term (z-tilde case only)
405	! ------------------
406	IF( ln_vvl_ztilde ) THEN
407	DO jk = 1, jpk
408	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
409	END DO
410	ENDIF
411
412	! 3 - Thickness diffusion term
413	! ----------------------------
414	zwu(:,:) = 0._wp
415	zwv(:,:) = 0._wp
416	DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes
417	DO jj = 1, jpjm1
418	DO ji = 1, fs_jpim1 ! vector opt.
419	un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) &
420	& * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) )
421	vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) &
422	& * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) )
423	zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
424	zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
425	END DO
426	END DO
427	END DO
428	DO jj = 1, jpj ! b - correction for last oceanic u-v points
429	DO ji = 1, jpi
430	un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
431	vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
432	END DO
433	END DO
434	DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes
435	DO jj = 2, jpjm1
436	DO ji = fs_2, fs_jpim1 ! vector opt.
437	tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) &
438	& + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) &
439	& ) * r1_e1e2t(ji,jj)
440	END DO
441	END DO
442	END DO
443	! ! d - thickness diffusion transport: boundary conditions
444	! (stored for tracer advction and continuity equation)
445	CALL lbc_lnk_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
446
447	! 4 - Time stepping of baroclinic scale factors
448	! ---------------------------------------------
449	! Leapfrog time stepping
450	! ~~~~~~~~~~~~~~~~~~~~~~
451	IF( neuler == 0 .AND. kt == nit000 ) THEN
452	z2dt = rdt
453	ELSE
454	z2dt = 2.0_wp * rdt
455	ENDIF
456	CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
457	tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
458
459	! Maximum deformation control
460	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
461	ze3t(:,:,jpk) = 0._wp
462	DO jk = 1, jpkm1
463	ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
464	END DO
465	z_tmax = MAXVAL( ze3t(:,:,:) )
466	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
467	z_tmin = MINVAL( ze3t(:,:,:) )
468	CALL mpp_min( 'domvvl', z_tmin ) ! min over the global domain
469	! - ML - test: for the moment, stop simulation for too large e3_t variations
470	IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
471	IF( lk_mpp ) THEN
472	CALL mpp_maxloc( 'domvvl', ze3t, tmask, z_tmax, ijk_max )
473	CALL mpp_minloc( 'domvvl', ze3t, tmask, z_tmin, ijk_min )
474	ELSE
475	ijk_max = MAXLOC( ze3t(:,:,:) )
476	ijk_max(1) = ijk_max(1) + nimpp - 1
477	ijk_max(2) = ijk_max(2) + njmpp - 1
478	ijk_min = MINLOC( ze3t(:,:,:) )
479	ijk_min(1) = ijk_min(1) + nimpp - 1
480	ijk_min(2) = ijk_min(2) + njmpp - 1
481	ENDIF
482	IF (lwp) THEN
483	WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
484	WRITE(numout, *) 'at i, j, k=', ijk_max
485	WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
486	WRITE(numout, *) 'at i, j, k=', ijk_min
487	CALL ctl_stop( 'STOP', 'MAX( ABS( tilde_e3t_a(:,:,: ) ) / e3t_0(:,:,:) ) too high')
488	ENDIF
489	ENDIF
490	! - ML - end test
491	! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
492	tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) )
493	tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
494
495	!
496	! "tilda" change in the after scale factor
497	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
498	DO jk = 1, jpkm1
499	dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
500	END DO
501	! III - Barotropic repartition of the sea surface height over the baroclinic profile
502	! ==================================================================================
503	! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
504	! - ML - baroclinicity error should be better treated in the future
505	! i.e. locally and not spread over the water column.
506	! (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
507	zht(:,:) = 0.
508	DO jk = 1, jpkm1
509	zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
510	END DO
511	z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
512	DO jk = 1, jpkm1
513	dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
514	END DO
515
516	ENDIF
517
518	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate !
519	! ! ---baroclinic part--------- !
520	DO jk = 1, jpkm1
521	e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
522	END DO
523	ENDIF
524
525	IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN ! - ML - test: control prints for debuging
526	!
527	IF( lwp ) WRITE(numout, *) 'kt =', kt
528	IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
529	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
530	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
531	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
532	END IF
533	!
534	zht(:,:) = 0.0_wp
535	DO jk = 1, jpkm1
536	zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
537	END DO
538	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
539	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
540	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax
541	!
542	zht(:,:) = 0.0_wp
543	DO jk = 1, jpkm1
544	zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
545	END DO
546	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
547	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
548	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax
549	!
550	zht(:,:) = 0.0_wp
551	DO jk = 1, jpkm1
552	zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
553	END DO
554	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
555	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
556	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax
557	!
558	z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) )
559	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
560	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
561	!
562	z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshn(:,:) ) )
563	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
564	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
565	!
566	z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssha(:,:) ) )
567	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
568	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
569	END IF
570
571	! *********************************** !
572	! After scale factors at u- v- points !
573	! *********************************** !
574
575	CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
576	CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
577
578	! *********************************** !
579	! After depths at u- v points !
580	! *********************************** !
581
582	hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
583	hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
584	DO jk = 2, jpkm1
585	hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
586	hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
587	END DO
588	! ! Inverse of the local depth
589	!!gm BUG ? don't understand the use of umask_i here .....
590	r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
591	r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
592	!
593	IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt')
594	!
595	END SUBROUTINE dom_vvl_sf_nxt
596
597
598	SUBROUTINE dom_vvl_sf_swp( kt )
599	!!----------------------------------------------------------------------
600	!! * ROUTINE dom_vvl_sf_swp *
601	!!
602	!! ** Purpose : compute time filter and swap of scale factors
603	!! compute all depths and related variables for next time step
604	!! write outputs and restart file
605	!!
606	!! ** Method : - swap of e3t with trick for volume/tracer conservation
607	!! - reconstruct scale factor at other grid points (interpolate)
608	!! - recompute depths and water height fields
609	!!
610	!! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
611	!! - Recompute:
612	!! e3(u/v)_b
613	!! e3w_n
614	!! e3(u/v)w_b
615	!! e3(u/v)w_n
616	!! gdept_n, gdepw_n and gde3w_n
617	!! h(u/v) and h(u/v)r
618	!!
619	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
620	!! Leclair, M., and G. Madec, 2011, Ocean Modelling.
621	!!----------------------------------------------------------------------
622	INTEGER, INTENT( in ) :: kt ! time step
623	!
624	INTEGER :: ji, jj, jk ! dummy loop indices
625	REAL(wp) :: zcoef ! local scalar
626	!!----------------------------------------------------------------------
627	!
628	IF( ln_linssh ) RETURN ! No calculation in linear free surface
629	!
630	IF( ln_timing ) CALL timing_start('dom_vvl_sf_swp')
631	!
632	IF( kt == nit000 ) THEN
633	IF(lwp) WRITE(numout,*)
634	IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
635	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step'
636	ENDIF
637	!
638	! Time filter and swap of scale factors
639	! =====================================
640	! - ML - e3(t/u/v)_b are allready computed in dynnxt.
641	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
642	IF( neuler == 0 .AND. kt == nit000 ) THEN
643	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
644	ELSE
645	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &
646	& + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
647	ENDIF
648	tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
649	ENDIF
650	gdept_b(:,:,:) = gdept_n(:,:,:)
651	gdepw_b(:,:,:) = gdepw_n(:,:,:)
652
653	e3t_n(:,:,:) = e3t_a(:,:,:)
654	e3u_n(:,:,:) = e3u_a(:,:,:)
655	e3v_n(:,:,:) = e3v_a(:,:,:)
656
657	! Compute all missing vertical scale factor and depths
658	! ====================================================
659	! Horizontal scale factor interpolations
660	! --------------------------------------
661	! - ML - e3u_b and e3v_b are allready computed in dynnxt
662	! - JC - hu_b, hv_b, hur_b, hvr_b also
663
664	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
665
666	! Vertical scale factor interpolations
667	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' )
668	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
669	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
670	CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b(:,:,:), 'W' )
671	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
672	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
673
674	! t- and w- points depth (set the isf depth as it is in the initial step)
675	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
676	gdepw_n(:,:,1) = 0.0_wp
677	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
678	DO jk = 2, jpk
679	DO jj = 1,jpj
680	DO ji = 1,jpi
681	! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
682	! 1 for jk = mikt
683	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
684	gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
685	gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) &
686	& + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) )
687	gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
688	END DO
689	END DO
690	END DO
691
692	! Local depth and Inverse of the local depth of the water
693	! -------------------------------------------------------
694	hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:)
695	hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:)
696	!
697	ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
698	DO jk = 2, jpkm1
699	ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
700	END DO
701
702	! write restart file
703	! ==================
704	IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
705	!
706	IF( ln_timing ) CALL timing_stop('dom_vvl_sf_swp')
707	!
708	END SUBROUTINE dom_vvl_sf_swp
709
710
711	SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
712	!!---------------------------------------------------------------------
713	!! * ROUTINE dom_vvl__interpol *
714	!!
715	!! ** Purpose : interpolate scale factors from one grid point to another
716	!!
717	!! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0)
718	!! - horizontal interpolation: grid cell surface averaging
719	!! - vertical interpolation: simple averaging
720	!!----------------------------------------------------------------------
721	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3_in ! input e3 to be interpolated
722	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3_out ! output interpolated e3
723	CHARACTER(LEN=*) , INTENT(in ) :: pout ! grid point of out scale factors
724	! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'
725	!
726	INTEGER :: ji, jj, jk ! dummy loop indices
727	INTEGER :: iku, ikum1, ikv, ikvm1, ikf, ikfm1 !
728	REAL(wp) :: zlnwd ! =1./0. when ln_wd_il = T/
729	REAL(wp), DIMENSION(jpi,jpj) :: zssh ! work array to retrieve ssh (nn_vvl_interp > 1)
730	!!----------------------------------------------------------------------
731	!
732	IF(ln_wd_il) THEN
733	zlnwd = 1.0_wp
734	ELSE
735	zlnwd = 0.0_wp
736	END IF
737	!
738	SELECT CASE ( pout ) !== type of interpolation ==!
739	!
740	CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean
741	SELECT CASE ( nn_vvl_interp )
742	CASE ( 0 )
743	!
744	DO jk = 1, jpk
745	DO jj = 1, jpjm1
746	DO ji = 1, fs_jpim1 ! vector opt.
747	pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) &
748	& * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
749	& + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
750	END DO
751	END DO
752	END DO
753	!
754	CASE ( 1 )
755	!
756	DO jk = 1, jpk
757	DO jj = 1, jpjm1
758	DO ji = 1, fs_jpim1 ! vector opt.
759	pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) &
760	& * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
761	& + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
762	END DO
763	END DO
764	END DO
765	!
766	! Bottom correction:
767	DO jj = 1, jpjm1
768	DO ji = 1, fs_jpim1 ! vector opt.
769	iku = mbku(ji ,jj)
770	ikum1 = iku - 1
771	pe3_out(ji,jj,iku) = ( umask(ji,jj,iku) * (1.0_wp - zlnwd) + zlnwd ) &
772	& ( 0.5_wp r1_e1e2u(ji,jj) &
773	& ( e1e2t(ji ,jj) ( SUM(tmask(ji ,jj,:)*(pe3_in(ji ,jj,:) - e3t_0(ji ,jj,:))) ) &
774	& + e1e2t(ji+1,jj) * ( SUM(tmask(ji+1,jj,:)*(pe3_in(ji+1,jj,:) - e3t_0(ji+1,jj,:))) ) ) &
775	& - SUM(pe3_out(ji,jj,1:ikum1)))
776	END DO
777	END DO
778	!
779	CASE ( 2 )
780	zssh(:,:) = SUM(tmask(:,:,:)*(pe3_in(:,:,:)-e3t_0(:,:,:)), DIM=3)
781	DO jk = 1, jpk
782	DO jj = 1, jpjm1
783	DO ji = 1, fs_jpim1 ! vector opt.
784	pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) &
785	& * ( e1e2t(ji ,jj) * zssh(ji ,jj) + e1e2t(ji+1,jj) * zssh(ji+1,jj)) &
786	& * e3u_0(ji,jj,jk) / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) )
787	END DO
788	END DO
789	END DO
790	!
791	END SELECT
792	!
793	CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
794	pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
795	!
796	CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean
797	SELECT CASE ( nn_vvl_interp )
798	CASE ( 0 )
799	!
800	DO jk = 1, jpk
801	DO jj = 1, jpjm1
802	DO ji = 1, fs_jpim1 ! vector opt.
803	pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) &
804	& * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
805	& + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
806	END DO
807	END DO
808	END DO
809	! !
810	CASE ( 1 )
811	!
812	DO jk = 1, jpk
813	DO jj = 1, jpjm1
814	DO ji = 1, fs_jpim1 ! vector opt.
815	pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) &
816	& * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
817	& + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
818	END DO
819	END DO
820	END DO
821	!
822	! Bottom correction:
823	DO jj = 1, jpjm1
824	DO ji = 1, fs_jpim1 ! vector opt.
825	ikv = mbkv(ji ,jj)
826	ikvm1 = ikv - 1
827	pe3_out(ji,jj,ikv) = ( vmask(ji,jj,ikv) * (1.0_wp - zlnwd) + zlnwd ) &
828	& ( 0.5_wp r1_e1e2v(ji,jj) &
829	& ( e1e2t(ji,jj ) ( SUM(tmask(ji,jj ,:)*(pe3_in(ji,jj ,:) - e3t_0(ji,jj ,:))) ) &
830	& + e1e2t(ji,jj+1) * ( SUM(tmask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3t_0(ji,jj+1,:))) ) ) &
831	& - SUM(pe3_out(ji,jj,1:ikvm1)))
832	END DO
833	END DO
834	!
835	CASE ( 2 )
836	zssh(:,:) = SUM(tmask(:,:,:)*(pe3_in(:,:,:)-e3t_0(:,:,:)), DIM=3)
837	DO jk = 1, jpk
838	DO jj = 1, jpjm1
839	DO ji = 1, fs_jpim1 ! vector opt.
840	pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) &
841	& * ( e1e2t(ji ,jj) * zssh(ji ,jj) + e1e2t(ji,jj+1) * zssh(ji,jj+1)) &
842	& * e3v_0(ji,jj,jk) / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) )
843	END DO
844	END DO
845	END DO
846	!
847	END SELECT
848	!
849	CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
850	pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
851	!
852	CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean
853	SELECT CASE ( nn_vvl_interp )
854	CASE ( 0 )
855	DO jk = 1, jpk
856	DO jj = 1, jpjm1
857	DO ji = 1, fs_jpim1 ! vector opt.
858	pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
859	& * r1_e1e2f(ji,jj) &
860	& * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) &
861	& + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
862	END DO
863	END DO
864	END DO
865	!
866	CASE ( 1 )
867	!
868	DO jk = 1, jpk
869	DO jj = 1, jpjm1
870	DO ji = 1, fs_jpim1 ! vector opt.
871	pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
872	& * r1_e1e2f(ji,jj) &
873	& * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) &
874	& + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
875	END DO
876	END DO
877	END DO
878	!
879	! Bottom correction:
880	DO jj = 1, jpjm1
881	DO ji = 1, fs_jpim1 ! vector opt.
882	ikf = MIN(mbku(ji ,jj),mbku(ji,jj+1))
883	ikfm1 = ikf - 1
884	pe3_out(ji,jj,ikf) = ( umask(ji,jj,ikf) * umask(ji,jj+1,ikf) * (1.0_wp - zlnwd) + zlnwd ) &
885	& * ( 0.5_wp * r1_e1e2f(ji,jj) &
886	& * ( e1e2u(ji,jj ) * ( SUM(umask(ji,jj ,:)*(pe3_in(ji,jj ,:) - e3u_0(ji,jj ,:))) ) &
887	& + e1e2u(ji,jj+1) * ( SUM(umask(ji,jj+1,:)*(pe3_in(ji,jj+1,:) - e3u_0(ji,jj+1,:))) ) ) &
888	& - SUM(pe3_out(ji,jj,1:ikfm1)))
889	END DO
890	END DO
891	!
892	CASE ( 2 )
893	zssh(:,:) = SUM(umask(:,:,:)*(pe3_in(:,:,:)-e3u_0(:,:,:)), DIM=3)
894	DO jk = 1, jpk
895	DO jj = 1, jpjm1
896	DO ji = 1, fs_jpim1 ! vector opt.
897	pe3_out(ji,jj,jk) = ( umask(ji,jj,jk)* umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
898	& * 0.5_wp * r1_e1e2f(ji,jj) &
899	& * (e1e2u(ji ,jj) * zssh(ji ,jj) + e1e2u(ji,jj+1) * zssh(ji,jj+1)) &
900	& * e3f_0(ji,jj,jk) / ( hf_0(ji,jj) + 1._wp - ssumask(ji,jj)*ssumask(ji,jj+1) )
901	END DO
902	END DO
903	END DO
904	!
905	END SELECT
906	!
907	CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
908	pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
909	!
910	CASE( 'W' ) !* from T- to W-point : vertical simple mean
911	!
912	pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
913	! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
914	!!gm BUG? use here wmask in case of ISF ? to be checked
915	DO jk = 2, jpk
916	pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
917	& * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) &
918	& + 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
919	& * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) )
920	END DO
921	!
922	CASE( 'UW' ) !* from U- to UW-point : vertical simple mean
923	!
924	pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
925	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
926	!!gm BUG? use here wumask in case of ISF ? to be checked
927	DO jk = 2, jpk
928	pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
929	& * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) &
930	& + 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
931	& * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) )
932	END DO
933	!
934	CASE( 'VW' ) !* from V- to VW-point : vertical simple mean
935	!
936	pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
937	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
938	!!gm BUG? use here wvmask in case of ISF ? to be checked
939	DO jk = 2, jpk
940	pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
941	& * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) &
942	& + 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
943	& * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) )
944	END DO
945	END SELECT
946	!
947	END SUBROUTINE dom_vvl_interpol
948
949
950	SUBROUTINE dom_vvl_rst( kt, cdrw )
951	!!---------------------------------------------------------------------
952	!! * ROUTINE dom_vvl_rst *
953	!!
954	!! ** Purpose : Read or write VVL file in restart file
955	!!
956	!! ** Method : use of IOM library
957	!! if the restart does not contain vertical scale factors,
958	!! they are set to the _0 values
959	!! if the restart does not contain vertical scale factors increments (z_tilde),
960	!! they are set to 0.
961	!!----------------------------------------------------------------------
962	INTEGER , INTENT(in) :: kt ! ocean time-step
963	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
964	!
965	INTEGER :: ji, jj, jk
966	INTEGER :: id1, id2, id3, id4, id5 ! local integers
967	!!----------------------------------------------------------------------
968	!
969	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
970	! ! ===============
971	IF( ln_rstart ) THEN !* Read the restart file
972	CALL rst_read_open ! open the restart file if necessary
973	CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn, ldxios = lrxios )
974	!
975	id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
976	id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
977	id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
978	id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
979	id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
980	! ! --------- !
981	! ! all cases !
982	! ! --------- !
983	IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist
984	CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
985	CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
986	! needed to restart if land processor not computed
987	IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
988	WHERE ( tmask(:,:,:) == 0.0_wp )
989	e3t_n(:,:,:) = e3t_0(:,:,:)
990	e3t_b(:,:,:) = e3t_0(:,:,:)
991	END WHERE
992	IF( neuler == 0 ) THEN
993	e3t_b(:,:,:) = e3t_n(:,:,:)
994	ENDIF
995	ELSE IF( id1 > 0 ) THEN
996	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
997	IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
998	IF(lwp) write(numout,*) 'neuler is forced to 0'
999	CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
1000	e3t_n(:,:,:) = e3t_b(:,:,:)
1001	neuler = 0
1002	ELSE IF( id2 > 0 ) THEN
1003	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
1004	IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
1005	IF(lwp) write(numout,*) 'neuler is forced to 0'
1006	CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
1007	e3t_b(:,:,:) = e3t_n(:,:,:)
1008	neuler = 0
1009	ELSE
1010	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
1011	IF(lwp) write(numout,*) 'Compute scale factor from sshn'
1012	IF(lwp) write(numout,*) 'neuler is forced to 0'
1013	DO jk = 1, jpk
1014	e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
1015	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
1016	& + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
1017	END DO
1018	e3t_b(:,:,:) = e3t_n(:,:,:)
1019	neuler = 0
1020	ENDIF
1021	! ! ----------- !
1022	IF( ln_vvl_zstar ) THEN ! z_star case !
1023	! ! ----------- !
1024	IF( MIN( id3, id4 ) > 0 ) THEN
1025	CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
1026	ENDIF
1027	! ! ----------------------- !
1028	ELSE ! z_tilde and layer cases !
1029	! ! ----------------------- !
1030	IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist
1031	CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )
1032	CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )
1033	ELSE ! one at least array is missing
1034	tilde_e3t_b(:,:,:) = 0.0_wp
1035	tilde_e3t_n(:,:,:) = 0.0_wp
1036	ENDIF
1037	! ! ------------ !
1038	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
1039	! ! ------------ !
1040	IF( id5 > 0 ) THEN ! required array exists
1041	CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )
1042	ELSE ! array is missing
1043	hdiv_lf(:,:,:) = 0.0_wp
1044	ENDIF
1045	ENDIF
1046	ENDIF
1047	!
1048	ELSE !* Initialize at "rest"
1049	!
1050
1051	IF( ll_wd ) THEN ! MJB ll_wd edits start here - these are essential
1052	!
1053	IF( cn_cfg == 'wad' ) THEN
1054	! Wetting and drying test case
1055	CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb )
1056	tsn (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones
1057	sshn (:,:) = sshb(:,:)
1058	un (:,:,:) = ub (:,:,:)
1059	vn (:,:,:) = vb (:,:,:)
1060	ELSE
1061	! if not test case
1062	sshn(:,:) = -ssh_ref
1063	sshb(:,:) = -ssh_ref
1064
1065	DO jj = 1, jpj
1066	DO ji = 1, jpi
1067	IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth
1068
1069	sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
1070	sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
1071	ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
1072	ENDIF
1073	ENDDO
1074	ENDDO
1075	ENDIF !If test case else
1076
1077	! Adjust vertical metrics for all wad
1078	DO jk = 1, jpk
1079	e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
1080	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
1081	& + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
1082	END DO
1083	e3t_b(:,:,:) = e3t_n(:,:,:)
1084
1085	DO ji = 1, jpi
1086	DO jj = 1, jpj
1087	IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
1088	CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
1089	ENDIF
1090	END DO
1091	END DO
1092	!
1093	ELSE
1094	!
1095	! usr_def_istate called here only to get sshb, that is needed to
1096	! initialize e3t_b and e3t_n
1097	CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb )
1098	! usr_def_istate will be called again in istate_init to
1099	! initialize ts(bn), ssh(bn), u(bn) and v(bn)
1100	!
1101	DO jk=1,jpk
1102	e3t_b(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &
1103	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
1104	& + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t_b != 0 on land points
1105	END DO
1106	e3t_n(:,:,:) = e3t_b(:,:,:)
1107	sshn(:,:) = sshb(:,:) ! needed later for gde3w
1108	!!$ e3t_n(:,:,:)=e3t_0(:,:,:)
1109	!!$ e3t_b(:,:,:)=e3t_0(:,:,:)
1110	!
1111	END IF ! end of ll_wd edits
1112
1113	IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
1114	tilde_e3t_b(:,:,:) = 0._wp
1115	tilde_e3t_n(:,:,:) = 0._wp
1116	IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
1117	END IF
1118	ENDIF
1119	!
1120	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
1121	! ! ===================
1122	IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
1123	IF( lwxios ) CALL iom_swap( cwxios_context )
1124	! ! --------- !
1125	! ! all cases !
1126	! ! --------- !
1127	CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios )
1128	CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )
1129	! ! ----------------------- !
1130	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
1131	! ! ----------------------- !
1132	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
1133	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
1134	END IF
1135	! ! -------------!
1136	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
1137	! ! ------------ !
1138	CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
1139	ENDIF
1140	!
1141	IF( lwxios ) CALL iom_swap( cxios_context )
1142	ENDIF
1143	!
1144	END SUBROUTINE dom_vvl_rst
1145
1146
1147	SUBROUTINE dom_vvl_ctl
1148	!!---------------------------------------------------------------------
1149	!! * ROUTINE dom_vvl_ctl *
1150	!!
1151	!! ** Purpose : Control the consistency between namelist options
1152	!! for vertical coordinate
1153	!!----------------------------------------------------------------------
1154	INTEGER :: ioptio, ios
1155	!!
1156	NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
1157	& ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , &
1158	& rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg , &! not yet implemented: ln_vvl_kepe
1159	& nn_vvl_interp
1160	!!----------------------------------------------------------------------
1161	!
1162	REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist :
1163	READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
1164	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
1165	REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run
1166	READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
1167	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
1168	IF(lwm) WRITE ( numond, nam_vvl )
1169	!
1170	IF(lwp) THEN ! Namelist print
1171	WRITE(numout,*)
1172	WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
1173	WRITE(numout,*) '~~~~~~~~~~~'
1174	WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate'
1175	WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar
1176	WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde
1177	WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer
1178	WRITE(numout,*) ' ztilde as zstar ln_vvl_ztilde_as_zstar = ', ln_vvl_ztilde_as_zstar
1179	WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor
1180	WRITE(numout,*) ' !'
1181	WRITE(numout,*) ' thickness diffusion coefficient rn_ahe3 = ', rn_ahe3
1182	WRITE(numout,*) ' maximum e3t deformation fractional change rn_zdef_max = ', rn_zdef_max
1183	IF( ln_vvl_ztilde_as_zstar ) THEN
1184	WRITE(numout,*) ' ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
1185	WRITE(numout,*) ' ignoring namelist timescale parameters and using:'
1186	WRITE(numout,*) ' hard-wired : z-tilde to zstar restoration timescale (days)'
1187	WRITE(numout,*) ' rn_rst_e3t = 0.e0'
1188	WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1189	WRITE(numout,*) ' rn_lf_cutoff = 1.0/rdt'
1190	ELSE
1191	WRITE(numout,*) ' z-tilde to zstar restoration timescale (days) rn_rst_e3t = ', rn_rst_e3t
1192	WRITE(numout,*) ' z-tilde cutoff frequency of low-pass filter (days) rn_lf_cutoff = ', rn_lf_cutoff
1193	ENDIF
1194	WRITE(numout,*) ' debug prints flag ln_vvl_dbg = ', ln_vvl_dbg
1195	WRITE(numout,*) ' Method to compute scale factors anomaly at U/V/F points nn_vvl_interp = ', nn_vvl_interp
1196	ENDIF
1197	!
1198	ioptio = 0 ! Parameter control
1199	IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
1200	IF( ln_vvl_zstar ) ioptio = ioptio + 1
1201	IF( ln_vvl_ztilde ) ioptio = ioptio + 1
1202	IF( ln_vvl_layer ) ioptio = ioptio + 1
1203	!
1204	IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1205	IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1206	!
1207	IF( .NOT. ln_vvl_zstar .AND. (nn_vvl_interp==2 ) ) CALL ctl_stop( 'nn_vvl_interp must be < 2 if ln_vvl_zstar=F' )
1208	IF(lwp) THEN ! Print the choice
1209	WRITE(numout,*)
1210	IF( ln_vvl_zstar ) WRITE(numout,*) ' ==>>> zstar vertical coordinate is used'
1211	IF( ln_vvl_ztilde ) WRITE(numout,*) ' ==>>> ztilde vertical coordinate is used'
1212	IF( ln_vvl_layer ) WRITE(numout,*) ' ==>>> layer vertical coordinate is used'
1213	IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) ' ==>>> to emulate a zstar coordinate'
1214	ENDIF
1215	!
1216	#if defined key_agrif
1217	IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) ) CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1218	#endif
1219	!
1220	END SUBROUTINE dom_vvl_ctl
1221
1222	!!======================================================================
1223	END MODULE domvvl

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

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

Context Navigation

source: NEMO/releases/r4.0/r4.0-HEAD/tests/CANAL/MY_SRC/domvvl.F90 @ 15610

Download in other formats: