Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

domvvl.F90 @ 9176

Last change on this file since 9176 was 9176, checked in by andmirek, 6 years ago
#2001: OMP directives
File size: 77.1 KB

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

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/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 9176

Download in other formats: