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

Context Navigation

source: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5910

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

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

Download in other formats: