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

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

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

Context Navigation

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

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

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

Download in other formats: