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 @ 6004

Last change on this file since 6004 was 6004, checked in by gm, 8 years ago
#1613: vvl by default, step III: Merge with the trunk (free surface simplification) (see wiki)
Property svn:keywords set to `Id`
File size: 54.9 KB

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

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

Download in other formats: