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.

dynspg_ts.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9256

Last change on this file since 9256 was 9169, checked in by gm, 6 years ago
dev_merge_2017: all SRC: finalize the removal of useless warning when reading namelist_cfg + remove all nn_closea + nn_msh replaced by a logical
Property svn:keywords set to `Id`
File size: 70.9 KB

Rev	Line
[358]	1	MODULE dynspg_ts
[9023]	2
	3	!! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out !
	4
[358]	5	!!======================================================================
[6140]	6	!! * MODULE dynspg_ts *
	7	!! Ocean dynamics: surface pressure gradient trend, split-explicit scheme
	8	!!======================================================================
[1502]	9	!! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code
	10	!! - ! 2005-11 (V. Garnier, G. Madec) optimization
	11	!! - ! 2006-08 (S. Masson) distributed restart using iom
	12	!! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines
	13	!! - ! 2008-01 (R. Benshila) change averaging method
	14	!! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
[2528]	15	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
[2724]	16	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
[4292]	17	!! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping
	18	!! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility
[5930]	19	!! 3.7 ! 2015-11 (J. Chanut) free surface simplification
[7646]	20	!! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence
[9019]	21	!! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90)
[2724]	22	!!---------------------------------------------------------------------
[6140]	23
[358]	24	!!----------------------------------------------------------------------
[6140]	25	!! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme
	26	!! dyn_spg_ts_init: initialisation of the time-splitting scheme
	27	!! ts_wgt : set time-splitting weights for temporal averaging (or not)
	28	!! ts_rst : read/write time-splitting fields in restart file
[358]	29	!!----------------------------------------------------------------------
	30	USE oce ! ocean dynamics and tracers
	31	USE dom_oce ! ocean space and time domain
[888]	32	USE sbc_oce ! surface boundary condition: ocean
[9019]	33	USE zdf_oce ! vertical physics: variables
	34	USE zdfdrg ! vertical physics: top/bottom drag coef.
[5120]	35	USE sbcisf ! ice shelf variable (fwfisf)
[6140]	36	USE sbcapr ! surface boundary condition: atmospheric pressure
	37	USE dynadv , ONLY: ln_dynadv_vec
[358]	38	USE phycst ! physical constants
	39	USE dynvor ! vorticity term
[6152]	40	USE wet_dry ! wetting/drying flux limter
[7646]	41	USE bdy_oce ! open boundary
[5930]	42	USE bdytides ! open boundary condition data
[3294]	43	USE bdydyn2d ! open boundary conditions on barotropic variables
[4292]	44	USE sbctide ! tides
	45	USE updtide ! tide potential
[7646]	46	USE sbcwave ! surface wave
[9019]	47	USE diatmb ! Top,middle,bottom output
	48	#if defined key_agrif
	49	USE agrif_opa_interp ! agrif
[9124]	50	USE agrif_oce
[9019]	51	#endif
	52	#if defined key_asminc
	53	USE asminc ! Assimilation increment
	54	#endif
[6140]	55	!
	56	USE in_out_manager ! I/O manager
[358]	57	USE lib_mpp ! distributed memory computing library
	58	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	59	USE prtctl ! Print control
[2715]	60	USE iom ! IOM library
[4292]	61	USE restart ! only for lrst_oce
[9023]	62	USE diatmb ! Top,middle,bottom output
[358]	63
	64	IMPLICIT NONE
	65	PRIVATE
	66
[9124]	67	PUBLIC dyn_spg_ts ! called by dyn_spg
	68	PUBLIC dyn_spg_ts_init ! - - dyn_spg_init
[358]	69
[9019]	70	!! Time filtered arrays at baroclinic time step:
	71	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step
[9124]	72	!
[9023]	73	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
	74	REAL(wp),SAVE :: rdtbt ! Barotropic time step
[9019]	75	!
	76	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields
[9124]	77	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points
	78	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
	79	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
[4292]	80
[9043]	81	REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios
	82	REAL(wp) :: r1_8 = 0.125_wp !
	83	REAL(wp) :: r1_4 = 0.25_wp !
	84	REAL(wp) :: r1_2 = 0.5_wp !
[508]	85
[358]	86	!! * Substitutions
	87	# include "vectopt_loop_substitute.h90"
[2715]	88	!!----------------------------------------------------------------------
[9019]	89	!! NEMO/OPA 4.0 , NEMO Consortium (2017)
[5217]	90	!! $Id$
[2715]	91	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	92	!!----------------------------------------------------------------------
[358]	93	CONTAINS
	94
[2715]	95	INTEGER FUNCTION dyn_spg_ts_alloc()
	96	!!----------------------------------------------------------------------
	97	!! * routine dyn_spg_ts_alloc *
	98	!!----------------------------------------------------------------------
[6140]	99	INTEGER :: ierr(3)
[4292]	100	!!----------------------------------------------------------------------
	101	ierr(:) = 0
[6140]	102	!
	103	ALLOCATE( wgtbtp1(3nn_baro), wgtbtp2(3nn_baro), zwz(jpi,jpj), STAT=ierr(1) )
	104	!
	105	IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &
	106	& ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) )
	107	!
	108	ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) )
	109	!
	110	dyn_spg_ts_alloc = MAXVAL( ierr(:) )
	111	!
[2715]	112	IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc )
[5930]	113	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays')
[2715]	114	!
	115	END FUNCTION dyn_spg_ts_alloc
	116
[5836]	117
[358]	118	SUBROUTINE dyn_spg_ts( kt )
	119	!!----------------------------------------------------------------------
	120	!!
[6140]	121	!! ** Purpose : - Compute the now trend due to the explicit time stepping
	122	!! of the quasi-linear barotropic system, and add it to the
	123	!! general momentum trend.
[358]	124	!!
[6140]	125	!! ** Method : - split-explicit schem (time splitting) :
[4374]	126	!! Barotropic variables are advanced from internal time steps
	127	!! "n" to "n+1" if ln_bt_fw=T
	128	!! or from
	129	!! "n-1" to "n+1" if ln_bt_fw=F
	130	!! thanks to a generalized forward-backward time stepping (see ref. below).
[358]	131	!!
[4374]	132	!! ** Action :
	133	!! -Update the filtered free surface at step "n+1" : ssha
	134	!! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
[9023]	135	!! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv
[4374]	136	!! These are used to advect tracers and are compliant with discrete
	137	!! continuity equation taken at the baroclinic time steps. This
	138	!! ensures tracers conservation.
[6140]	139	!! - (ua, va) momentum trend updated with barotropic component.
[358]	140	!!
[6140]	141	!! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
[358]	142	!!---------------------------------------------------------------------
[1502]	143	INTEGER, INTENT(in) :: kt ! ocean time-step index
[2715]	144	!
[4292]	145	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[9019]	146	LOGICAL :: ll_fw_start ! =T : forward integration
	147	LOGICAL :: ll_init ! =T : special startup of 2d equations
	148	LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D
	149	INTEGER :: ikbu, iktu, noffset ! local integers
	150	INTEGER :: ikbv, iktv ! - -
[9043]	151	REAL(wp) :: r1_2dt_b, z2dt_bf ! local scalars
[9019]	152	REAL(wp) :: zx1, zx2, zu_spg, zhura ! - -
	153	REAL(wp) :: zy1, zy2, zv_spg, zhvra ! - -
	154	REAL(wp) :: za0, za1, za2, za3 ! - -
	155	REAL(wp) :: zmdi, zztmp ! - -
	156	REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf
	157	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc
	158	REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv
	159	REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e
	160	REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e
	161	REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points
[3294]	162	!
[9023]	163	REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True.
	164
	165	INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point
	166
	167	REAL(wp) :: zepsilon, zgamma ! - -
[9019]	168	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef.
[9023]	169	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
	170	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask
[358]	171	!!----------------------------------------------------------------------
[3294]	172	!
[9023]	173	IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) )
	174	! !* Allocate temporary arrays
	175	IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
[3294]	176	!
[6140]	177	zmdi=1.e+20 ! missing data indicator for masking
[9019]	178	!
[9023]	179	zwdramp = r_rn_wdmin1 ! simplest ramp
	180	! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
	181	! ! reciprocal of baroclinic time step
[6140]	182	IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt
	183	ELSE ; z2dt_bf = 2.0_wp * rdt
[4292]	184	ENDIF
[9043]	185	r1_2dt_b = 1.0_wp / z2dt_bf
[4292]	186	!
[9023]	187	ll_init = ln_bt_av ! if no time averaging, then no specific restart
[4292]	188	ll_fw_start = .FALSE.
[9023]	189	! ! time offset in steps for bdy data update
[6140]	190	IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro
	191	ELSE ; noffset = 0
	192	ENDIF
[4292]	193	!
[9023]	194	IF( kt == nit000 ) THEN !* initialisation
[508]	195	!
[358]	196	IF(lwp) WRITE(numout,*)
	197	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
	198	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
[4354]	199	IF(lwp) WRITE(numout,*)
[1502]	200	!
[6140]	201	IF( neuler == 0 ) ll_init=.TRUE.
[1502]	202	!
[6140]	203	IF( ln_bt_fw .OR. neuler == 0 ) THEN
	204	ll_fw_start =.TRUE.
	205	noffset = 0
[4292]	206	ELSE
[6140]	207	ll_fw_start =.FALSE.
[4292]	208	ENDIF
	209	!
	210	! Set averaging weights and cycle length:
[6140]	211	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
[4292]	212	!
	213	ENDIF
	214	!
[9019]	215	IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)
	216	DO jj = 2, jpjm1
	217	DO ji = fs_2, fs_jpim1 ! vector opt.
[9112]	218	zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
	219	zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
[9019]	220	END DO
	221	END DO
	222	ELSE ! bottom friction only
	223	DO jj = 2, jpjm1
	224	DO ji = fs_2, fs_jpim1 ! vector opt.
[9112]	225	zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
	226	zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
[9019]	227	END DO
	228	END DO
	229	ENDIF
	230	!
[4292]	231	! Set arrays to remove/compute coriolis trend.
	232	! Do it once at kt=nit000 if volume is fixed, else at each long time step.
	233	! Note that these arrays are also used during barotropic loop. These are however frozen
[4374]	234	! although they should be updated in the variable volume case. Not a big approximation.
[4292]	235	! To remove this approximation, copy lines below inside barotropic loop
[4374]	236	! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
[4292]	237	!
[6140]	238	IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN
	239	IF( ln_dynvor_een ) THEN !== EEN scheme ==!
[7646]	240	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
[5836]	241	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
	242	DO jj = 1, jpjm1
	243	DO ji = 1, jpim1
[6140]	244	zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &
	245	& ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp
[7646]	246	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
[5836]	247	END DO
[5032]	248	END DO
[5836]	249	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
	250	DO jj = 1, jpjm1
	251	DO ji = 1, jpim1
[6140]	252	zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &
	253	& ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) &
[5836]	254	& / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + &
[4292]	255	& tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) )
[7646]	256	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
[5836]	257	END DO
[4292]	258	END DO
[5836]	259	END SELECT
[4292]	260	CALL lbc_lnk( zwz, 'F', 1._wp )
[5836]	261	!
[7753]	262	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[358]	263	DO jj = 2, jpj
[5836]	264	DO ji = 2, jpi
[4292]	265	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	266	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	267	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	268	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
[358]	269	END DO
	270	END DO
[5836]	271	!
	272	ELSE !== all other schemes (ENE, ENS, MIX)
[7753]	273	zwz(:,:) = 0._wp
	274	zhf(:,:) = 0._wp
[7646]	275
[9124]	276	!!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
[7646]	277	!!gm A priori a better value should be something like :
	278	!!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
	279	!!gm divided by the sum of the corresponding mask
	280	!!gm
	281	!!
[9019]	282	IF( .NOT.ln_sco ) THEN
[7646]	283
	284	!!gm agree the JC comment : this should be done in a much clear way
	285
	286	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
	287	! Set it to zero for the time being
	288	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
	289	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
	290	! ENDIF
	291	! zhf(:,:) = gdepw_0(:,:,jk+1)
	292	ELSE
	293	!zhf(:,:) = hbatf(:,:)
	294	DO jj = 1, jpjm1
	295	DO ji = 1, jpim1
	296	zhf(ji,jj) = MAX( 0._wp, &
	297	& ( ht_0(ji ,jj )*tmask(ji ,jj ,1) + &
	298	& ht_0(ji+1,jj )*tmask(ji+1,jj ,1) + &
	299	& ht_0(ji ,jj+1)*tmask(ji ,jj+1,1) + &
	300	& ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) / &
	301	& ( tmask(ji ,jj ,1) + tmask(ji+1,jj ,1) +&
	302	& tmask(ji ,jj+1,1) + tmask(ji+1,jj+1,1) +&
	303	& rsmall ) )
	304	END DO
	305	END DO
	306	END IF
	307
	308	DO jj = 1, jpjm1
[7753]	309	zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
[7646]	310	END DO
	311	!!gm end
[5836]	312
[4292]	313	DO jk = 1, jpkm1
	314	DO jj = 1, jpjm1
[7753]	315	zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
[4292]	316	END DO
	317	END DO
[4370]	318	CALL lbc_lnk( zhf, 'F', 1._wp )
[4292]	319	! JC: TBC. hf should be greater than 0
	320	DO jj = 1, jpj
	321	DO ji = 1, jpi
[4370]	322	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
[4292]	323	END DO
	324	END DO
[7753]	325	zwz(:,:) = ff_f(:,:) * zwz(:,:)
[358]	326	ENDIF
[508]	327	ENDIF
[1502]	328	!
[4292]	329	! If forward start at previous time step, and centered integration,
	330	! then update averaging weights:
[5836]	331	IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
[4292]	332	ll_fw_start=.FALSE.
[9019]	333	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
[4292]	334	ENDIF
	335
[358]	336	! -----------------------------------------------------------------------------
	337	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
	338	! -----------------------------------------------------------------------------
[1502]	339	!
[4292]	340	!
[4354]	341	! !* e3*d/dt(Ua) (Vertically integrated)
[4292]	342	! ! --------------------------------------------------
[7753]	343	zu_frc(:,:) = 0._wp
	344	zv_frc(:,:) = 0._wp
[1502]	345	!
	346	DO jk = 1, jpkm1
[7753]	347	zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
	348	zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
[1502]	349	END DO
[4292]	350	!
[7753]	351	zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:)
	352	zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:)
[4292]	353	!
[7753]	354	!
[1502]	355	! !* baroclinic momentum trend (remove the vertical mean trend)
[4292]	356	DO jk = 1, jpkm1 ! -----------------------------------------------------------
[1502]	357	DO jj = 2, jpjm1
	358	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	359	ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
	360	va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
[1502]	361	END DO
[358]	362	END DO
[1502]	363	END DO
[7646]	364
	365	!!gm Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
	366	!!gm Is it correct to do so ? I think so...
	367
	368
[4292]	369	! !* barotropic Coriolis trends (vorticity scheme dependent)
	370	! ! --------------------------------------------------------
[7753]	371	zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes
	372	zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)
[1502]	373	!
[358]	374	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
	375	DO jj = 2, jpjm1
	376	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	377	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
	378	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
	379	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
	380	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
[358]	381	! energy conserving formulation for planetary vorticity term
[9043]	382	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	383	zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
[358]	384	END DO
	385	END DO
[508]	386	!
[4374]	387	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
[358]	388	DO jj = 2, jpjm1
	389	DO ji = fs_2, fs_jpim1 ! vector opt.
[9043]	390	zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
[5836]	391	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
[9043]	392	zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
[5836]	393	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
[4292]	394	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	395	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
[358]	396	END DO
	397	END DO
[508]	398	!
[5836]	399	ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme
[358]	400	DO jj = 2, jpjm1
	401	DO ji = fs_2, fs_jpim1 ! vector opt.
[9043]	402	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
[5836]	403	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
	404	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
	405	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
[9043]	406	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
[5836]	407	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
	408	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
	409	& + ftne(ji,jj ) * zwx(ji ,jj ) )
[358]	410	END DO
	411	END DO
[508]	412	!
[4292]	413	ENDIF
	414	!
[1502]	415	! !* Right-Hand-Side of the barotropic momentum equation
	416	! ! ----------------------------------------------------
[6140]	417	IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient
[9023]	418	IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters
	419	DO jj = 2, jpjm1
	420	DO ji = 2, jpim1
	421	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > &
	422	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
	423	& MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) &
	424	& > rn_wdmin1 + rn_wdmin2
	425	ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( &
	426	& MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > &
	427	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
	428
	429	IF(ll_tmp1) THEN
	430	zcpx(ji,jj) = 1.0_wp
	431	ELSE IF(ll_tmp2) THEN
	432	! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here
	433	zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
	434	& / (sshn(ji+1,jj) - sshn(ji ,jj)) )
	435	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
	436
	437	ELSE
	438	zcpx(ji,jj) = 0._wp
	439	END IF
	440
	441	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > &
	442	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
	443	& MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) &
	444	& > rn_wdmin1 + rn_wdmin2
	445	ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( &
	446	& MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > &
	447	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
	448
	449	IF(ll_tmp1) THEN
	450	zcpy(ji,jj) = 1.0_wp
	451	ELSE IF(ll_tmp2) THEN
	452	! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here
	453	zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
	454	& / (sshn(ji,jj+1) - sshn(ji,jj )) )
	455	zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp)
	456
	457	ELSE
	458	zcpy(ji,jj) = 0._wp
	459	END IF
	460	END DO
	461	END DO
	462
	463	DO jj = 2, jpjm1
	464	DO ji = 2, jpim1
	465	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) &
	466	& * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth
	467	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) &
	468	& * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth
	469
	470	END DO
	471	END DO
	472
[6152]	473	ELSE
[9019]	474	!
	475	DO jj = 2, jpjm1
	476	DO ji = fs_2, fs_jpim1 ! vector opt.
	477	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj)
	478	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj)
	479	END DO
	480	END DO
	481	ENDIF
	482	!
[1502]	483	ENDIF
[9019]	484	!
[4292]	485	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
[358]	486	DO ji = fs_2, fs_jpim1
[6140]	487	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
	488	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
[3294]	489	END DO
[4292]	490	END DO
	491	!
[9023]	492	! ! Add bottom stress contribution from baroclinic velocities:
	493	IF (ln_bt_fw) THEN
[4292]	494	DO jj = 2, jpjm1
	495	DO ji = fs_2, fs_jpim1 ! vector opt.
	496	ikbu = mbku(ji,jj)
	497	ikbv = mbkv(ji,jj)
	498	zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
	499	zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
	500	END DO
	501	END DO
[3294]	502	ELSE
[4292]	503	DO jj = 2, jpjm1
	504	DO ji = fs_2, fs_jpim1 ! vector opt.
	505	ikbu = mbku(ji,jj)
	506	ikbv = mbkv(ji,jj)
	507	zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
	508	zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
	509	END DO
	510	END DO
	511	ENDIF
[1502]	512	!
[4292]	513	! Note that the "unclipped" bottom friction parameter is used even with explicit drag
[9023]	514	IF( ln_wd_il ) THEN
[9075]	515	zztmp = -1._wp / rdtbt
[9045]	516	DO jj = 2, jpjm1
	517	DO ji = fs_2, fs_jpim1 ! vector opt.
[9112]	518	zu_frc(ji,jj) = zu_frc(ji,jj) + &
	519	& MAX(r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ), zztmp ) * zwx(ji,jj) * wdrampu(ji,jj)
	520	zv_frc(ji,jj) = zv_frc(ji,jj) + &
	521	& MAX(r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ), zztmp ) * zwy(ji,jj) * wdrampv(ji,jj)
[9045]	522	END DO
	523	END DO
[9023]	524	ELSE
[9045]	525	DO jj = 2, jpjm1
	526	DO ji = fs_2, fs_jpim1 ! vector opt.
[9112]	527	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj)
	528	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj)
[9045]	529	END DO
	530	END DO
[9023]	531	END IF
	532	!
[9019]	533	IF( ln_isfcav ) THEN ! Add TOP stress contribution from baroclinic velocities:
	534	IF( ln_bt_fw ) THEN
	535	DO jj = 2, jpjm1
	536	DO ji = fs_2, fs_jpim1 ! vector opt.
	537	iktu = miku(ji,jj)
	538	iktv = mikv(ji,jj)
	539	zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
	540	zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
	541	END DO
	542	END DO
	543	ELSE
	544	DO jj = 2, jpjm1
	545	DO ji = fs_2, fs_jpim1 ! vector opt.
	546	iktu = miku(ji,jj)
	547	iktv = mikv(ji,jj)
	548	zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
	549	zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
	550	END DO
	551	END DO
	552	ENDIF
	553	!
	554	! Note that the "unclipped" top friction parameter is used even with explicit drag
	555	DO jj = 2, jpjm1
	556	DO ji = fs_2, fs_jpim1 ! vector opt.
[9112]	557	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj)
	558	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj)
[9019]	559	END DO
	560	END DO
	561	ENDIF
[6140]	562	!
[9019]	563	IF( ln_bt_fw ) THEN ! Add wind forcing
	564	DO jj = 2, jpjm1
	565	DO ji = fs_2, fs_jpim1 ! vector opt.
	566	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)
	567	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)
	568	END DO
	569	END DO
[2724]	570	ELSE
[9043]	571	zztmp = r1_rau0 * r1_2
[9019]	572	DO jj = 2, jpjm1
	573	DO ji = fs_2, fs_jpim1 ! vector opt.
	574	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)
	575	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
	576	END DO
	577	END DO
[4292]	578	ENDIF
	579	!
[9019]	580	IF( ln_apr_dyn ) THEN ! Add atm pressure forcing
	581	IF( ln_bt_fw ) THEN
[4292]	582	DO jj = 2, jpjm1
	583	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	584	zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
	585	zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
[4292]	586	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
	587	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
	588	END DO
	589	END DO
	590	ELSE
[9043]	591	zztmp = grav * r1_2
[4292]	592	DO jj = 2, jpjm1
	593	DO ji = fs_2, fs_jpim1 ! vector opt.
[9019]	594	zu_spg = zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
	595	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
	596	zv_spg = zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
	597	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
[4292]	598	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
	599	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
	600	END DO
	601	END DO
	602	ENDIF
[2724]	603	ENDIF
[4292]	604	! !* Right-Hand-Side of the barotropic ssh equation
	605	! ! -----------------------------------------------
	606	! ! Surface net water flux and rivers
	607	IF (ln_bt_fw) THEN
[9019]	608	zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
[4292]	609	ELSE
[9043]	610	zztmp = r1_rau0 * r1_2
[9019]	611	zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) &
	612	& + fwfisf(:,:) + fwfisf_b(:,:) )
[4292]	613	ENDIF
[7646]	614	!
	615	IF( ln_sdw ) THEN ! Stokes drift divergence added if necessary
[7753]	616	zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
[7646]	617	ENDIF
	618	!
[4292]	619	#if defined key_asminc
	620	! ! Include the IAU weighted SSH increment
	621	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
[7753]	622	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
[4292]	623	ENDIF
	624	#endif
[5656]	625	! !* Fill boundary data arrays for AGRIF
	626	! ! ------------------------------------
[4486]	627	#if defined key_agrif
	628	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
	629	#endif
[4292]	630	!
[358]	631	! -----------------------------------------------------------------------
[4292]	632	! Phase 2 : Integration of the barotropic equations
[358]	633	! -----------------------------------------------------------------------
[1502]	634	!
	635	! ! ==================== !
	636	! ! Initialisations !
[4292]	637	! ! ==================== !
[4370]	638	! Initialize barotropic variables:
[4770]	639	IF( ll_init )THEN
[7753]	640	sshbb_e(:,:) = 0._wp
	641	ubb_e (:,:) = 0._wp
	642	vbb_e (:,:) = 0._wp
	643	sshb_e (:,:) = 0._wp
	644	ub_e (:,:) = 0._wp
	645	vb_e (:,:) = 0._wp
[4700]	646	ENDIF
[6152]	647
[4700]	648	!
[4370]	649	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
[7753]	650	sshn_e(:,:) = sshn(:,:)
	651	un_e (:,:) = un_b(:,:)
	652	vn_e (:,:) = vn_b(:,:)
	653	!
	654	hu_e (:,:) = hu_n(:,:)
	655	hv_e (:,:) = hv_n(:,:)
	656	hur_e (:,:) = r1_hu_n(:,:)
	657	hvr_e (:,:) = r1_hv_n(:,:)
[4370]	658	ELSE ! CENTRED integration: start from BEFORE fields
[7753]	659	sshn_e(:,:) = sshb(:,:)
	660	un_e (:,:) = ub_b(:,:)
	661	vn_e (:,:) = vb_b(:,:)
	662	!
	663	hu_e (:,:) = hu_b(:,:)
	664	hv_e (:,:) = hv_b(:,:)
	665	hur_e (:,:) = r1_hu_b(:,:)
	666	hvr_e (:,:) = r1_hv_b(:,:)
[4292]	667	ENDIF
	668	!
	669	!
[4370]	670	!
[4292]	671	! Initialize sums:
[7753]	672	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
	673	va_b (:,:) = 0._wp
	674	ssha (:,:) = 0._wp ! Sum for after averaged sea level
	675	un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop
	676	vn_adv(:,:) = 0._wp
[1502]	677	! ! ==================== !
[9023]	678
	679	IF (ln_wd_dl) THEN
	680	zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary)
	681	zvwdmask(:,:) = 0._wp !
	682	zuwdav2(:,:) = 0._wp
	683	zvwdav2(:,:) = 0._wp
	684	END IF
	685
	686
[4292]	687	DO jn = 1, icycle ! sub-time-step loop !
[1502]	688	! ! ==================== !
[3294]	689	! !* Update the forcing (BDY and tides)
[1502]	690	! ! ------------------
[4292]	691	! Update only tidal forcing at open boundaries
[7646]	692	IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
	693	IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset )
[4292]	694	!
	695	! Set extrapolation coefficients for predictor step:
	696	IF ((jn<3).AND.ll_init) THEN ! Forward
	697	za1 = 1._wp
	698	za2 = 0._wp
	699	za3 = 0._wp
	700	ELSE ! AB3-AM4 Coefficients: bet=0.281105
	701	za1 = 1.781105_wp ! za1 = 3/2 + bet
	702	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
	703	za3 = 0.281105_wp ! za3 = bet
	704	ENDIF
[367]	705
[4292]	706	! Extrapolate barotropic velocities at step jit+0.5:
[7753]	707	ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
	708	va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
[4292]	709
[6140]	710	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
[4292]	711	! ! ------------------
	712	! Extrapolate Sea Level at step jit+0.5:
[7753]	713	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
[9023]	714
	715	! set wetting & drying mask at tracer points for this barotropic sub-step
	716	IF ( ln_wd_dl ) THEN
	717
	718	IF ( ln_wd_dl_rmp ) THEN
	719	DO jj = 1, jpj
	720	DO ji = 1, jpi ! vector opt.
	721	IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
	722	! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
	723	ztwdmask(ji,jj) = 1._wp
	724	ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
	725	ztwdmask(ji,jj) = (tanh(50._wp( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1)) )
	726	ELSE
	727	ztwdmask(ji,jj) = 0._wp
	728	END IF
	729	END DO
	730	END DO
	731	ELSE
	732	DO jj = 1, jpj
	733	DO ji = 1, jpi ! vector opt.
	734	IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
	735	ztwdmask(ji,jj) = 1._wp
	736	ELSE
	737	ztwdmask(ji,jj) = 0._wp
	738	END IF
	739	END DO
	740	END DO
	741	END IF
	742
	743	END IF
	744
	745
[4292]	746	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
	747	DO ji = 2, fs_jpim1 ! Vector opt.
[9043]	748	zwx(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
[5836]	749	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
	750	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
[9043]	751	zwy(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
[5836]	752	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
	753	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
[4292]	754	END DO
	755	END DO
[5429]	756	CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
[4292]	757	!
[7753]	758	zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
	759	zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
[4370]	760	ELSE
[7753]	761	zhup2_e (:,:) = hu_n(:,:)
	762	zhvp2_e (:,:) = hv_n(:,:)
[4292]	763	ENDIF
	764	! !* after ssh
[1502]	765	! ! -----------
[4292]	766	! One should enforce volume conservation at open boundaries here
	767	! considering fluxes below:
	768	!
[7753]	769	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
	770	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
[4486]	771	!
	772	#if defined key_agrif
[6140]	773	! Set fluxes during predictor step to ensure volume conservation
	774	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
[4486]	775	IF((nbondi == -1).OR.(nbondi == 2)) THEN
[9019]	776	DO jj = 1, jpj
[9116]	777	zwx(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)
[4486]	778	END DO
	779	ENDIF
	780	IF((nbondi == 1).OR.(nbondi == 2)) THEN
	781	DO jj=1,jpj
[9116]	782	zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
[4486]	783	END DO
	784	ENDIF
	785	IF((nbondj == -1).OR.(nbondj == 2)) THEN
	786	DO ji=1,jpi
[9116]	787	zwy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)
[4486]	788	END DO
	789	ENDIF
	790	IF((nbondj == 1).OR.(nbondj == 2)) THEN
	791	DO ji=1,jpi
[9116]	792	zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
[4486]	793	END DO
	794	ENDIF
	795	ENDIF
	796	#endif
[9023]	797	IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
	798
	799	IF ( ln_wd_dl ) THEN
	800
	801
	802	! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells
	803
	804	DO jj = 1, jpjm1
	805	DO ji = 1, jpim1
	806	IF ( zwx(ji,jj) > 0.0 ) THEN
	807	zuwdmask(ji, jj) = ztwdmask(ji ,jj)
	808	ELSE
	809	zuwdmask(ji, jj) = ztwdmask(ji+1,jj)
	810	END IF
	811	zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj)
	812	un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj)
	813
	814	IF ( zwy(ji,jj) > 0.0 ) THEN
	815	zvwdmask(ji, jj) = ztwdmask(ji, jj )
	816	ELSE
	817	zvwdmask(ji, jj) = ztwdmask(ji, jj+1)
	818	END IF
	819	zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj)
	820	vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj)
	821	END DO
	822	END DO
	823
	824
	825	END IF
	826
[4486]	827	! Sum over sub-time-steps to compute advective velocities
	828	za2 = wgtbtp2(jn)
[7753]	829	un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
	830	vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
[9023]	831
	832	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True)
	833	IF ( ln_wd_dl_bc ) THEN
	834	zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:)
	835	zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:)
	836	END IF
	837
[4486]	838	! Set next sea level:
[4292]	839	DO jj = 2, jpjm1
[358]	840	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	841	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
[5836]	842	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj)
[358]	843	END DO
	844	END DO
[7753]	845	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
	846
[4292]	847	CALL lbc_lnk( ssha_e, 'T', 1._wp )
	848
[6140]	849	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
[7646]	850	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
[4292]	851	#if defined key_agrif
[6140]	852	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
[4292]	853	#endif
	854	!
	855	! Sea Surface Height at u-,v-points (vvl case only)
[6140]	856	IF( .NOT.ln_linssh ) THEN
[4292]	857	DO jj = 2, jpjm1
	858	DO ji = 2, jpim1 ! NO Vector Opt.
[9043]	859	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
[6140]	860	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
	861	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
[9043]	862	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
[6140]	863	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
	864	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
[4292]	865	END DO
[358]	866	END DO
[5429]	867	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
[4292]	868	ENDIF
	869	!
	870	! Half-step back interpolation of SSH for surface pressure computation:
	871	!----------------------------------------------------------------------
	872	IF ((jn==1).AND.ll_init) THEN
	873	za0=1._wp ! Forward-backward
	874	za1=0._wp
	875	za2=0._wp
	876	za3=0._wp
	877	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
	878	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
	879	za1=-0.1666666666666_wp ! za1 = gam
	880	za2= 0.0833333333333_wp ! za2 = eps
	881	za3= 0._wp
	882	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
[9023]	883	IF (rn_bt_alpha==0._wp) THEN
	884	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
	885	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
	886	za2=0.088_wp ! za2 = gam
	887	za3=0.013_wp ! za3 = eps
	888	ELSE
	889	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
	890	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
	891	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
	892	za1 = 1._wp - za0 - zgamma - zepsilon
	893	za2 = zgamma
	894	za3 = zepsilon
	895	ENDIF
[4292]	896	ENDIF
[6140]	897	!
[7753]	898	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
	899	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
[9023]	900	IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters
[6152]	901	DO jj = 2, jpjm1
[7646]	902	DO ji = 2, jpim1
	903	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
[9023]	904	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
	905	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) &
[7646]	906	& > rn_wdmin1 + rn_wdmin2
	907	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( &
	908	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
[9023]	909	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
[7646]	910
	911	IF(ll_tmp1) THEN
	912	zcpx(ji,jj) = 1.0_wp
	913	ELSE IF(ll_tmp2) THEN
	914	! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here
[9023]	915	zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
[7646]	916	& / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) )
	917	ELSE
	918	zcpx(ji,jj) = 0._wp
	919	END IF
	920
	921	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
[9023]	922	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
	923	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) &
[7646]	924	& > rn_wdmin1 + rn_wdmin2
	925	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( &
	926	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
[9023]	927	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
[7646]	928
	929	IF(ll_tmp1) THEN
	930	zcpy(ji,jj) = 1.0_wp
	931	ELSE IF(ll_tmp2) THEN
	932	! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here
[9023]	933	zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
[7646]	934	& / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) )
	935	ELSE
	936	zcpy(ji,jj) = 0._wp
	937	END IF
[6152]	938	END DO
[7646]	939	END DO
	940	END IF
[1502]	941	!
[4292]	942	! Compute associated depths at U and V points:
[9023]	943	IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form
[4292]	944	!
	945	DO jj = 2, jpjm1
	946	DO ji = 2, jpim1
[9043]	947	zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) &
[5836]	948	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) &
	949	& + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
[9043]	950	zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) &
[5836]	951	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) &
	952	& + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
[4292]	953	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
	954	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
	955	END DO
	956	END DO
[6152]	957
[4292]	958	ENDIF
	959	!
	960	! Add Coriolis trend:
[6140]	961	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
[4292]	962	! at each time step. We however keep them constant here for optimization.
	963	! Recall that zwx and zwy arrays hold fluxes at this stage:
	964	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
	965	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
	966	!
[6140]	967	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
[358]	968	DO jj = 2, jpjm1
	969	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	970	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
	971	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
	972	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
	973	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
[9043]	974	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	975	zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
[358]	976	END DO
	977	END DO
[508]	978	!
[6140]	979	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
[358]	980	DO jj = 2, jpjm1
	981	DO ji = fs_2, fs_jpim1 ! vector opt.
[9043]	982	zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
[5836]	983	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
[9043]	984	zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
[5836]	985	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
[4292]	986	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	987	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
[358]	988	END DO
	989	END DO
[508]	990	!
[6140]	991	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!
[358]	992	DO jj = 2, jpjm1
	993	DO ji = fs_2, fs_jpim1 ! vector opt.
[9043]	994	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
[5836]	995	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
	996	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
	997	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
[9043]	998	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
[5836]	999	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
	1000	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
	1001	& + ftne(ji,jj ) * zwx(ji ,jj ) )
[358]	1002	END DO
	1003	END DO
[508]	1004	!
[358]	1005	ENDIF
[4292]	1006	!
	1007	! Add tidal astronomical forcing if defined
[7646]	1008	IF ( ln_tide .AND. ln_tide_pot ) THEN
[4292]	1009	DO jj = 2, jpjm1
	1010	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	1011	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
	1012	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
[4292]	1013	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
	1014	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
	1015	END DO
	1016	END DO
	1017	ENDIF
	1018	!
[9023]	1019	! Add bottom stresses:
	1020	!jth do implicitly instead
	1021	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
[9045]	1022	DO jj = 2, jpjm1
	1023	DO ji = fs_2, fs_jpim1 ! vector opt.
	1024	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
	1025	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
	1026	END DO
	1027	END DO
[9023]	1028	ENDIF
[4292]	1029	!
	1030	! Surface pressure trend:
[9023]	1031	IF( ln_wd_il ) THEN
[6152]	1032	DO jj = 2, jpjm1
	1033	DO ji = 2, jpim1
	1034	! Add surface pressure gradient
	1035	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
	1036	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
[9023]	1037	zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj)
	1038	zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj)
[6152]	1039	END DO
	1040	END DO
	1041	ELSE
	1042	DO jj = 2, jpjm1
	1043	DO ji = fs_2, fs_jpim1 ! vector opt.
	1044	! Add surface pressure gradient
	1045	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
	1046	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
[9023]	1047	zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg
	1048	zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg
[6152]	1049	END DO
	1050	END DO
	1051	END IF
	1052
[4292]	1053	!
	1054	! Set next velocities:
[9023]	1055	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
[4292]	1056	DO jj = 2, jpjm1
	1057	DO ji = fs_2, fs_jpim1 ! vector opt.
[5930]	1058	ua_e(ji,jj) = ( un_e(ji,jj) &
[4292]	1059	& + rdtbt * ( zwx(ji,jj) &
	1060	& + zu_trd(ji,jj) &
	1061	& + zu_frc(ji,jj) ) &
[6140]	1062	& ) * ssumask(ji,jj)
[358]	1063
[5930]	1064	va_e(ji,jj) = ( vn_e(ji,jj) &
[4292]	1065	& + rdtbt * ( zwy(ji,jj) &
	1066	& + zv_trd(ji,jj) &
	1067	& + zv_frc(ji,jj) ) &
[6140]	1068	& ) * ssvmask(ji,jj)
[9023]	1069
	1070	!jth implicit bottom friction:
	1071	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
[9036]	1072	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
	1073	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
[9023]	1074	ENDIF
	1075
[4292]	1076	END DO
	1077	END DO
[6140]	1078	!
[9023]	1079	ELSE !* Flux form
[4292]	1080	DO jj = 2, jpjm1
	1081	DO ji = fs_2, fs_jpim1 ! vector opt.
[3294]	1082
[9023]	1083	zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
	1084	zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
	1085
[6152]	1086	zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
	1087	zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
	1088
[5930]	1089	ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) &
[4292]	1090	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
	1091	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
[6140]	1092	& + hu_n(ji,jj) * zu_frc(ji,jj) ) &
[4292]	1093	& ) * zhura
[358]	1094
[5930]	1095	va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) &
[4292]	1096	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
	1097	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
[6140]	1098	& + hv_n(ji,jj) * zv_frc(ji,jj) ) &
[4292]	1099	& ) * zhvra
[592]	1100	END DO
	1101	END DO
[4292]	1102	ENDIF
[9023]	1103
	1104
[6140]	1105	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
[9023]	1106	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
	1107	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
[7753]	1108	hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
	1109	hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
[1502]	1110	!
[1438]	1111	ENDIF
[6140]	1112	! !* domain lateral boundary
	1113	CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
[4292]	1114	!
[6140]	1115	! ! open boundaries
[7646]	1116	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
[4486]	1117	#if defined key_agrif
	1118	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
[4292]	1119	#endif
	1120	! !* Swap
	1121	! ! ----
[7753]	1122	ubb_e (:,:) = ub_e (:,:)
	1123	ub_e (:,:) = un_e (:,:)
	1124	un_e (:,:) = ua_e (:,:)
	1125	!
	1126	vbb_e (:,:) = vb_e (:,:)
	1127	vb_e (:,:) = vn_e (:,:)
	1128	vn_e (:,:) = va_e (:,:)
	1129	!
	1130	sshbb_e(:,:) = sshb_e(:,:)
	1131	sshb_e (:,:) = sshn_e(:,:)
	1132	sshn_e (:,:) = ssha_e(:,:)
[4292]	1133
	1134	! !* Sum over whole bt loop
	1135	! ! ----------------------
	1136	za1 = wgtbtp1(jn)
[6140]	1137	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
[7753]	1138	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
	1139	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
[9023]	1140	ELSE ! Sum transports
	1141	IF ( .NOT.ln_wd_dl ) THEN
	1142	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
	1143	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
	1144	ELSE
	1145	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
	1146	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
	1147	END IF
[4292]	1148	ENDIF
[9023]	1149	! ! Sum sea level
[7753]	1150	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
[9023]	1151
[358]	1152	! ! ==================== !
	1153	END DO ! end loop !
	1154	! ! ==================== !
[1438]	1155	! -----------------------------------------------------------------------------
[1502]	1156	! Phase 3. update the general trend with the barotropic trend
[1438]	1157	! -----------------------------------------------------------------------------
[1502]	1158	!
[4292]	1159	! Set advection velocity correction:
[9023]	1160	IF (ln_bt_fw) THEN
	1161	zwx(:,:) = un_adv(:,:)
	1162	zwy(:,:) = vn_adv(:,:)
	1163	IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
[9043]	1164	un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
	1165	vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
[9023]	1166	!
	1167	! Update corrective fluxes for next time step:
	1168	un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
	1169	vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
	1170	ELSE
	1171	un_bf(:,:) = 0._wp
	1172	vn_bf(:,:) = 0._wp
	1173	END IF
	1174	! Save integrated transport for next computation
[7753]	1175	ub2_b(:,:) = zwx(:,:)
	1176	vb2_b(:,:) = zwy(:,:)
[4292]	1177	ENDIF
[9023]	1178
	1179
[4292]	1180	!
	1181	! Update barotropic trend:
[6140]	1182	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
[4292]	1183	DO jk=1,jpkm1
[9043]	1184	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
	1185	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
[4292]	1186	END DO
	1187	ELSE
[5930]	1188	! At this stage, ssha has been corrected: compute new depths at velocity points
	1189	DO jj = 1, jpjm1
	1190	DO ji = 1, jpim1 ! NO Vector Opt.
[9043]	1191	zsshu_a(ji,jj) = r1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &
[5930]	1192	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
	1193	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
[9043]	1194	zsshv_a(ji,jj) = r1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &
[5930]	1195	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
	1196	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
	1197	END DO
	1198	END DO
	1199	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
	1200	!
[4292]	1201	DO jk=1,jpkm1
[9043]	1202	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
	1203	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
[4292]	1204	END DO
	1205	! Save barotropic velocities not transport:
[7753]	1206	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
	1207	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
[4292]	1208	ENDIF
[9023]	1209
	1210
	1211	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
[4292]	1212	DO jk = 1, jpkm1
[9023]	1213	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)r1_hu_n(:,:) - un_b(:,:) ) umask(:,:,jk)
	1214	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)r1_hv_n(:,:) - vn_b(:,:) ) vmask(:,:,jk)
[358]	1215	END DO
[9023]	1216
	1217	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
	1218	DO jk = 1, jpkm1
[9109]	1219	un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
	1220	& + zuwdav2(:,:)(un(:,:,jk) - un_adv(:,:)r1_hu_n(:,:)) ) * umask(:,:,jk)
	1221	vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) &
	1222	& + zvwdav2(:,:)(vn(:,:,jk) - vn_adv(:,:)r1_hv_n(:,:)) ) * vmask(:,:,jk)
[9023]	1223	END DO
	1224	END IF
	1225
	1226
	1227	CALL iom_put( "ubar", un_adv(:,:)*r1_hu_n(:,:) ) ! barotropic i-current
	1228	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv_n(:,:) ) ! barotropic i-current
[1502]	1229	!
[4486]	1230	#if defined key_agrif
	1231	! Save time integrated fluxes during child grid integration
[5656]	1232	! (used to update coarse grid transports at next time step)
[4486]	1233	!
[6140]	1234	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
	1235	IF( Agrif_NbStepint() == 0 ) THEN
[7753]	1236	ub2_i_b(:,:) = 0._wp
	1237	vb2_i_b(:,:) = 0._wp
[4486]	1238	END IF
	1239	!
	1240	za1 = 1._wp / REAL(Agrif_rhot(), wp)
[7753]	1241	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
	1242	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
[4486]	1243	ENDIF
	1244	#endif
[1502]	1245	! !* write time-spliting arrays in the restart
[6140]	1246	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
[508]	1247	!
[9023]	1248	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
	1249	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
[1662]	1250	!
[9019]	1251	IF( ln_diatmb ) THEN
[6140]	1252	CALL iom_put( "baro_u" , un_bumask(:,:,1)+zmdi(1-umask(:,:,1 ) ) ) ! Barotropic U Velocity
	1253	CALL iom_put( "baro_v" , vn_bvmask(:,:,1)+zmdi(1-vmask(:,:,1 ) ) ) ! Barotropic V Velocity
	1254	ENDIF
[2715]	1255	!
[508]	1256	END SUBROUTINE dyn_spg_ts
	1257
[6140]	1258
[4292]	1259	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
	1260	!!---------------------------------------------------------------------
	1261	!! * ROUTINE ts_wgt *
	1262	!!
	1263	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
	1264	!!----------------------------------------------------------------------
	1265	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
	1266	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
	1267	INTEGER, INTENT(inout) :: jpit ! cycle length
	1268	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
	1269	zwgt2 ! Secondary weights
	1270
	1271	INTEGER :: jic, jn, ji ! temporary integers
	1272	REAL(wp) :: za1, za2
	1273	!!----------------------------------------------------------------------
[508]	1274
[4292]	1275	zwgt1(:) = 0._wp
	1276	zwgt2(:) = 0._wp
	1277
	1278	! Set time index when averaged value is requested
	1279	IF (ll_fw) THEN
	1280	jic = nn_baro
	1281	ELSE
	1282	jic = 2 * nn_baro
	1283	ENDIF
	1284
	1285	! Set primary weights:
	1286	IF (ll_av) THEN
	1287	! Define simple boxcar window for primary weights
	1288	! (width = nn_baro, centered around jic)
	1289	SELECT CASE ( nn_bt_flt )
	1290	CASE( 0 ) ! No averaging
	1291	zwgt1(jic) = 1._wp
	1292	jpit = jic
	1293
	1294	CASE( 1 ) ! Boxcar, width = nn_baro
	1295	DO jn = 1, 3*nn_baro
	1296	za1 = ABS(float(jn-jic))/float(nn_baro)
	1297	IF (za1 < 0.5_wp) THEN
	1298	zwgt1(jn) = 1._wp
	1299	jpit = jn
	1300	ENDIF
	1301	ENDDO
	1302
	1303	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
	1304	DO jn = 1, 3*nn_baro
	1305	za1 = ABS(float(jn-jic))/float(nn_baro)
	1306	IF (za1 < 1._wp) THEN
	1307	zwgt1(jn) = 1._wp
	1308	jpit = jn
	1309	ENDIF
	1310	ENDDO
	1311	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
	1312	END SELECT
	1313
	1314	ELSE ! No time averaging
	1315	zwgt1(jic) = 1._wp
	1316	jpit = jic
	1317	ENDIF
	1318
	1319	! Set secondary weights
	1320	DO jn = 1, jpit
	1321	DO ji = jn, jpit
	1322	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
	1323	END DO
	1324	END DO
	1325
	1326	! Normalize weigths:
	1327	za1 = 1._wp / SUM(zwgt1(1:jpit))
	1328	za2 = 1._wp / SUM(zwgt2(1:jpit))
	1329	DO jn = 1, jpit
	1330	zwgt1(jn) = zwgt1(jn) * za1
	1331	zwgt2(jn) = zwgt2(jn) * za2
	1332	END DO
	1333	!
	1334	END SUBROUTINE ts_wgt
	1335
[6140]	1336
[508]	1337	SUBROUTINE ts_rst( kt, cdrw )
	1338	!!---------------------------------------------------------------------
	1339	!! * ROUTINE ts_rst *
	1340	!!
	1341	!! ** Purpose : Read or write time-splitting arrays in restart file
	1342	!!----------------------------------------------------------------------
	1343	INTEGER , INTENT(in) :: kt ! ocean time-step
	1344	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	1345	!
	1346	!!----------------------------------------------------------------------
	1347	!
	1348	IF( TRIM(cdrw) == 'READ' ) THEN
[4292]	1349	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) )
	1350	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) )
[9023]	1351	CALL iom_get( numror, jpdom_autoglo, 'un_bf' , un_bf (:,:) )
	1352	CALL iom_get( numror, jpdom_autoglo, 'vn_bf' , vn_bf (:,:) )
[4370]	1353	IF( .NOT.ln_bt_av ) THEN
[4292]	1354	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) )
	1355	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) )
	1356	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) )
	1357	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) )
	1358	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) )
	1359	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) )
[508]	1360	ENDIF
[4486]	1361	#if defined key_agrif
	1362	! Read time integrated fluxes
	1363	IF ( .NOT.Agrif_Root() ) THEN
	1364	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:) )
	1365	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:) )
	1366	ENDIF
	1367	#endif
[4292]	1368	!
	1369	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
	1370	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
	1371	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
[9023]	1372	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:) )
	1373	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:) )
[4292]	1374	!
	1375	IF (.NOT.ln_bt_av) THEN
	1376	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
	1377	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
	1378	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
	1379	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
	1380	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
	1381	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
	1382	ENDIF
[4486]	1383	#if defined key_agrif
	1384	! Save time integrated fluxes
	1385	IF ( .NOT.Agrif_Root() ) THEN
	1386	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
	1387	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
	1388	ENDIF
	1389	#endif
[4292]	1390	ENDIF
	1391	!
	1392	END SUBROUTINE ts_rst
[2528]	1393
[6140]	1394
	1395	SUBROUTINE dyn_spg_ts_init
[4292]	1396	!!---------------------------------------------------------------------
	1397	!! * ROUTINE dyn_spg_ts_init *
	1398	!!
	1399	!! ** Purpose : Set time splitting options
	1400	!!----------------------------------------------------------------------
[6140]	1401	INTEGER :: ji ,jj ! dummy loop indices
	1402	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
[9019]	1403	REAL(wp), DIMENSION(jpi,jpj) :: zcu
[4292]	1404	!!----------------------------------------------------------------------
[4370]	1405	!
[5930]	1406	! Max courant number for ext. grav. waves
[4370]	1407	!
[5930]	1408	DO jj = 1, jpj
	1409	DO ji =1, jpi
	1410	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
	1411	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
[7646]	1412	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
[4370]	1413	END DO
[5930]	1414	END DO
	1415	!
[5836]	1416	zcmax = MAXVAL( zcu(:,:) )
[4292]	1417	IF( lk_mpp ) CALL mpp_max( zcmax )
[2528]	1418
[4370]	1419	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
[6140]	1420	IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
[4292]	1421
[5836]	1422	rdtbt = rdt / REAL( nn_baro , wp )
[4292]	1423	zcmax = zcmax * rdtbt
[9023]	1424	! Print results
[4292]	1425	IF(lwp) WRITE(numout,*)
[9169]	1426	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
	1427	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
[5930]	1428	IF( ln_bt_auto ) THEN
[9169]	1429	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_baro '
[4370]	1430	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
[4292]	1431	ELSE
[9169]	1432	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist nn_baro = ', nn_baro
[358]	1433	ENDIF
[4292]	1434
	1435	IF(ln_bt_av) THEN
[9169]	1436	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_baro time steps is on '
[4292]	1437	ELSE
[9169]	1438	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
[4292]	1439	ENDIF
[508]	1440	!
[4292]	1441	!
	1442	IF(ln_bt_fw) THEN
[4370]	1443	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
[4292]	1444	ELSE
[4370]	1445	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
[4292]	1446	ENDIF
	1447	!
[4486]	1448	#if defined key_agrif
	1449	! Restrict the use of Agrif to the forward case only
[9023]	1450	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
[4486]	1451	#endif
	1452	!
[4370]	1453	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
[4292]	1454	SELECT CASE ( nn_bt_flt )
[6140]	1455	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
	1456	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
	1457	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
[9169]	1458	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
[4292]	1459	END SELECT
	1460	!
[4370]	1461	IF(lwp) WRITE(numout,*) ' '
	1462	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
	1463	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
	1464	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
	1465	!
[9023]	1466	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
	1467	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
	1468	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
	1469	ENDIF
	1470	!
[6140]	1471	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
[4292]	1472	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
	1473	ENDIF
[6140]	1474	IF( zcmax>0.9_wp ) THEN
[4292]	1475	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
	1476	ENDIF
	1477	!
[9124]	1478	! ! Allocate time-splitting arrays
	1479	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
	1480	!
	1481	! ! read restart when needed
	1482	IF( neuler /= 0 .AND. ln_bt_fw ) CALL ts_rst( nit000, 'READ' )
	1483	!
[4292]	1484	END SUBROUTINE dyn_spg_ts_init
[508]	1485
[358]	1486	!!======================================================================
	1487	END MODULE dynspg_ts

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

Download in other formats: