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

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

domvvl.F90 in branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 5014

Last change on this file since 5014 was 5014, checked in by hliu, 9 years ago
upload the modifications for W/D based on r:4826
Property svn:keywords set to `Id`
File size: 73.1 KB

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

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

Download in other formats: