Context Navigation

dynspg_ts.F90 @ 8888

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

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

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

Context Navigation

source: branches/2017/dev_r8624_ENHANCE4_FREESURFACE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8888

Download in other formats: