Context Navigation

domvvl.F90 @ 5727

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

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

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

Context Navigation

source: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 5727

Download in other formats: