Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

domvvl.F90 @ 8011

Last change on this file since 8011 was 8011, checked in by jchanut, 7 years ago
VORTEX test case for AGRIF
File size: 76.7 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	!! 'key_vvl' variable volume
13	!!----------------------------------------------------------------------
14	!!----------------------------------------------------------------------
15	!! dom_vvl_init : define initial vertical scale factors, depths and column thickness
16	!! dom_vvl_sf_nxt : Compute next vertical scale factors
17	!! dom_vvl_sf_swp : Swap vertical scale factors and update the vertical grid
18	!! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
19	!! dom_vvl_rst : read/write restart file
20	!! dom_vvl_ctl : Check the vvl options
21	!! dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors
22	!! : to account for manual changes to e[1,2][u,v] in some Straits
23	!!----------------------------------------------------------------------
24	!! * Modules used
25	USE oce ! ocean dynamics and tracers
26	USE dom_oce ! ocean space and time domain
27	USE sbc_oce ! ocean surface boundary condition
28	USE in_out_manager ! I/O manager
29	USE iom ! I/O manager library
30	USE restart ! ocean restart
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	!! * Routine accessibility
40	PUBLIC dom_vvl_init ! called by domain.F90
41	PUBLIC dom_vvl_sf_nxt ! called by step.F90
42	PUBLIC dom_vvl_sf_swp ! called by step.F90
43	PUBLIC dom_vvl_interpol ! called by dynnxt.F90
44	PRIVATE dom_vvl_orca_fix ! called by dom_vvl_interpol
45
46	!!* Namelist nam_vvl
47	LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate
48	LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate
49	LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate
50	LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate
51	LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate
52	LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer
53	! ! conservation: not used yet
54	REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient
55	REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days]
56	REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days]
57	REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation
58	LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints
59
60	!! * Module variables
61	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport
62	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence
63	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors
64	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors
65	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors
66	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence
67
68	!! * Substitutions
69	# include "domzgr_substitute.h90"
70	# include "vectopt_loop_substitute.h90"
71	!!----------------------------------------------------------------------
72	!! NEMO/OPA 3.3 , NEMO-Consortium (2010)
73	!! $Id: domvvl.F90 7971 2017-04-26 10:03:26Z jchanut $
74	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
75	!!----------------------------------------------------------------------
76
77	CONTAINS
78
79	INTEGER FUNCTION dom_vvl_alloc()
80	!!----------------------------------------------------------------------
81	!! * FUNCTION dom_vvl_alloc *
82	!!----------------------------------------------------------------------
83	IF( ln_vvl_zstar ) dom_vvl_alloc = 0
84	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
85	ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , &
86	& dtilde_e3t_a(jpi,jpj,jpk) , un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , &
87	& STAT = dom_vvl_alloc )
88	IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc )
89	IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
90	un_td = 0.0_wp
91	vn_td = 0.0_wp
92	ENDIF
93	IF( ln_vvl_ztilde ) THEN
94	ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
95	IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc )
96	IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
97	ENDIF
98
99	END FUNCTION dom_vvl_alloc
100
101
102	SUBROUTINE dom_vvl_init
103	!!----------------------------------------------------------------------
104	!! * ROUTINE dom_vvl_init *
105	!!
106	!! ** Purpose : Initialization of all scale factors, depths
107	!! and water column heights
108	!!
109	!! ** Method : - use restart file and/or initialize
110	!! - interpolate scale factors
111	!!
112	!! ** Action : - fse3t_(n/b) and tilde_e3t_(n/b)
113	!! - Regrid: fse3(u/v)_n
114	!! fse3(u/v)_b
115	!! fse3w_n
116	!! fse3(u/v)w_b
117	!! fse3(u/v)w_n
118	!! fsdept_n, fsdepw_n and fsde3w_n
119	!! - h(t/u/v)_0
120	!! - frq_rst_e3t and frq_rst_hdv
121	!!
122	!! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling.
123	!!----------------------------------------------------------------------
124	USE phycst, ONLY : rpi, rsmall, rad
125	!! * Local declarations
126	INTEGER :: ji,jj,jk
127	INTEGER :: ii0, ii1, ij0, ij1
128	REAL(wp):: zcoef
129	!!----------------------------------------------------------------------
130	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init')
131
132	IF(lwp) WRITE(numout,*)
133	IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
134	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
135
136	! choose vertical coordinate (z_star, z_tilde or layer)
137	! ==========================
138	CALL dom_vvl_ctl
139
140	! Allocate module arrays
141	! ======================
142	IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
143
144	! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk))
145	! ============================================================================
146	CALL dom_vvl_rst( nit000, 'READ' )
147	fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
148
149	! Reconstruction of all vertical scale factors at now and before time steps
150	! =============================================================================
151	! Horizontal scale factor interpolations
152	! --------------------------------------
153	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
154	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
155	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
156	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
157	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
158	! Vertical scale factor interpolations
159	! ------------------------------------
160	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )
161	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
162	CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
163	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' )
164	CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
165	CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
166	! t- and w- points depth
167	! ----------------------
168	! set the isf depth as it is in the initial step
169	fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
170	fsdepw_n(:,:,1) = 0.0_wp
171	fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
172	fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1)
173	fsdepw_b(:,:,1) = 0.0_wp
174
175	DO jk = 2, jpk
176	DO jj = 1,jpj
177	DO ji = 1,jpi
178	! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
179	! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
180	! 0.5 where jk = mikt
181	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
182	fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
183	fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) &
184	& + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))
185	fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
186	fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1)
187	fsdept_b(ji,jj,jk) = zcoef * ( fsdepw_b(ji,jj,jk ) + 0.5 * fse3w_b(ji,jj,jk)) &
188	& + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk))
189	END DO
190	END DO
191	END DO
192
193	! Before depth and Inverse of the local depth of the water column at u- and v- points
194	! -----------------------------------------------------------------------------------
195	hu_b(:,:) = 0.
196	hv_b(:,:) = 0.
197	DO jk = 1, jpkm1
198	hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
199	hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
200	END DO
201	hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) )
202	hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) )
203
204	! Restoring frequencies for z_tilde coordinate
205	! ============================================
206	IF( ln_vvl_ztilde ) THEN
207	! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings
208	frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp )
209	frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
210	IF( ln_vvl_ztilde_as_zstar ) THEN
211	! Ignore namelist settings and use these next two to emulate z-star using z-tilde
212	frq_rst_e3t(:,:) = 0.0_wp
213	frq_rst_hdv(:,:) = 1.0_wp / rdt
214	ENDIF
215	IF ( ln_vvl_zstar_at_eqtor ) THEN
216	DO jj = 1, jpj
217	DO ji = 1, jpi
218	IF( ABS(gphit(ji,jj)) >= 6.) THEN
219	! values outside the equatorial band and transition zone (ztilde)
220	frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp )
221	frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
222	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
223	! values inside the equatorial band (ztilde as zstar)
224	frq_rst_e3t(ji,jj) = 0.0_wp
225	frq_rst_hdv(ji,jj) = 1.0_wp / rdt
226	ELSE
227	! values in the transition band (linearly vary from ztilde to ztilde as zstar values)
228	frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp &
229	& * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
230	& * 180._wp / 3.5_wp ) )
231	frq_rst_hdv(ji,jj) = (1.0_wp / rdt) &
232	& + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp &
233	& * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
234	& * 180._wp / 3.5_wp ) )
235	ENDIF
236	END DO
237	END DO
238	IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN
239	ii0 = 103 ; ii1 = 111 ! Suppress ztilde in the Foxe Basin for ORCA2
240	ij0 = 128 ; ij1 = 135 ;
241	frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp
242	frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rdt
243	ENDIF
244	ENDIF
245	ENDIF
246
247	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_init')
248
249	END SUBROUTINE dom_vvl_init
250
251
252	SUBROUTINE dom_vvl_sf_nxt( kt, kcall )
253	!!----------------------------------------------------------------------
254	!! * ROUTINE dom_vvl_sf_nxt *
255	!!
256	!! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt,
257	!! tranxt and dynspg routines
258	!!
259	!! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness.
260	!! - z_tilde_case: after scale factor increment =
261	!! high frequency part of horizontal divergence
262	!! + retsoring towards the background grid
263	!! + thickness difusion
264	!! Then repartition of ssh INCREMENT proportionnaly
265	!! to the "baroclinic" level thickness.
266	!!
267	!! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case
268	!! - tilde_e3t_a: after increment of vertical scale factor
269	!! in z_tilde case
270	!! - fse3(t/u/v)_a
271	!!
272	!! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling.
273	!!----------------------------------------------------------------------
274	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t
275	REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv
276	!! * Arguments
277	INTEGER, INTENT( in ) :: kt ! time step
278	INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence
279	!! * Local declarations
280	INTEGER :: ji, jj, jk ! dummy loop indices
281	INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers
282	REAL(wp) :: z2dt ! temporary scalars
283	REAL(wp) :: z_tmin, z_tmax ! temporary scalars
284	LOGICAL :: ll_do_bclinic ! temporary logical
285	!!----------------------------------------------------------------------
286	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt')
287	CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
288	CALL wrk_alloc( jpi, jpj, jpk, ze3t )
289
290	IF(kt == nit000) THEN
291	IF(lwp) WRITE(numout,*)
292	IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
293	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
294	ENDIF
295
296	ll_do_bclinic = .TRUE.
297	IF( PRESENT(kcall) ) THEN
298	IF ( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.
299	ENDIF
300
301	! ******************************* !
302	! After acale factors at t-points !
303	! ******************************* !
304
305	! ! --------------------------------------------- !
306	! z_star coordinate and barotropic z-tilde part !
307	! ! --------------------------------------------- !
308
309	z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
310	DO jk = 1, jpkm1
311	! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
312	fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
313	END DO
314
315	IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !
316	! ! ------baroclinic part------ !
317
318	! I - initialization
319	! ==================
320
321	! 1 - barotropic divergence
322	! -------------------------
323	zhdiv(:,:) = 0.
324	zht(:,:) = 0.
325	DO jk = 1, jpkm1
326	zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
327	zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
328	END DO
329	zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
330
331	! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only)
332	! --------------------------------------------------
333	IF( ln_vvl_ztilde ) THEN
334	IF( kt .GT. nit000 ) THEN
335	DO jk = 1, jpkm1
336	hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) &
337	& * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
338	END DO
339	ENDIF
340	END IF
341
342	! II - after z_tilde increments of vertical scale factors
343	! =======================================================
344	tilde_e3t_a(:,:,:) = 0.0_wp ! tilde_e3t_a used to store tendency terms
345
346	! 1 - High frequency divergence term
347	! ----------------------------------
348	IF( ln_vvl_ztilde ) THEN ! z_tilde case
349	DO jk = 1, jpkm1
350	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
351	END DO
352	ELSE ! layer case
353	DO jk = 1, jpkm1
354	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
355	END DO
356	END IF
357
358	! 2 - Restoring term (z-tilde case only)
359	! ------------------
360	IF( ln_vvl_ztilde ) THEN
361	DO jk = 1, jpk
362	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
363	END DO
364	END IF
365
366	! 3 - Thickness diffusion term
367	! ----------------------------
368	zwu(:,:) = 0.0_wp
369	zwv(:,:) = 0.0_wp
370	! a - first derivative: diffusive fluxes
371	DO jk = 1, jpkm1
372	DO jj = 1, jpjm1
373	DO ji = 1, fs_jpim1 ! vector opt.
374	un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) &
375	& * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) )
376	vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) &
377	& * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) )
378	zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
379	zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
380	END DO
381	END DO
382	END DO
383	! b - correction for last oceanic u-v points
384	DO jj = 1, jpj
385	DO ji = 1, jpi
386	un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
387	vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
388	END DO
389	END DO
390	! c - second derivative: divergence of diffusive fluxes
391	DO jk = 1, jpkm1
392	DO jj = 2, jpjm1
393	DO ji = fs_2, fs_jpim1 ! vector opt.
394	tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) &
395	& + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) &
396	& ) * r1_e12t(ji,jj)
397	END DO
398	END DO
399	END DO
400	! d - thickness diffusion transport: boundary conditions
401	! (stored for tracer advction and continuity equation)
402	CALL lbc_lnk( un_td , 'U' , -1._wp)
403	CALL lbc_lnk( vn_td , 'V' , -1._wp)
404
405	! 4 - Time stepping of baroclinic scale factors
406	! ---------------------------------------------
407	! Leapfrog time stepping
408	! ~~~~~~~~~~~~~~~~~~~~~~
409	IF( neuler == 0 .AND. kt == nit000 ) THEN
410	z2dt = rdt
411	ELSE
412	z2dt = 2.0_wp * rdt
413	ENDIF
414	CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )
415	tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
416
417	! Maximum deformation control
418	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
419	ze3t(:,:,jpk) = 0.0_wp
420	DO jk = 1, jpkm1
421	ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
422	END DO
423	z_tmax = MAXVAL( ze3t(:,:,:) )
424	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
425	z_tmin = MINVAL( ze3t(:,:,:) )
426	IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain
427	! - ML - test: for the moment, stop simulation for too large e3_t variations
428	IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN
429	IF( lk_mpp ) THEN
430	CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
431	CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
432	ELSE
433	ijk_max = MAXLOC( ze3t(:,:,:) )
434	ijk_max(1) = ijk_max(1) + nimpp - 1
435	ijk_max(2) = ijk_max(2) + njmpp - 1
436	ijk_min = MINLOC( ze3t(:,:,:) )
437	ijk_min(1) = ijk_min(1) + nimpp - 1
438	ijk_min(2) = ijk_min(2) + njmpp - 1
439	ENDIF
440	IF (lwp) THEN
441	WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
442	WRITE(numout, *) 'at i, j, k=', ijk_max
443	WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
444	WRITE(numout, *) 'at i, j, k=', ijk_min
445	CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
446	ENDIF
447	ENDIF
448	! - ML - end test
449	! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
450	tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) )
451	tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
452
453	!
454	! "tilda" change in the after scale factor
455	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
456	DO jk = 1, jpkm1
457	dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
458	END DO
459	! III - Barotropic repartition of the sea surface height over the baroclinic profile
460	! ==================================================================================
461	! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
462	! - ML - baroclinicity error should be better treated in the future
463	! i.e. locally and not spread over the water column.
464	! (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
465	zht(:,:) = 0.
466	DO jk = 1, jpkm1
467	zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
468	END DO
469	z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
470	DO jk = 1, jpkm1
471	dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
472	END DO
473
474	ENDIF
475
476	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate !
477	! ! ---baroclinic part--------- !
478	DO jk = 1, jpkm1
479	fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
480	END DO
481	ENDIF
482
483	IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN ! - ML - test: control prints for debuging
484	!
485	IF( lwp ) WRITE(numout, *) 'kt =', kt
486	IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
487	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
488	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
489	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
490	END IF
491	!
492	zht(:,:) = 0.0_wp
493	DO jk = 1, jpkm1
494	zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
495	END DO
496	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
497	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
498	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', z_tmax
499	!
500	zht(:,:) = 0.0_wp
501	DO jk = 1, jpkm1
502	zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
503	END DO
504	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
505	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
506	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', z_tmax
507	!
508	zht(:,:) = 0.0_wp
509	DO jk = 1, jpkm1
510	zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)
511	END DO
512	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
513	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
514	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(fse3t_b))) =', z_tmax
515	!
516	z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) )
517	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
518	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
519	!
520	z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshn(:,:) ) )
521	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
522	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
523	!
524	z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssha(:,:) ) )
525	IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain
526	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
527	END IF
528
529	! *********************************** !
530	! After scale factors at u- v- points !
531	! *********************************** !
532
533	CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )
534	CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
535
536	! *********************************** !
537	! After depths at u- v points !
538	! *********************************** !
539
540	hu_a(:,:) = 0._wp ! Ocean depth at U-points
541	hv_a(:,:) = 0._wp ! Ocean depth at V-points
542	DO jk = 1, jpkm1
543	hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
544	hv_a(:,:) = hv_a(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
545	END DO
546	! ! Inverse of the local depth
547	hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
548	hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
549
550	CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
551	CALL wrk_dealloc( jpi, jpj, jpk, ze3t )
552
553	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_nxt')
554
555	END SUBROUTINE dom_vvl_sf_nxt
556
557
558	SUBROUTINE dom_vvl_sf_swp( kt )
559	!!----------------------------------------------------------------------
560	!! * ROUTINE dom_vvl_sf_swp *
561	!!
562	!! ** Purpose : compute time filter and swap of scale factors
563	!! compute all depths and related variables for next time step
564	!! write outputs and restart file
565	!!
566	!! ** Method : - swap of e3t with trick for volume/tracer conservation
567	!! - reconstruct scale factor at other grid points (interpolate)
568	!! - recompute depths and water height fields
569	!!
570	!! ** Action : - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step
571	!! - Recompute:
572	!! fse3(u/v)_b
573	!! fse3w_n
574	!! fse3(u/v)w_b
575	!! fse3(u/v)w_n
576	!! fsdept_n, fsdepw_n and fsde3w_n
577	!! h(u/v) and h(u/v)r
578	!!
579	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
580	!! Leclair, M., and G. Madec, 2011, Ocean Modelling.
581	!!----------------------------------------------------------------------
582	!! * Arguments
583	INTEGER, INTENT( in ) :: kt ! time step
584	!! * Local declarations
585	INTEGER :: ji,jj,jk ! dummy loop indices
586	REAL(wp) :: zcoef
587	!!----------------------------------------------------------------------
588
589	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp')
590	!
591	IF( kt == nit000 ) THEN
592	IF(lwp) WRITE(numout,*)
593	IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
594	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step'
595	ENDIF
596	!
597	! Time filter and swap of scale factors
598	! =====================================
599	! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
600	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
601	IF( neuler == 0 .AND. kt == nit000 ) THEN
602	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
603	ELSE
604	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &
605	& + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
606	ENDIF
607	tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
608	ENDIF
609	fsdept_b(:,:,:) = fsdept_n(:,:,:)
610	fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
611
612	fse3t_n(:,:,:) = fse3t_a(:,:,:)
613	fse3u_n(:,:,:) = fse3u_a(:,:,:)
614	fse3v_n(:,:,:) = fse3v_a(:,:,:)
615
616	! Compute all missing vertical scale factor and depths
617	! ====================================================
618	! Horizontal scale factor interpolations
619	! --------------------------------------
620	! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
621	! - JC - hu_b, hv_b, hur_b, hvr_b also
622	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F' )
623	! Vertical scale factor interpolations
624	! ------------------------------------
625	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )
626	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
627	CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
628	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' )
629	CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
630	CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
631	! t- and w- points depth
632	! ----------------------
633	! set the isf depth as it is in the initial step
634	fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
635	fsdepw_n(:,:,1) = 0.0_wp
636	fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
637
638	DO jk = 2, jpk
639	DO jj = 1,jpj
640	DO ji = 1,jpi
641	! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
642	! 1 for jk = mikt
643	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
644	fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
645	fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) &
646	& + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))
647	fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
648	END DO
649	END DO
650	END DO
651
652	! Local depth and Inverse of the local depth of the water column at u- and v- points
653	! ----------------------------------------------------------------------------------
654	hu (:,:) = hu_a (:,:)
655	hv (:,:) = hv_a (:,:)
656
657	! Inverse of the local depth
658	hur(:,:) = hur_a(:,:)
659	hvr(:,:) = hvr_a(:,:)
660
661	! Local depth of the water column at t- points
662	! --------------------------------------------
663	ht(:,:) = 0.
664	DO jk = 1, jpkm1
665	ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
666	END DO
667
668	! write restart file
669	! ==================
670	IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
671	!
672	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp')
673
674	END SUBROUTINE dom_vvl_sf_swp
675
676
677	SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
678	!!---------------------------------------------------------------------
679	!! * ROUTINE dom_vvl__interpol *
680	!!
681	!! ** Purpose : interpolate scale factors from one grid point to another
682	!!
683	!! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0)
684	!! - horizontal interpolation: grid cell surface averaging
685	!! - vertical interpolation: simple averaging
686	!!----------------------------------------------------------------------
687	!! * Arguments
688	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated
689	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3
690	CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors
691	! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'
692	!! * Local declarations
693	INTEGER :: ji, jj, jk ! dummy loop indices
694	LOGICAL :: l_is_orca ! local logical
695	!!----------------------------------------------------------------------
696	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol')
697	!
698	l_is_orca = .FALSE.
699	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE. ! ORCA R2 configuration - will need to correct some locations
700
701	SELECT CASE ( pout )
702	! ! ------------------------------------- !
703	CASE( 'U' ) ! interpolation from T-point to U-point !
704	! ! ------------------------------------- !
705	! horizontal surface weighted interpolation
706	DO jk = 1, jpk
707	DO jj = 1, jpjm1
708	DO ji = 1, fs_jpim1 ! vector opt.
709	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) &
710	& * ( e12t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
711	& + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
712	END DO
713	END DO
714	END DO
715	!
716	IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
717	! boundary conditions
718	CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
719	pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
720	! ! ------------------------------------- !
721	CASE( 'V' ) ! interpolation from T-point to V-point !
722	! ! ------------------------------------- !
723	! horizontal surface weighted interpolation
724	DO jk = 1, jpk
725	DO jj = 1, jpjm1
726	DO ji = 1, fs_jpim1 ! vector opt.
727	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) &
728	& * ( e12t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
729	& + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
730	END DO
731	END DO
732	END DO
733	!
734	IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
735	! boundary conditions
736	CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
737	pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
738	! ! ------------------------------------- !
739	CASE( 'F' ) ! interpolation from U-point to F-point !
740	! ! ------------------------------------- !
741	! horizontal surface weighted interpolation
742	DO jk = 1, jpk
743	DO jj = 1, jpjm1
744	DO ji = 1, fs_jpim1 ! vector opt.
745	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) &
746	& * ( e12u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) &
747	& + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
748	END DO
749	END DO
750	END DO
751	!
752	IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
753	! boundary conditions
754	CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
755	pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
756	! ! ------------------------------------- !
757	CASE( 'W' ) ! interpolation from T-point to W-point !
758	! ! ------------------------------------- !
759	! vertical simple interpolation
760	pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
761	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
762	DO jk = 2, jpk
763	pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) &
764	& + 0.5_wp * tmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) )
765	END DO
766	! ! -------------------------------------- !
767	CASE( 'UW' ) ! interpolation from U-point to UW-point !
768	! ! -------------------------------------- !
769	! vertical simple interpolation
770	pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
771	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
772	DO jk = 2, jpk
773	pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) &
774	& + 0.5_wp * umask(:,:,jk) * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) )
775	END DO
776	! ! -------------------------------------- !
777	CASE( 'VW' ) ! interpolation from V-point to VW-point !
778	! ! -------------------------------------- !
779	! vertical simple interpolation
780	pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
781	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
782	DO jk = 2, jpk
783	pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) &
784	& + 0.5_wp * vmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) )
785	END DO
786	END SELECT
787	!
788
789	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol')
790
791	END SUBROUTINE dom_vvl_interpol
792
793	SUBROUTINE dom_vvl_rst( kt, cdrw )
794	!!---------------------------------------------------------------------
795	!! * ROUTINE dom_vvl_rst *
796	!!
797	!! ** Purpose : Read or write VVL file in restart file
798	!!
799	!! ** Method : use of IOM library
800	!! if the restart does not contain vertical scale factors,
801	!! they are set to the _0 values
802	!! if the restart does not contain vertical scale factors increments (z_tilde),
803	!! they are set to 0.
804	!!----------------------------------------------------------------------
805	!! * Arguments
806	INTEGER , INTENT(in) :: kt ! ocean time-step
807	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
808	!! * Local declarations
809	INTEGER :: jk
810	INTEGER :: id1, id2, id3, id4, id5 ! local integers
811	!!----------------------------------------------------------------------
812	!
813	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_rst')
814	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
815	! ! ===============
816	IF( ln_rstart ) THEN !* Read the restart file
817	CALL rst_read_open ! open the restart file if necessary
818	CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn )
819	!
820	id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
821	id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
822	id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
823	id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
824	id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
825	! ! --------- !
826	! ! all cases !
827	! ! --------- !
828	IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist
829	CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
830	CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
831	! needed to restart if land processor not computed
832	IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'
833	WHERE ( tmask(:,:,:) == 0.0_wp )
834	fse3t_n(:,:,:) = e3t_0(:,:,:)
835	fse3t_b(:,:,:) = e3t_0(:,:,:)
836	END WHERE
837	IF( neuler == 0 ) THEN
838	fse3t_b(:,:,:) = fse3t_n(:,:,:)
839	ENDIF
840	ELSE IF( id1 > 0 ) THEN
841	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
842	IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
843	IF(lwp) write(numout,*) 'neuler is forced to 0'
844	CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
845	fse3t_n(:,:,:) = fse3t_b(:,:,:)
846	neuler = 0
847	ELSE IF( id2 > 0 ) THEN
848	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
849	IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
850	IF(lwp) write(numout,*) 'neuler is forced to 0'
851	CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
852	fse3t_b(:,:,:) = fse3t_n(:,:,:)
853	neuler = 0
854	ELSE
855	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'
856	IF(lwp) write(numout,*) 'Compute scale factor from sshn'
857	IF(lwp) write(numout,*) 'neuler is forced to 0'
858	DO jk=1,jpk
859	fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
860	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
861	& + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
862	END DO
863	fse3t_b(:,:,:) = fse3t_n(:,:,:)
864	neuler = 0
865	ENDIF
866	! ! ----------- !
867	IF( ln_vvl_zstar ) THEN ! z_star case !
868	! ! ----------- !
869	IF( MIN( id3, id4 ) > 0 ) THEN
870	CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
871	ENDIF
872	! ! ----------------------- !
873	ELSE ! z_tilde and layer cases !
874	! ! ----------------------- !
875	IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist
876	CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
877	CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
878	ELSE ! one at least array is missing
879	tilde_e3t_b(:,:,:) = 0.0_wp
880	tilde_e3t_n(:,:,:) = 0.0_wp
881	ENDIF
882	! ! ------------ !
883	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
884	! ! ------------ !
885	IF( id5 > 0 ) THEN ! required array exists
886	CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
887	ELSE ! array is missing
888	hdiv_lf(:,:,:) = 0.0_wp
889	ENDIF
890	ENDIF
891	ENDIF
892	!
893	ELSE !* Initialize at "rest" or from analytical solution
894	sshn(:,:) = 0.0_wp
895	IF ( cp_cfg == 'vortex' ) THEN
896	CALL istate_vortex_ssh
897	ENDIF
898	DO jk=1,jpk
899	fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
900	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
901	& + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
902	END DO
903	fse3t_b(:,:,:) = fse3t_n(:,:,:)
904
905	IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
906	tilde_e3t_b(:,:,:) = 0.0_wp
907	tilde_e3t_n(:,:,:) = 0.0_wp
908	IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
909	END IF
910	ENDIF
911
912	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
913	! ! ===================
914	IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
915	! ! --------- !
916	! ! all cases !
917	! ! --------- !
918	CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
919	CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
920	! ! ----------------------- !
921	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
922	! ! ----------------------- !
923	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
924	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
925	END IF
926	! ! -------------!
927	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
928	! ! ------------ !
929	CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
930	ENDIF
931
932	ENDIF
933	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_rst')
934
935	END SUBROUTINE dom_vvl_rst
936
937
938	SUBROUTINE dom_vvl_ctl
939	!!---------------------------------------------------------------------
940	!! * ROUTINE dom_vvl_ctl *
941	!!
942	!! ** Purpose : Control the consistency between namelist options
943	!! for vertical coordinate
944	!!----------------------------------------------------------------------
945	INTEGER :: ioptio
946	INTEGER :: ios
947
948	NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
949	& ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , &
950	& rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg ! not yet implemented: ln_vvl_kepe
951	!!----------------------------------------------------------------------
952
953	REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist :
954	READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
955	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
956
957	REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run
958	READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
959	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
960	IF(lwm) WRITE ( numond, nam_vvl )
961
962	IF(lwp) THEN ! Namelist print
963	WRITE(numout,*)
964	WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
965	WRITE(numout,*) '~~~~~~~~~~~'
966	WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate'
967	WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar
968	WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde
969	WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer
970	WRITE(numout,*) ' ztilde as zstar ln_vvl_ztilde_as_zstar = ', ln_vvl_ztilde_as_zstar
971	WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor
972	! WRITE(numout,*) ' Namelist nam_vvl : chose kinetic-to-potential energy conservation'
973	! WRITE(numout,*) ' ln_vvl_kepe = ', ln_vvl_kepe
974	WRITE(numout,*) ' Namelist nam_vvl : thickness diffusion coefficient'
975	WRITE(numout,*) ' rn_ahe3 = ', rn_ahe3
976	WRITE(numout,*) ' Namelist nam_vvl : maximum e3t deformation fractional change'
977	WRITE(numout,*) ' rn_zdef_max = ', rn_zdef_max
978	IF( ln_vvl_ztilde_as_zstar ) THEN
979	WRITE(numout,*) ' ztilde running in zstar emulation mode; '
980	WRITE(numout,*) ' ignoring namelist timescale parameters and using:'
981	WRITE(numout,*) ' hard-wired : z-tilde to zstar restoration timescale (days)'
982	WRITE(numout,*) ' rn_rst_e3t = 0.0'
983	WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
984	WRITE(numout,*) ' rn_lf_cutoff = 1.0/rdt'
985	ELSE
986	WRITE(numout,*) ' Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
987	WRITE(numout,*) ' rn_rst_e3t = ', rn_rst_e3t
988	WRITE(numout,*) ' Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
989	WRITE(numout,*) ' rn_lf_cutoff = ', rn_lf_cutoff
990	ENDIF
991	WRITE(numout,*) ' Namelist nam_vvl : debug prints'
992	WRITE(numout,*) ' ln_vvl_dbg = ', ln_vvl_dbg
993	ENDIF
994
995	ioptio = 0 ! Parameter control
996	IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
997	IF( ln_vvl_zstar ) ioptio = ioptio + 1
998	IF( ln_vvl_ztilde ) ioptio = ioptio + 1
999	IF( ln_vvl_layer ) ioptio = ioptio + 1
1000
1001	IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1002	IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1003
1004	IF(lwp) THEN ! Print the choice
1005	WRITE(numout,*)
1006	IF( ln_vvl_zstar ) WRITE(numout,*) ' zstar vertical coordinate is used'
1007	IF( ln_vvl_ztilde ) WRITE(numout,*) ' ztilde vertical coordinate is used'
1008	IF( ln_vvl_layer ) WRITE(numout,*) ' layer vertical coordinate is used'
1009	IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) ' to emulate a zstar coordinate'
1010	! - ML - Option not developed yet
1011	! IF( ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option used'
1012	! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option not used'
1013	ENDIF
1014
1015	#if defined key_agrif
1016	IF ((.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar)) CALL ctl_stop( 'AGRIF implemented with zstar coordinate only (key_vvl)' )
1017	#endif
1018
1019	END SUBROUTINE dom_vvl_ctl
1020
1021	SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1022	!!---------------------------------------------------------------------
1023	!! * ROUTINE dom_vvl_orca_fix *
1024	!!
1025	!! ** Purpose : Correct surface weighted, horizontally interpolated,
1026	!! scale factors at locations that have been individually
1027	!! modified in domhgr. Such modifications break the
1028	!! relationship between e12t and e1u*e2u etc.
1029	!! Recompute some scale factors ignoring the modified metric.
1030	!!----------------------------------------------------------------------
1031	!! * Arguments
1032	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated
1033	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3
1034	CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors
1035	! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'
1036	!! * Local declarations
1037	INTEGER :: ji, jj, jk ! dummy loop indices
1038	INTEGER :: ij0, ij1, ii0, ii1 ! dummy loop indices
1039	INTEGER :: isrow ! index for ORCA1 starting row
1040	!! acc
1041	!! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
1042	!! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
1043	!!
1044	! ! =====================
1045	IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN ! ORCA R2 configuration
1046	! ! =====================
1047	!! acc
1048	IF( nn_cla == 0 ) THEN
1049	!
1050	ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u was modified)
1051	ij0 = 102 ; ij1 = 102
1052	DO jk = 1, jpkm1
1053	DO jj = mj0(ij0), mj1(ij1)
1054	DO ji = mi0(ii0), mi1(ii1)
1055	SELECT CASE ( pout )
1056	CASE( 'U' )
1057	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1058	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1059	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1060	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1061	CASE( 'F' )
1062	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1063	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1064	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1065	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1066	END SELECT
1067	END DO
1068	END DO
1069	END DO
1070	!
1071	ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u and e1v were modified)
1072	ij0 = 88 ; ij1 = 88
1073	DO jk = 1, jpkm1
1074	DO jj = mj0(ij0), mj1(ij1)
1075	DO ji = mi0(ii0), mi1(ii1)
1076	SELECT CASE ( pout )
1077	CASE( 'U' )
1078	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1079	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1080	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1081	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1082	CASE( 'V' )
1083	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1084	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1085	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1086	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1087	CASE( 'F' )
1088	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1089	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1090	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1091	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1092	END SELECT
1093	END DO
1094	END DO
1095	END DO
1096	ENDIF
1097
1098	ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u was modified)
1099	ij0 = 116 ; ij1 = 116
1100	DO jk = 1, jpkm1
1101	DO jj = mj0(ij0), mj1(ij1)
1102	DO ji = mi0(ii0), mi1(ii1)
1103	SELECT CASE ( pout )
1104	CASE( 'U' )
1105	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1106	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1107	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1108	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1109	CASE( 'F' )
1110	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1111	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1112	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1113	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1114	END SELECT
1115	END DO
1116	END DO
1117	END DO
1118	ENDIF
1119	!
1120	! ! =====================
1121	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
1122	! ! =====================
1123	! This dirty section will be suppressed by simplification process:
1124	! all this will come back in input files
1125	! Currently these hard-wired indices relate to configuration with
1126	! extend grid (jpjglo=332)
1127	! which had a grid-size of 362x292.
1128	isrow = 332 - jpjglo
1129	!
1130	ii0 = 282 ; ii1 = 283 ! Gibraltar Strait (e2u was modified)
1131	ij0 = 241 - isrow ; ij1 = 241 - isrow
1132	DO jk = 1, jpkm1
1133	DO jj = mj0(ij0), mj1(ij1)
1134	DO ji = mi0(ii0), mi1(ii1)
1135	SELECT CASE ( pout )
1136	CASE( 'U' )
1137	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1138	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1139	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1140	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1141	CASE( 'F' )
1142	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1143	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1144	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1145	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1146	END SELECT
1147	END DO
1148	END DO
1149	END DO
1150	!
1151	ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u was modified)
1152	ij0 = 248 - isrow ; ij1 = 248 - isrow
1153	DO jk = 1, jpkm1
1154	DO jj = mj0(ij0), mj1(ij1)
1155	DO ji = mi0(ii0), mi1(ii1)
1156	SELECT CASE ( pout )
1157	CASE( 'U' )
1158	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1159	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1160	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1161	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1162	CASE( 'F' )
1163	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1164	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1165	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1166	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1167	END SELECT
1168	END DO
1169	END DO
1170	END DO
1171	!
1172	ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v was modified)
1173	ij0 = 164 - isrow ; ij1 = 165 - isrow
1174	DO jk = 1, jpkm1
1175	DO jj = mj0(ij0), mj1(ij1)
1176	DO ji = mi0(ii0), mi1(ii1)
1177	SELECT CASE ( pout )
1178	CASE( 'V' )
1179	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1180	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1181	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1182	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1183	END SELECT
1184	END DO
1185	END DO
1186	END DO
1187	!
1188	ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
1189	ij0 = 164 - isrow ; ij1 = 165 - isrow
1190	DO jk = 1, jpkm1
1191	DO jj = mj0(ij0), mj1(ij1)
1192	DO ji = mi0(ii0), mi1(ii1)
1193	SELECT CASE ( pout )
1194	CASE( 'V' )
1195	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1196	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1197	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1198	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1199	END SELECT
1200	END DO
1201	END DO
1202	END DO
1203	!
1204	ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v was modified)
1205	ij0 = 164 - isrow ; ij1 = 165 - isrow
1206	DO jk = 1, jpkm1
1207	DO jj = mj0(ij0), mj1(ij1)
1208	DO ji = mi0(ii0), mi1(ii1)
1209	SELECT CASE ( pout )
1210	CASE( 'V' )
1211	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1212	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1213	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1214	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1215	END SELECT
1216	END DO
1217	END DO
1218	END DO
1219	!
1220	ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v was modified)
1221	ij0 = 164 - isrow ; ij1 = 165 - isrow
1222	DO jk = 1, jpkm1
1223	DO jj = mj0(ij0), mj1(ij1)
1224	DO ji = mi0(ii0), mi1(ii1)
1225	SELECT CASE ( pout )
1226	CASE( 'V' )
1227	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1228	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1229	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1230	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1231	END SELECT
1232	END DO
1233	END DO
1234	END DO
1235	!
1236	ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v was modified)
1237	ij0 = 181 - isrow ; ij1 = 182 - isrow
1238	DO jk = 1, jpkm1
1239	DO jj = mj0(ij0), mj1(ij1)
1240	DO ji = mi0(ii0), mi1(ii1)
1241	SELECT CASE ( pout )
1242	CASE( 'V' )
1243	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1244	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1245	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1246	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1247	END SELECT
1248	END DO
1249	END DO
1250	END DO
1251	!
1252	ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v was modified)
1253	ij0 = 181 - isrow ; ij1 = 182 - isrow
1254	DO jk = 1, jpkm1
1255	DO jj = mj0(ij0), mj1(ij1)
1256	DO ji = mi0(ii0), mi1(ii1)
1257	SELECT CASE ( pout )
1258	CASE( 'V' )
1259	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1260	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1261	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1262	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1263	END SELECT
1264	END DO
1265	END DO
1266	END DO
1267	ENDIF
1268	! ! =====================
1269	IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration
1270	! ! =====================
1271	!
1272	ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u was modified)
1273	ij0 = 327 ; ij1 = 327
1274	DO jk = 1, jpkm1
1275	DO jj = mj0(ij0), mj1(ij1)
1276	DO ji = mi0(ii0), mi1(ii1)
1277	SELECT CASE ( pout )
1278	CASE( 'U' )
1279	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1280	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1281	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1282	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1283	CASE( 'F' )
1284	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1285	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1286	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1287	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1288	END SELECT
1289	END DO
1290	END DO
1291	END DO
1292	!
1293	ii0 = 627 ; ii1 = 628 ! Bosphorus Strait (e2u was modified)
1294	ij0 = 343 ; ij1 = 343
1295	DO jk = 1, jpkm1
1296	DO jj = mj0(ij0), mj1(ij1)
1297	DO ji = mi0(ii0), mi1(ii1)
1298	SELECT CASE ( pout )
1299	CASE( 'U' )
1300	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1301	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1302	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1303	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1304	CASE( 'F' )
1305	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1306	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1307	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1308	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1309	END SELECT
1310	END DO
1311	END DO
1312	END DO
1313	!
1314	ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u was modified)
1315	ij0 = 232 ; ij1 = 232
1316	DO jk = 1, jpkm1
1317	DO jj = mj0(ij0), mj1(ij1)
1318	DO ji = mi0(ii0), mi1(ii1)
1319	SELECT CASE ( pout )
1320	CASE( 'U' )
1321	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1322	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1323	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1324	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1325	CASE( 'F' )
1326	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1327	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1328	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1329	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1330	END SELECT
1331	END DO
1332	END DO
1333	END DO
1334	!
1335	ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u was modified)
1336	ij0 = 232 ; ij1 = 232
1337	DO jk = 1, jpkm1
1338	DO jj = mj0(ij0), mj1(ij1)
1339	DO ji = mi0(ii0), mi1(ii1)
1340	SELECT CASE ( pout )
1341	CASE( 'U' )
1342	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1343	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1344	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1345	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1346	CASE( 'F' )
1347	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1348	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1349	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1350	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1351	END SELECT
1352	END DO
1353	END DO
1354	END DO
1355	!
1356	ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u was modified)
1357	ij0 = 270 ; ij1 = 270
1358	DO jk = 1, jpkm1
1359	DO jj = mj0(ij0), mj1(ij1)
1360	DO ji = mi0(ii0), mi1(ii1)
1361	SELECT CASE ( pout )
1362	CASE( 'U' )
1363	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
1364	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
1365	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1366	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
1367	CASE( 'F' )
1368	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
1369	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
1370	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1371	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
1372	END SELECT
1373	END DO
1374	END DO
1375	END DO
1376	!
1377	ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v was modified)
1378	ij0 = 232 ; ij1 = 233
1379	DO jk = 1, jpkm1
1380	DO jj = mj0(ij0), mj1(ij1)
1381	DO ji = mi0(ii0), mi1(ii1)
1382	SELECT CASE ( pout )
1383	CASE( 'V' )
1384	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1385	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1386	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1387	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1388	END SELECT
1389	END DO
1390	END DO
1391	END DO
1392	!
1393	ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v was modified)
1394	ij0 = 276 ; ij1 = 276
1395	DO jk = 1, jpkm1
1396	DO jj = mj0(ij0), mj1(ij1)
1397	DO ji = mi0(ii0), mi1(ii1)
1398	SELECT CASE ( pout )
1399	CASE( 'V' )
1400	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
1401	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
1402	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1403	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
1404	END SELECT
1405	END DO
1406	END DO
1407	END DO
1408	ENDIF
1409	END SUBROUTINE dom_vvl_orca_fix
1410
1411	SUBROUTINE istate_vortex_ssh
1412	!!----------------------------------------------------------------------
1413	!! * ROUTINE istate_vortex_ssh *
1414	!!
1415	!! ** Purpose : Initialization of baroclinic vortex ssh
1416	!!
1417	!! ** Method : - set ssh
1418	!!----------------------------------------------------------------------
1419	INTEGER :: ji, jj, jk ! dummy loop indices
1420	REAL(wp) :: zx, zy, zP0, zumax, zlambda, zf0, zH, zrho1, za
1421	!!----------------------------------------------------------------------
1422	!
1423	!
1424	zf0 = 2._wp * omega * SIN( rad * ppgphi0 )
1425	zumax = 1._wp
1426	zlambda = SQRT(2._wp)*60.e3
1427	zH = 0.5_wp * rn_bathy
1428	!
1429	zP0 = rau0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
1430	! Sea level:
1431	za = -zP0 * (1._wp-EXP(-zH)) / (grav*(zH-1._wp + exp(-zH)))
1432	DO ji=1, jpi
1433	DO jj=1, jpj
1434	zx = glamt(ji,jj) * 1.e3
1435	zy = gphit(ji,jj) * 1.e3
1436	zrho1 = rau0 + za * exp(-(zx2+zy2)/zlambda**2)
1437	sshn(ji,jj) = zP0 * exp(-(zx2+zy2)/zlambda*2)/(zrho1grav) * tmask(ji,jj,1)
1438	END DO
1439	END DO
1440	sshb(:,:) = sshn(:,:)
1441	!
1442	END SUBROUTINE istate_vortex_ssh
1443
1444	!!======================================================================
1445	END MODULE domvvl
1446
1447
1448

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

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

Context Navigation

source: branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/MY_SRC/domvvl.F90 @ 8011

Download in other formats: