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 trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 7698

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

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

Download in other formats: