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/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg_ts.F90 @ 10919

Last change on this file since 10919 was 10919, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps :

Introduce ssh, uu_b, vv_b variables with time dimension (and pointers).
Convert dynspg modules (but no change to subcycling in dynspg_ts.F90).

Property svn:keywords set to Id

File size: 78.4 KB

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

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

Download in other formats: