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

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4374

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

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

Download in other formats: