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

Context Navigation

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7037

Last change on this file since 7037 was 7037, checked in by mocavero, 8 years ago
ORCA2_LIM_PISCES hybrid version update
Property svn:keywords set to `Id`
File size: 68.1 KB

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

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

Download in other formats: