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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

domvvl.F90 in branches/UKMO/dev_r5518_nemovar_init_vvl/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_nemovar_init_vvl/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 6666

Last change on this file since 6666 was 6666, checked in by blemieux, 8 years ago

T/S profile assimilation on the Shelf

Updating the previous commit

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

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

Download in other formats: