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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4729

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

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

Download in other formats: