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

Context Navigation

source: branches/UKMO/antarctic_partial_slip/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5958

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

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

Download in other formats: