Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

dynspg_ts.F90 @ 12359

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_PS44/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 12359

Download in other formats: