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 NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynspg_ts.F90 @ 11541

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

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

Download in other formats: