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

Context Navigation

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 6667

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

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

Download in other formats: