Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

domvvl.F90 @ 11738

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

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

Download in other formats: