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

Context Navigation

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 6348

Last change on this file since 6348 was 6348, checked in by cetlod, 8 years ago
bugfix: move the output of scale factor before time swapping and output some variables needs for offline, see ticket #1682
Property svn:keywords set to `Id`
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
[4292]	596	!
	597	! Time filter and swap of scale factors
	598	! =====================================
	599	! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
	600	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
	601	IF( neuler == 0 .AND. kt == nit000 ) THEN
	602	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
	603	ELSE
	604	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &
	605	& + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
	606	ENDIF
	607	tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
	608	ENDIF
[4488]	609	fsdept_b(:,:,:) = fsdept_n(:,:,:)
	610	fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
	611
[4292]	612	fse3t_n(:,:,:) = fse3t_a(:,:,:)
	613	fse3u_n(:,:,:) = fse3u_a(:,:,:)
	614	fse3v_n(:,:,:) = fse3v_a(:,:,:)
	615
	616	! Compute all missing vertical scale factor and depths
	617	! ====================================================
	618	! Horizontal scale factor interpolations
	619	! --------------------------------------
	620	! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
[4370]	621	! - JC - hu_b, hv_b, hur_b, hvr_b also
[4292]	622	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F' )
	623	! Vertical scale factor interpolations
	624	! ------------------------------------
	625	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )
	626	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
	627	CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
[4488]	628	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' )
[4292]	629	CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
	630	CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
	631	! t- and w- points depth
	632	! ----------------------
[5120]	633	! set the isf depth as it is in the initial step
[4292]	634	fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
	635	fsdepw_n(:,:,1) = 0.0_wp
	636	fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
[5120]	637
	638	DO jk = 2, jpk
	639	DO jj = 1,jpj
	640	DO ji = 1,jpi
	641	! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
	642	! 1 for jk = mikt
	643	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
[4990]	644	fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
[5120]	645	fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) &
	646	& + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))
	647	fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
[4990]	648	END DO
	649	END DO
[4292]	650	END DO
[5120]	651
[4292]	652	! Local depth and Inverse of the local depth of the water column at u- and v- points
	653	! ----------------------------------------------------------------------------------
[4370]	654	hu (:,:) = hu_a (:,:)
	655	hv (:,:) = hv_a (:,:)
	656
[4292]	657	! Inverse of the local depth
[4370]	658	hur(:,:) = hur_a(:,:)
	659	hvr(:,:) = hvr_a(:,:)
[4292]	660
[4370]	661	! Local depth of the water column at t- points
	662	! --------------------------------------------
	663	ht(:,:) = 0.
	664	DO jk = 1, jpkm1
	665	ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
	666	END DO
	667
[4292]	668	! write restart file
	669	! ==================
	670	IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
	671	!
	672	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp')
	673
	674	END SUBROUTINE dom_vvl_sf_swp
	675
	676
	677	SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
	678	!!---------------------------------------------------------------------
	679	!! * ROUTINE dom_vvl__interpol *
	680	!!
	681	!! ** Purpose : interpolate scale factors from one grid point to another
	682	!!
	683	!! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0)
	684	!! - horizontal interpolation: grid cell surface averaging
	685	!! - vertical interpolation: simple averaging
	686	!!----------------------------------------------------------------------
	687	!! * Arguments
	688	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated
	689	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3
	690	CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors
	691	! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'
	692	!! * Local declarations
	693	INTEGER :: ji, jj, jk ! dummy loop indices
	694	LOGICAL :: l_is_orca ! local logical
	695	!!----------------------------------------------------------------------
	696	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol')
	697	!
	698	l_is_orca = .FALSE.
	699	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE. ! ORCA R2 configuration - will need to correct some locations
	700
	701	SELECT CASE ( pout )
	702	! ! ------------------------------------- !
	703	CASE( 'U' ) ! interpolation from T-point to U-point !
	704	! ! ------------------------------------- !
	705	! horizontal surface weighted interpolation
	706	DO jk = 1, jpk
	707	DO jj = 1, jpjm1
	708	DO ji = 1, fs_jpim1 ! vector opt.
	709	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) &
	710	& * ( e12t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	711	& + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
	712	END DO
[2528]	713	END DO
	714	END DO
[4292]	715	!
	716	IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
	717	! boundary conditions
[4990]	718	CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
[4292]	719	pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
	720	! ! ------------------------------------- !
	721	CASE( 'V' ) ! interpolation from T-point to V-point !
	722	! ! ------------------------------------- !
	723	! horizontal surface weighted interpolation
	724	DO jk = 1, jpk
	725	DO jj = 1, jpjm1
	726	DO ji = 1, fs_jpim1 ! vector opt.
	727	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) &
	728	& * ( e12t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
	729	& + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
	730	END DO
	731	END DO
	732	END DO
	733	!
	734	IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
	735	! boundary conditions
[4990]	736	CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
[4292]	737	pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
	738	! ! ------------------------------------- !
	739	CASE( 'F' ) ! interpolation from U-point to F-point !
	740	! ! ------------------------------------- !
	741	! horizontal surface weighted interpolation
	742	DO jk = 1, jpk
	743	DO jj = 1, jpjm1
	744	DO ji = 1, fs_jpim1 ! vector opt.
	745	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) &
	746	& * ( e12u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) &
	747	& + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
	748	END DO
	749	END DO
	750	END DO
	751	!
	752	IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
	753	! boundary conditions
[4990]	754	CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
[4292]	755	pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
	756	! ! ------------------------------------- !
	757	CASE( 'W' ) ! interpolation from T-point to W-point !
	758	! ! ------------------------------------- !
	759	! vertical simple interpolation
	760	pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
	761	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
	762	DO jk = 2, jpk
	763	pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) &
	764	& + 0.5_wp * tmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) )
	765	END DO
	766	! ! -------------------------------------- !
	767	CASE( 'UW' ) ! interpolation from U-point to UW-point !
	768	! ! -------------------------------------- !
	769	! vertical simple interpolation
	770	pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
	771	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
	772	DO jk = 2, jpk
	773	pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) &
	774	& + 0.5_wp * umask(:,:,jk) * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) )
	775	END DO
	776	! ! -------------------------------------- !
	777	CASE( 'VW' ) ! interpolation from V-point to VW-point !
	778	! ! -------------------------------------- !
	779	! vertical simple interpolation
	780	pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
	781	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
	782	DO jk = 2, jpk
	783	pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) &
	784	& + 0.5_wp * vmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) )
	785	END DO
	786	END SELECT
	787	!
[3294]	788
[4292]	789	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol')
	790
	791	END SUBROUTINE dom_vvl_interpol
	792
	793	SUBROUTINE dom_vvl_rst( kt, cdrw )
	794	!!---------------------------------------------------------------------
	795	!! * ROUTINE dom_vvl_rst *
	796	!!
	797	!! ** Purpose : Read or write VVL file in restart file
	798	!!
	799	!! ** Method : use of IOM library
	800	!! if the restart does not contain vertical scale factors,
	801	!! they are set to the _0 values
	802	!! if the restart does not contain vertical scale factors increments (z_tilde),
	803	!! they are set to 0.
	804	!!----------------------------------------------------------------------
	805	!! * Arguments
	806	INTEGER , INTENT(in) :: kt ! ocean time-step
	807	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	808	!! * Local declarations
[4490]	809	INTEGER :: jk
[4292]	810	INTEGER :: id1, id2, id3, id4, id5 ! local integers
	811	!!----------------------------------------------------------------------
	812	!
	813	IF( nn_timing == 1 ) CALL timing_start('dom_vvl_rst')
	814	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
	815	! ! ===============
	816	IF( ln_rstart ) THEN !* Read the restart file
	817	CALL rst_read_open ! open the restart file if necessary
[4366]	818	CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn )
	819	!
[4292]	820	id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
	821	id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
	822	id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
	823	id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
[4795]	824	id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
[4292]	825	! ! --------- !
	826	! ! all cases !
	827	! ! --------- !
	828	IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist
	829	CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
	830	CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
[4990]	831	! needed to restart if land processor not computed
	832	IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'
	833	WHERE ( tmask(:,:,:) == 0.0_wp )
	834	fse3t_n(:,:,:) = e3t_0(:,:,:)
	835	fse3t_b(:,:,:) = e3t_0(:,:,:)
	836	END WHERE
[4292]	837	IF( neuler == 0 ) THEN
	838	fse3t_b(:,:,:) = fse3t_n(:,:,:)
	839	ENDIF
	840	ELSE IF( id1 > 0 ) THEN
[4990]	841	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
	842	IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
	843	IF(lwp) write(numout,*) 'neuler is forced to 0'
	844	CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
	845	fse3t_n(:,:,:) = fse3t_b(:,:,:)
	846	neuler = 0
	847	ELSE IF( id2 > 0 ) THEN
[4490]	848	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
	849	IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
	850	IF(lwp) write(numout,*) 'neuler is forced to 0'
[4990]	851	CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
[4292]	852	fse3t_b(:,:,:) = fse3t_n(:,:,:)
[4490]	853	neuler = 0
	854	ELSE
	855	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'
	856	IF(lwp) write(numout,*) 'Compute scale factor from sshn'
	857	IF(lwp) write(numout,*) 'neuler is forced to 0'
	858	DO jk=1,jpk
	859	fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
[4990]	860	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
[4490]	861	& + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
	862	END DO
	863	fse3t_b(:,:,:) = fse3t_n(:,:,:)
	864	neuler = 0
[4292]	865	ENDIF
	866	! ! ----------- !
	867	IF( ln_vvl_zstar ) THEN ! z_star case !
	868	! ! ----------- !
	869	IF( MIN( id3, id4 ) > 0 ) THEN
	870	CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
	871	ENDIF
	872	! ! ----------------------- !
	873	ELSE ! z_tilde and layer cases !
	874	! ! ----------------------- !
	875	IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist
	876	CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
	877	CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
	878	ELSE ! one at least array is missing
	879	tilde_e3t_b(:,:,:) = 0.0_wp
	880	tilde_e3t_n(:,:,:) = 0.0_wp
	881	ENDIF
	882	! ! ------------ !
	883	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
	884	! ! ------------ !
	885	IF( id5 > 0 ) THEN ! required array exists
	886	CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
	887	ELSE ! array is missing
	888	hdiv_lf(:,:,:) = 0.0_wp
	889	ENDIF
	890	ENDIF
	891	ENDIF
	892	!
	893	ELSE !* Initialize at "rest"
	894	fse3t_b(:,:,:) = e3t_0(:,:,:)
	895	fse3t_n(:,:,:) = e3t_0(:,:,:)
[4366]	896	sshn(:,:) = 0.0_wp
[4292]	897	IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
	898	tilde_e3t_b(:,:,:) = 0.0_wp
	899	tilde_e3t_n(:,:,:) = 0.0_wp
	900	IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
	901	END IF
	902	ENDIF
	903
	904	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
	905	! ! ===================
	906	IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
	907	! ! --------- !
	908	! ! all cases !
	909	! ! --------- !
	910	CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
	911	CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
	912	! ! ----------------------- !
	913	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
	914	! ! ----------------------- !
	915	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
	916	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
	917	END IF
	918	! ! -------------!
	919	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
	920	! ! ------------ !
	921	CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
	922	ENDIF
	923
	924	ENDIF
	925	IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_rst')
	926
	927	END SUBROUTINE dom_vvl_rst
	928
	929
	930	SUBROUTINE dom_vvl_ctl
	931	!!---------------------------------------------------------------------
	932	!! * ROUTINE dom_vvl_ctl *
	933	!!
	934	!! ** Purpose : Control the consistency between namelist options
	935	!! for vertical coordinate
	936	!!----------------------------------------------------------------------
	937	INTEGER :: ioptio
[4294]	938	INTEGER :: ios
[4292]	939
	940	NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
	941	& ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , &
	942	& rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg ! not yet implemented: ln_vvl_kepe
	943	!!----------------------------------------------------------------------
	944
[4294]	945	REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist :
	946	READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
	947	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
[4292]	948
[4294]	949	REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run
	950	READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
	951	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
[4624]	952	IF(lwm) WRITE ( numond, nam_vvl )
[4294]	953
[4292]	954	IF(lwp) THEN ! Namelist print
	955	WRITE(numout,*)
	956	WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
	957	WRITE(numout,*) '~~~~~~~~~~~'
	958	WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate'
	959	WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar
	960	WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde
	961	WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer
	962	WRITE(numout,*) ' ztilde as zstar ln_vvl_ztilde_as_zstar = ', ln_vvl_ztilde_as_zstar
	963	WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor
	964	! WRITE(numout,*) ' Namelist nam_vvl : chose kinetic-to-potential energy conservation'
	965	! WRITE(numout,*) ' ln_vvl_kepe = ', ln_vvl_kepe
	966	WRITE(numout,*) ' Namelist nam_vvl : thickness diffusion coefficient'
	967	WRITE(numout,*) ' rn_ahe3 = ', rn_ahe3
	968	WRITE(numout,*) ' Namelist nam_vvl : maximum e3t deformation fractional change'
	969	WRITE(numout,*) ' rn_zdef_max = ', rn_zdef_max
	970	IF( ln_vvl_ztilde_as_zstar ) THEN
	971	WRITE(numout,*) ' ztilde running in zstar emulation mode; '
	972	WRITE(numout,*) ' ignoring namelist timescale parameters and using:'
	973	WRITE(numout,*) ' hard-wired : z-tilde to zstar restoration timescale (days)'
	974	WRITE(numout,*) ' rn_rst_e3t = 0.0'
	975	WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
	976	WRITE(numout,*) ' rn_lf_cutoff = 1.0/rdt'
	977	ELSE
	978	WRITE(numout,*) ' Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
	979	WRITE(numout,*) ' rn_rst_e3t = ', rn_rst_e3t
	980	WRITE(numout,*) ' Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
	981	WRITE(numout,*) ' rn_lf_cutoff = ', rn_lf_cutoff
	982	ENDIF
	983	WRITE(numout,*) ' Namelist nam_vvl : debug prints'
	984	WRITE(numout,*) ' ln_vvl_dbg = ', ln_vvl_dbg
	985	ENDIF
	986
	987	ioptio = 0 ! Parameter control
	988	IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
	989	IF( ln_vvl_zstar ) ioptio = ioptio + 1
	990	IF( ln_vvl_ztilde ) ioptio = ioptio + 1
	991	IF( ln_vvl_layer ) ioptio = ioptio + 1
	992
	993	IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
[4990]	994	IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
[4292]	995
	996	IF(lwp) THEN ! Print the choice
	997	WRITE(numout,*)
	998	IF( ln_vvl_zstar ) WRITE(numout,*) ' zstar vertical coordinate is used'
	999	IF( ln_vvl_ztilde ) WRITE(numout,*) ' ztilde vertical coordinate is used'
	1000	IF( ln_vvl_layer ) WRITE(numout,*) ' layer vertical coordinate is used'
	1001	IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) ' to emulate a zstar coordinate'
	1002	! - ML - Option not developed yet
	1003	! IF( ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option used'
	1004	! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option not used'
	1005	ENDIF
	1006
[4486]	1007	#if defined key_agrif
	1008	IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )
	1009	#endif
	1010
[4292]	1011	END SUBROUTINE dom_vvl_ctl
	1012
	1013	SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
	1014	!!---------------------------------------------------------------------
	1015	!! * ROUTINE dom_vvl_orca_fix *
	1016	!!
	1017	!! ** Purpose : Correct surface weighted, horizontally interpolated,
	1018	!! scale factors at locations that have been individually
	1019	!! modified in domhgr. Such modifications break the
	1020	!! relationship between e12t and e1u*e2u etc.
	1021	!! Recompute some scale factors ignoring the modified metric.
	1022	!!----------------------------------------------------------------------
	1023	!! * Arguments
	1024	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated
	1025	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3
	1026	CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors
	1027	! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'
	1028	!! * Local declarations
	1029	INTEGER :: ji, jj, jk ! dummy loop indices
	1030	INTEGER :: ij0, ij1, ii0, ii1 ! dummy loop indices
[5385]	1031	INTEGER :: isrow ! index for ORCA1 starting row
[4292]	1032	!! acc
	1033	!! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
	1034	!! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
	1035	!!
[3294]	1036	! ! =====================
[4296]	1037	IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN ! ORCA R2 configuration
[3294]	1038	! ! =====================
[4292]	1039	!! acc
[3294]	1040	IF( nn_cla == 0 ) THEN
	1041	!
	1042	ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u was modified)
[4292]	1043	ij0 = 102 ; ij1 = 102
	1044	DO jk = 1, jpkm1
[3294]	1045	DO jj = mj0(ij0), mj1(ij1)
	1046	DO ji = mi0(ii0), mi1(ii1)
[4292]	1047	SELECT CASE ( pout )
	1048	CASE( 'U' )
	1049	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
	1050	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	1051	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
	1052	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
	1053	CASE( 'F' )
	1054	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
	1055	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
	1056	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
	1057	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
	1058	END SELECT
[3294]	1059	END DO
	1060	END DO
	1061	END DO
	1062	!
	1063	ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u and e1v were modified)
[4292]	1064	ij0 = 88 ; ij1 = 88
	1065	DO jk = 1, jpkm1
[3294]	1066	DO jj = mj0(ij0), mj1(ij1)
	1067	DO ji = mi0(ii0), mi1(ii1)
[4292]	1068	SELECT CASE ( pout )
	1069	CASE( 'U' )
	1070	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
	1071	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	1072	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
	1073	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
	1074	CASE( 'V' )
	1075	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
	1076	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
	1077	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
	1078	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
	1079	CASE( 'F' )
	1080	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
	1081	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
	1082	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
	1083	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
	1084	END SELECT
[3294]	1085	END DO
	1086	END DO
	1087	END DO
	1088	ENDIF
	1089
	1090	ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u was modified)
[4292]	1091	ij0 = 116 ; ij1 = 116
	1092	DO jk = 1, jpkm1
[3294]	1093	DO jj = mj0(ij0), mj1(ij1)
	1094	DO ji = mi0(ii0), mi1(ii1)
[4292]	1095	SELECT CASE ( pout )
	1096	CASE( 'U' )
	1097	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
	1098	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	1099	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
	1100	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
	1101	CASE( 'F' )
	1102	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
	1103	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
	1104	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
	1105	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
	1106	END SELECT
[3294]	1107	END DO
	1108	END DO
	1109	END DO
	1110	ENDIF
[4292]	1111	!
[3294]	1112	! ! =====================
	1113	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
	1114	! ! =====================
[5506]	1115	! This dirty section will be suppressed by simplification process:
	1116	! all this will come back in input files
	1117	! Currently these hard-wired indices relate to configuration with
	1118	! extend grid (jpjglo=332)
[5385]	1119	! which had a grid-size of 362x292.
[5506]	1120	isrow = 332 - jpjglo
[4292]	1121	!
[5385]	1122	ii0 = 282 ; ii1 = 283 ! Gibraltar Strait (e2u was modified)
[5506]	1123	ij0 = 241 - isrow ; ij1 = 241 - isrow
[4292]	1124	DO jk = 1, jpkm1
[3294]	1125	DO jj = mj0(ij0), mj1(ij1)
	1126	DO ji = mi0(ii0), mi1(ii1)
[4292]	1127	SELECT CASE ( pout )
	1128	CASE( 'U' )
	1129	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
	1130	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	1131	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
	1132	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
	1133	CASE( 'F' )
	1134	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
	1135	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
	1136	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
	1137	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
	1138	END SELECT
[3294]	1139	END DO
	1140	END DO
	1141	END DO
[4292]	1142	!
[5385]	1143	ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u was modified)
[5506]	1144	ij0 = 248 - isrow ; ij1 = 248 - isrow
[4292]	1145	DO jk = 1, jpkm1
[3294]	1146	DO jj = mj0(ij0), mj1(ij1)
	1147	DO ji = mi0(ii0), mi1(ii1)
[4292]	1148	SELECT CASE ( pout )
	1149	CASE( 'U' )
	1150	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
	1151	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	1152	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
	1153	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
	1154	CASE( 'F' )
	1155	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
	1156	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
	1157	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
	1158	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
	1159	END SELECT
[3294]	1160	END DO
	1161	END DO
	1162	END DO
[4292]	1163	!
[5385]	1164	ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v was modified)
[5506]	1165	ij0 = 164 - isrow ; ij1 = 165 - isrow
[4292]	1166	DO jk = 1, jpkm1
[3294]	1167	DO jj = mj0(ij0), mj1(ij1)
	1168	DO ji = mi0(ii0), mi1(ii1)
[4292]	1169	SELECT CASE ( pout )
	1170	CASE( 'V' )
	1171	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
	1172	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
	1173	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
	1174	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
	1175	END SELECT
[3294]	1176	END DO
	1177	END DO
	1178	END DO
[4292]	1179	!
[5385]	1180	ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
[5506]	1181	ij0 = 164 - isrow ; ij1 = 165 - isrow
[4292]	1182	DO jk = 1, jpkm1
[3294]	1183	DO jj = mj0(ij0), mj1(ij1)
	1184	DO ji = mi0(ii0), mi1(ii1)
[4292]	1185	SELECT CASE ( pout )
	1186	CASE( 'V' )
	1187	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
	1188	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
	1189	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
	1190	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
	1191	END SELECT
[3294]	1192	END DO
	1193	END DO
	1194	END DO
[4292]	1195	!
[5385]	1196	ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v was modified)
[5506]	1197	ij0 = 164 - isrow ; ij1 = 165 - isrow
[4292]	1198	DO jk = 1, jpkm1
[3294]	1199	DO jj = mj0(ij0), mj1(ij1)
	1200	DO ji = mi0(ii0), mi1(ii1)
[4292]	1201	SELECT CASE ( pout )
	1202	CASE( 'V' )
	1203	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
	1204	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
	1205	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
	1206	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
	1207	END SELECT
[3294]	1208	END DO
	1209	END DO
	1210	END DO
[4292]	1211	!
[5506]	1212	ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v was modified)
	1213	ij0 = 164 - isrow ; ij1 = 165 - isrow
[4292]	1214	DO jk = 1, jpkm1
[3294]	1215	DO jj = mj0(ij0), mj1(ij1)
	1216	DO ji = mi0(ii0), mi1(ii1)
[4292]	1217	SELECT CASE ( pout )
	1218	CASE( 'V' )
	1219	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
	1220	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
	1221	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
	1222	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
	1223	END SELECT
[3294]	1224	END DO
	1225	END DO
	1226	END DO
[4292]	1227	!
[5506]	1228	ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v was modified)
	1229	ij0 = 181 - isrow ; ij1 = 182 - isrow
[4292]	1230	DO jk = 1, jpkm1
[3294]	1231	DO jj = mj0(ij0), mj1(ij1)
	1232	DO ji = mi0(ii0), mi1(ii1)
[4292]	1233	SELECT CASE ( pout )
	1234	CASE( 'V' )
	1235	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
	1236	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
	1237	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
	1238	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
	1239	END SELECT
[3294]	1240	END DO
	1241	END DO
	1242	END DO
[4292]	1243	!
[5506]	1244	ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v was modified)
	1245	ij0 = 181 - isrow ; ij1 = 182 - isrow
[4292]	1246	DO jk = 1, jpkm1
[3294]	1247	DO jj = mj0(ij0), mj1(ij1)
	1248	DO ji = mi0(ii0), mi1(ii1)
[4292]	1249	SELECT CASE ( pout )
	1250	CASE( 'V' )
	1251	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
	1252	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
	1253	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
	1254	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
	1255	END SELECT
[3294]	1256	END DO
	1257	END DO
	1258	END DO
	1259	ENDIF
[4292]	1260	! ! =====================
[3294]	1261	IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration
[4292]	1262	! ! =====================
	1263	!
[3294]	1264	ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u was modified)
[4292]	1265	ij0 = 327 ; ij1 = 327
	1266	DO jk = 1, jpkm1
[3294]	1267	DO jj = mj0(ij0), mj1(ij1)
	1268	DO ji = mi0(ii0), mi1(ii1)
[4292]	1269	SELECT CASE ( pout )
	1270	CASE( 'U' )
	1271	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
	1272	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	1273	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
	1274	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
	1275	CASE( 'F' )
	1276	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
	1277	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
	1278	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
	1279	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
	1280	END SELECT
[3294]	1281	END DO
	1282	END DO
	1283	END DO
	1284	!
[4292]	1285	ii0 = 627 ; ii1 = 628 ! Bosphorus Strait (e2u was modified)
	1286	ij0 = 343 ; ij1 = 343
	1287	DO jk = 1, jpkm1
[3294]	1288	DO jj = mj0(ij0), mj1(ij1)
	1289	DO ji = mi0(ii0), mi1(ii1)
[4292]	1290	SELECT CASE ( pout )
	1291	CASE( 'U' )
	1292	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
	1293	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	1294	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
	1295	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
	1296	CASE( 'F' )
	1297	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
	1298	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
	1299	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
	1300	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
	1301	END SELECT
[3294]	1302	END DO
	1303	END DO
	1304	END DO
	1305	!
	1306	ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u was modified)
[4292]	1307	ij0 = 232 ; ij1 = 232
	1308	DO jk = 1, jpkm1
[3294]	1309	DO jj = mj0(ij0), mj1(ij1)
	1310	DO ji = mi0(ii0), mi1(ii1)
[4292]	1311	SELECT CASE ( pout )
	1312	CASE( 'U' )
	1313	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
	1314	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	1315	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
	1316	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
	1317	CASE( 'F' )
	1318	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
	1319	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
	1320	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
	1321	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
	1322	END SELECT
[3294]	1323	END DO
	1324	END DO
	1325	END DO
	1326	!
	1327	ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u was modified)
[4292]	1328	ij0 = 232 ; ij1 = 232
	1329	DO jk = 1, jpkm1
[3294]	1330	DO jj = mj0(ij0), mj1(ij1)
	1331	DO ji = mi0(ii0), mi1(ii1)
[4292]	1332	SELECT CASE ( pout )
	1333	CASE( 'U' )
	1334	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
	1335	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	1336	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
	1337	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
	1338	CASE( 'F' )
	1339	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
	1340	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
	1341	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
	1342	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
	1343	END SELECT
[3294]	1344	END DO
	1345	END DO
	1346	END DO
	1347	!
	1348	ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u was modified)
[4292]	1349	ij0 = 270 ; ij1 = 270
	1350	DO jk = 1, jpkm1
[3294]	1351	DO jj = mj0(ij0), mj1(ij1)
	1352	DO ji = mi0(ii0), mi1(ii1)
[4292]	1353	SELECT CASE ( pout )
	1354	CASE( 'U' )
	1355	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &
	1356	& * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
	1357	& + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
	1358	& ) / e1u(ji,jj) + e3u_0(ji,jj,jk)
	1359	CASE( 'F' )
	1360	pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &
	1361	& * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &
	1362	& + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
	1363	& ) / e1f(ji,jj) + e3f_0(ji,jj,jk)
	1364	END SELECT
[3294]	1365	END DO
	1366	END DO
	1367	END DO
	1368	!
	1369	ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v was modified)
[4292]	1370	ij0 = 232 ; ij1 = 233
	1371	DO jk = 1, jpkm1
[3294]	1372	DO jj = mj0(ij0), mj1(ij1)
	1373	DO ji = mi0(ii0), mi1(ii1)
[4292]	1374	SELECT CASE ( pout )
	1375	CASE( 'V' )
	1376	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
	1377	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
	1378	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
	1379	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
	1380	END SELECT
[3294]	1381	END DO
	1382	END DO
	1383	END DO
	1384	!
	1385	ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v was modified)
[4292]	1386	ij0 = 276 ; ij1 = 276
	1387	DO jk = 1, jpkm1
[3294]	1388	DO jj = mj0(ij0), mj1(ij1)
	1389	DO ji = mi0(ii0), mi1(ii1)
[4292]	1390	SELECT CASE ( pout )
	1391	CASE( 'V' )
	1392	pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &
	1393	& * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
	1394	& + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
	1395	& ) / e2v(ji,jj) + e3v_0(ji,jj,jk)
	1396	END SELECT
[3294]	1397	END DO
	1398	END DO
	1399	END DO
	1400	ENDIF
[4292]	1401	END SUBROUTINE dom_vvl_orca_fix
[3294]	1402
[592]	1403	!!======================================================================
	1404	END MODULE domvvl
[4370]	1405
	1406
	1407

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

Download in other formats: