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

Context Navigation

source: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 5208

Last change on this file since 5208 was 5208, checked in by davestorkey, 9 years ago
Merge in changes from trunk up to 5021.
Property svn:keywords set to `Id`
File size: 75.0 KB

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

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

Download in other formats: