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

Context Navigation

source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5211

Last change on this file since 5211 was 5211, checked in by hliu, 9 years ago
add bottom friction limiter in dynspg_ts.F90
Property svn:keywords set to `Id`
File size: 61.0 KB

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

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

Download in other formats: