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_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8214

Last change on this file since 8214 was 8214, checked in by gm, 7 years ago
#1883 (HPC-09) - correction of minor issues
Property svn:keywords set to `Id`
File size: 64.6 KB

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

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

Download in other formats: