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

Context Navigation

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5845

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

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

Download in other formats: