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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynspg_ts.F90 in branches/UKMO/r5936_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/UKMO/r5936_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7113

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

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

Download in other formats: