Context Navigation

domvvl.F90 @ 6491

Last change on this file since 6491 was 6491, checked in by davestorkey, 8 years ago

Commit remaining changes to UKMO/r5518_GO6_package branch from following branches:
UKMO/dev_r5021_nn_etau_revision@6238
UKMO/dev_r5107_mld_zint@5534
UKMO/dev_r5107_eorca025_closea@6390
UKMO/restart_datestamp@5539
UKMO/icebergs_latent_heat@5821
UKMO/icebergs_restart_single_file_corrected@6480
UKMO/product_diagnostics@5971
UKMO/antarctic_partial_slip@5961
Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r5961 cf. r5958 of /branches/UKMO/antarctic_partial_slip/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r6349 cf. r5962 of /branches/UKMO/product_diagnostics/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r6480 cf. r6479 of /branches/UKMO/icebergs_restart_single_file_corrected/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r5986 cf. r5852 of /branches/UKMO/icebergs_restart_single_file/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r5821 cf. r5808 of /branches/UKMO/icebergs_latent_heat/NEMOGCM@6490

File size: 75.5 KB

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 6491

Download in other formats: