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 NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DYN/dynspg_ts.F90 @ 11380

Last change on this file since 11380 was 11380, checked in by girrmann, 5 years ago
dev_r10984_HPC-13 : adding extra halos in dyn_spg_ts is now possible, only works with a single halo when used with tide or bdy, see #2308
Property svn:keywords set to `Id`
File size: 96.4 KB

Rev	Line
[358]	1	MODULE dynspg_ts
[9023]	2
	3	!! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out !
	4
[358]	5	!!======================================================================
[6140]	6	!! * MODULE dynspg_ts *
	7	!! Ocean dynamics: surface pressure gradient trend, split-explicit scheme
	8	!!======================================================================
[1502]	9	!! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code
	10	!! - ! 2005-11 (V. Garnier, G. Madec) optimization
	11	!! - ! 2006-08 (S. Masson) distributed restart using iom
	12	!! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines
	13	!! - ! 2008-01 (R. Benshila) change averaging method
	14	!! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
[2528]	15	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
[2724]	16	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
[4292]	17	!! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping
	18	!! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility
[5930]	19	!! 3.7 ! 2015-11 (J. Chanut) free surface simplification
[7646]	20	!! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence
[9019]	21	!! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90)
[2724]	22	!!---------------------------------------------------------------------
[6140]	23
[358]	24	!!----------------------------------------------------------------------
[6140]	25	!! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme
	26	!! dyn_spg_ts_init: initialisation of the time-splitting scheme
	27	!! ts_wgt : set time-splitting weights for temporal averaging (or not)
	28	!! ts_rst : read/write time-splitting fields in restart file
[358]	29	!!----------------------------------------------------------------------
	30	USE oce ! ocean dynamics and tracers
	31	USE dom_oce ! ocean space and time domain
[888]	32	USE sbc_oce ! surface boundary condition: ocean
[9019]	33	USE zdf_oce ! vertical physics: variables
	34	USE zdfdrg ! vertical physics: top/bottom drag coef.
[5120]	35	USE sbcisf ! ice shelf variable (fwfisf)
[6140]	36	USE sbcapr ! surface boundary condition: atmospheric pressure
	37	USE dynadv , ONLY: ln_dynadv_vec
[9528]	38	USE dynvor ! vortivity scheme indicators
[358]	39	USE phycst ! physical constants
	40	USE dynvor ! vorticity term
[6152]	41	USE wet_dry ! wetting/drying flux limter
[7646]	42	USE bdy_oce ! open boundary
[10481]	43	USE bdyvol ! open boundary volume conservation
[5930]	44	USE bdytides ! open boundary condition data
[3294]	45	USE bdydyn2d ! open boundary conditions on barotropic variables
[4292]	46	USE sbctide ! tides
	47	USE updtide ! tide potential
[7646]	48	USE sbcwave ! surface wave
[9019]	49	USE diatmb ! Top,middle,bottom output
	50	#if defined key_agrif
[9570]	51	USE agrif_oce_interp ! agrif
[9124]	52	USE agrif_oce
[9019]	53	#endif
	54	#if defined key_asminc
	55	USE asminc ! Assimilation increment
	56	#endif
[6140]	57	!
	58	USE in_out_manager ! I/O manager
[358]	59	USE lib_mpp ! distributed memory computing library
	60	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	61	USE prtctl ! Print control
[2715]	62	USE iom ! IOM library
[4292]	63	USE restart ! only for lrst_oce
[9023]	64	USE diatmb ! Top,middle,bottom output
[358]	65
[11234]	66
[358]	67	IMPLICIT NONE
	68	PRIVATE
	69
[9124]	70	PUBLIC dyn_spg_ts ! called by dyn_spg
	71	PUBLIC dyn_spg_ts_init ! - - dyn_spg_init
[358]	72
[9019]	73	!! Time filtered arrays at baroclinic time step:
	74	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step
[9124]	75	!
[9023]	76	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
	77	REAL(wp),SAVE :: rdtbt ! Barotropic time step
[9019]	78	!
[11380]	79	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields
	80	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points
	81	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
	82	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
[4292]	83
[11380]	84	!! Arrays at barotropic time step: ! befbefore! before ! now ! after !
	85	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubb_e , ub_e , un_e , ua_e !: u-external velocity
	86	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vbb_e , vb_e , vn_e , va_e !: v-external velocity
	87	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, sshn_e, ssha_e !: external ssh
	88	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_e !: external u-depth
	89	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_e !: external v-depth
	90	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hur_e !: inverse of u-depth
	91	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hvr_e !: inverse of v-depth
	92	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_b , vb2_b !: Half step fluxes (ln_bt_fw=T)
	93	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_bf , vn_bf !: Asselin filtered half step fluxes (ln_bt_fw=T)
	94
	95	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0_xtd , hu_0_xtd , hv_0_xtd
	96	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_n_xtd , hv_n_xtd
	97	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_e1e2t_xtd, r1_e1e2u_xtd, r1_e1e2v_xtd
	98	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2t_xtd
	99	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: e2u_xtd , e1v_xtd
	100	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_e1u_xtd, r1_e2v_xtd
	101	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask_xtd, ssumask_xtd, ssvmask_xtd
	102	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff_t_xtd ! used in ENT scheme
	103
	104	#if defined key_agrif
	105	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_i_b, vb2_i_b !: Half step time integrated fluxes
	106	#endif
	107
[9043]	108	REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios
	109	REAL(wp) :: r1_8 = 0.125_wp !
	110	REAL(wp) :: r1_4 = 0.25_wp !
	111	REAL(wp) :: r1_2 = 0.5_wp !
[508]	112
[358]	113	!! * Substitutions
	114	# include "vectopt_loop_substitute.h90"
[2715]	115	!!----------------------------------------------------------------------
[9598]	116	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5217]	117	!! $Id$
[10068]	118	!! Software governed by the CeCILL license (see ./LICENSE)
[2715]	119	!!----------------------------------------------------------------------
[358]	120	CONTAINS
	121
[2715]	122	INTEGER FUNCTION dyn_spg_ts_alloc()
	123	!!----------------------------------------------------------------------
	124	!! * routine dyn_spg_ts_alloc *
	125	!!----------------------------------------------------------------------
[11380]	126	INTEGER :: ierr(6)
	127	INTEGER :: idbi, idei, idbj, idej ! lower/upper bounds of extended arrays
[4292]	128	!!----------------------------------------------------------------------
	129	ierr(:) = 0
[11380]	130	idbi = 1 - nn_hlts ; idbj = 1 - nn_hlts
	131	idei = jpi + nn_hlts ; idej = jpj + nn_hlts
[6140]	132	!
[11380]	133	ALLOCATE( wgtbtp1(3nn_baro), wgtbtp2(3nn_baro), zwz(idbi:idei,idbj:idej), STAT=ierr(1) )
	134	IF( ln_dynvor_een .OR. ln_dynvor_eeT ) THEN
	135	ALLOCATE( ftnw(idbi:idei,idbj:idej) , ftne(idbi:idei,idbj:idej) , ftsw(idbi:idei,idbj:idej) , ftse(idbi:idei,idbj:idej) &
	136	& , STAT=ierr(2) )
	137	ELSEIF( ln_dynvor_enT ) THEN
	138	ALLOCATE( ff_t_xtd(idbi:idei,idbj:idej), STAT=ierr(2) )
	139	END IF
[6140]	140	!
	141	ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) )
[11380]	142	ALLOCATE( ssha_e(idbi:idei,idbj:idej), sshn_e(idbi:idei,idbj:idej), sshb_e(idbi:idei,idbj:idej), sshbb_e(idbi:idei,idbj:idej) &
	143	& , ua_e(idbi:idei,idbj:idej), un_e(idbi:idei,idbj:idej), ub_e(idbi:idei,idbj:idej), ubb_e(idbi:idei,idbj:idej) &
	144	& , va_e(idbi:idei,idbj:idej), vn_e(idbi:idei,idbj:idej), vb_e(idbi:idei,idbj:idej), vbb_e(idbi:idei,idbj:idej) &
	145	& , hu_e(idbi:idei,idbj:idej), hur_e(idbi:idei,idbj:idej), hv_e(idbi:idei,idbj:idej), hvr_e(idbi:idei,idbj:idej) &
	146	& , STAT=ierr(4) )
	147	!
	148	ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_bf(jpi,jpj), vn_bf(jpi,jpj) , STAT=ierr(5) )
	149
	150	#if defined key_agrif
	151	ALLOCATE( ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj) , STAT=ierr(5) )
	152	#endif
	153	ALLOCATE( ht_0_xtd (idbi:idei,idbj:idej), hu_0_xtd (idbi:idei,idbj:idej), hv_0_xtd (idbi:idei,idbj:idej) &
	154	& , r1_e1e2t_xtd(idbi:idei,idbj:idej), r1_e1e2u_xtd(idbi:idei,idbj:idej), r1_e1e2v_xtd(idbi:idei,idbj:idej) &
	155	& , ssmask_xtd (idbi:idei,idbj:idej), ssumask_xtd (idbi:idei,idbj:idej), ssvmask_xtd (idbi:idei,idbj:idej) &
	156	& , e1e2t_xtd (idbi:idei,idbj:idej), e2u_xtd (idbi:idei,idbj:idej), e1v_xtd (idbi:idei,idbj:idej) &
	157	& , r1_e1u_xtd (idbi:idei,idbj:idej), r1_e2v_xtd (idbi:idei,idbj:idej) &
	158	& , STAT=ierr(6) )
[6140]	159	!
	160	dyn_spg_ts_alloc = MAXVAL( ierr(:) )
	161	!
[10425]	162	CALL mpp_sum( 'dynspg_ts', dyn_spg_ts_alloc )
	163	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_stop( 'STOP', 'dyn_spg_ts_alloc: failed to allocate arrays' )
[2715]	164	!
	165	END FUNCTION dyn_spg_ts_alloc
	166
[5836]	167
[358]	168	SUBROUTINE dyn_spg_ts( kt )
	169	!!----------------------------------------------------------------------
	170	!!
[6140]	171	!! ** Purpose : - Compute the now trend due to the explicit time stepping
	172	!! of the quasi-linear barotropic system, and add it to the
	173	!! general momentum trend.
[358]	174	!!
[6140]	175	!! ** Method : - split-explicit schem (time splitting) :
[4374]	176	!! Barotropic variables are advanced from internal time steps
	177	!! "n" to "n+1" if ln_bt_fw=T
	178	!! or from
	179	!! "n-1" to "n+1" if ln_bt_fw=F
	180	!! thanks to a generalized forward-backward time stepping (see ref. below).
[358]	181	!!
[4374]	182	!! ** Action :
	183	!! -Update the filtered free surface at step "n+1" : ssha
	184	!! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
[9023]	185	!! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv
[4374]	186	!! These are used to advect tracers and are compliant with discrete
	187	!! continuity equation taken at the baroclinic time steps. This
	188	!! ensures tracers conservation.
[6140]	189	!! - (ua, va) momentum trend updated with barotropic component.
[358]	190	!!
[6140]	191	!! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
[358]	192	!!---------------------------------------------------------------------
[1502]	193	INTEGER, INTENT(in) :: kt ! ocean time-step index
[2715]	194	!
[11380]	195	INTEGER :: ji, jj, jk, jm ! dummy loop indices
[9019]	196	LOGICAL :: ll_fw_start ! =T : forward integration
[9554]	197	LOGICAL :: ll_init ! =T : special startup of 2d equations
[11240]	198	INTEGER :: noffset ! local integers : time offset for bdy update
	199	REAL(wp) :: r1_2dt_b, z1_hu, z1_hv ! local scalars
[9528]	200	REAL(wp) :: za0, za1, za2, za3 ! - -
[11240]	201	REAL(wp) :: zmdi, zztmp, zldg ! - -
	202	REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - -
[11234]	203	REAL(wp) :: zun_save, zvn_save ! - -
[3294]	204	!
[11380]	205	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zu_trd, zssh_frc
	206	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zv_trd
	207	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zhU , zhV
	208	!
	209	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: hu_n_xtd, hv_n_xtd
	210	!
	211	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zu_spg , zv_spg
	212	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zu_frc , zv_frc
	213	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsshu_a, zhup2_e, zhtp2_e
	214	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsshv_a, zhvp2_e, zsshp2_e
	215	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zCdU_u , zCdU_v ! top/bottom stress at u- & v-points
	216
	217	!
[9023]	218	REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True.
	219
	220	INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point
	221
	222	REAL(wp) :: zepsilon, zgamma ! - -
[9019]	223	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef.
[9023]	224	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
	225	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask
[11380]	226
	227
	228	INTEGER :: idbi, idei, idbj, idej ! lower/upper bounds of extended arrays
	229	INTEGER :: ixtd ! number of halos over which the solution is currently correct
	230	INTEGER :: ibi, iei, ibj, iej ! lower and upper bounds over which the solution is currently correct
[358]	231	!!----------------------------------------------------------------------
[3294]	232	!
[11380]	233
	234
[9023]	235	IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) )
	236	! !* Allocate temporary arrays
[11380]	237	IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj) )
	238
	239	idbi = 1 - nn_hlts ; idbj = 1 - nn_hlts
	240	idei = jpi + nn_hlts ; idej = jpj + nn_hlts
	241	! ! allocate local arrays
	242	ALLOCATE( zu_spg (idbi:idei,idbj:idej), zv_spg (idbi:idei,idbj:idej) &
	243	& , zsshu_a (idbi:idei,idbj:idej), zsshv_a (idbi:idei,idbj:idej) &
	244	& , zhup2_e (idbi:idei,idbj:idej), zhvp2_e (idbi:idei,idbj:idej) &
	245	& , zCdU_u (idbi:idei,idbj:idej), zCdU_v (idbi:idei,idbj:idej) &
	246	& , zhtp2_e (idbi:idei,idbj:idej), zsshp2_e(idbi:idei,idbj:idej) &
	247	& , zu_trd (idbi:idei,idbj:idej), zu_frc (idbi:idei,idbj:idej) &
	248	& , zv_trd (idbi:idei,idbj:idej), zv_frc (idbi:idei,idbj:idej) &
	249	& , zhU (idbi:idei,idbj:idej), zhV (idbi:idei,idbj:idej) &
	250	& , zssh_frc(idbi:idei,idbj:idej) )
	251	! ! allocate redundant arrays
	252	ALLOCATE( hu_n_xtd(idbi:idei,idbj:idej), hv_n_xtd(idbi:idei,idbj:idej) )
[3294]	253	!
[6140]	254	zmdi=1.e+20 ! missing data indicator for masking
[9019]	255	!
[9023]	256	zwdramp = r_rn_wdmin1 ! simplest ramp
	257	! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
[11234]	258	! ! inverse of baroclinic time step
	259	IF( kt == nit000 .AND. neuler == 0 ) THEN ; r1_2dt_b = 1._wp / ( rdt )
	260	ELSE ; r1_2dt_b = 1._wp / ( 2._wp * rdt )
[4292]	261	ENDIF
	262	!
[9023]	263	ll_init = ln_bt_av ! if no time averaging, then no specific restart
[4292]	264	ll_fw_start = .FALSE.
[9023]	265	! ! time offset in steps for bdy data update
[6140]	266	IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro
	267	ELSE ; noffset = 0
	268	ENDIF
[4292]	269	!
[9023]	270	IF( kt == nit000 ) THEN !* initialisation
[508]	271	!
[358]	272	IF(lwp) WRITE(numout,*)
	273	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
	274	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
[4354]	275	IF(lwp) WRITE(numout,*)
[1502]	276	!
[6140]	277	IF( neuler == 0 ) ll_init=.TRUE.
[1502]	278	!
[6140]	279	IF( ln_bt_fw .OR. neuler == 0 ) THEN
	280	ll_fw_start =.TRUE.
	281	noffset = 0
[4292]	282	ELSE
[6140]	283	ll_fw_start =.FALSE.
[4292]	284	ENDIF
[11234]	285	! ! Set averaging weights and cycle length:
[6140]	286	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
[4292]	287	!
	288	ENDIF
	289	!
	290	! If forward start at previous time step, and centered integration,
	291	! then update averaging weights:
[5836]	292	IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
[4292]	293	ll_fw_start=.FALSE.
[9019]	294	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
[4292]	295	ENDIF
[11234]	296	!
[4292]	297
[358]	298	! -----------------------------------------------------------------------------
	299	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
	300	! -----------------------------------------------------------------------------
[1502]	301	!
[4292]	302	!
[11234]	303	! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends)
	304	! ! --------------------------- !
[11380]	305	zu_frc(1:jpi,1:jpj) = SUM( e3u_n(:,:,:) * ua(:,:,:) * umask(:,:,:) , DIM=3 ) * r1_hu_n(:,:)
	306	zv_frc(1:jpi,1:jpj) = SUM( e3v_n(:,:,:) * va(:,:,:) * vmask(:,:,:) , DIM=3 ) * r1_hv_n(:,:)
[1502]	307	!
[4292]	308	!
[11234]	309	! != Ua => baroclinic trend =! (remove its vertical mean)
	310	DO jk = 1, jpkm1 ! ------------------------ !
[11380]	311	ua(:,:,jk) = ( ua(:,:,jk) - zu_frc(1:jpi,1:jpj) ) * umask(:,:,jk)
	312	va(:,:,jk) = ( va(:,:,jk) - zv_frc(1:jpi,1:jpj) ) * vmask(:,:,jk)
[1502]	313	END DO
[7646]	314
	315	!!gm Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
	316	!!gm Is it correct to do so ? I think so...
	317
[11234]	318	! != remove 2D Coriolis and pressure gradient trends =!
	319	! ! ------------------------------------------------- !
[9528]	320	!
[11380]	321	! Set zwz, the barotropic Coriolis force coefficient
	322	! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes
	323	IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2d_init( idbi, idei, idbj, idej )
[1502]	324	!
[11234]	325	! !* 2D Coriolis trends
[11380]	326	zhU(1:jpi,1:jpj) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes
	327	zhV(1:jpi,1:jpj) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) ! NB: FULL domain : put a value in last row and column
[11234]	328	!
[11380]	329	! ! ht_n, hu_n, hv_n, un_b, vn_b are of size 1:jpi 1:jpj
	330	! ! zhU, zhV, zu_trd, zv_trd are of size idbi:idei idbj:idej
	331	CALL dyn_cor_2d( ht_n, hu_n, hv_n, un_b, vn_b, zhU, zhV, 1 , jpi , 1 , jpj & ! <<== in
	332	& , zu_trd, zv_trd, idbi, idei, idbj, idej ) ! ==>> out
[11234]	333	!
	334	IF( .NOT.ln_linssh ) THEN !* surface pressure gradient (variable volume only)
[508]	335	!
[11234]	336	IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg
	337	CALL wad_spg( sshn, zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy
[9528]	338	DO jj = 2, jpjm1
[11234]	339	DO ji = 2, jpim1 ! SPG with the application of W/D gravity filters
[9528]	340	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) &
	341	& * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth
	342	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) &
	343	& * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth
	344	END DO
	345	END DO
[11234]	346	ELSE ! now suface pressure gradient
[9019]	347	DO jj = 2, jpjm1
	348	DO ji = fs_2, fs_jpim1 ! vector opt.
	349	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj)
	350	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj)
	351	END DO
	352	END DO
	353	ENDIF
	354	!
[1502]	355	ENDIF
[9019]	356	!
[4292]	357	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
[358]	358	DO ji = fs_2, fs_jpim1
[6140]	359	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
	360	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
[3294]	361	END DO
[4292]	362	END DO
	363	!
[11234]	364	! != Add bottom stress contribution from baroclinic velocities =!
	365	! ! ----------------------------------------------------------- !
[11380]	366	CALL drg_init( idbi, idei, idbj, idej, zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients
	367	! ! arrays are computed on inner domain
[1502]	368	!
[11234]	369	! != Add atmospheric pressure forcing =!
	370	! ! ---------------------------------- !
	371	IF( ln_apr_dyn ) THEN
	372	IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
	373	DO jj = 2, jpjm1
[9019]	374	DO ji = fs_2, fs_jpim1 ! vector opt.
[11234]	375	zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
	376	zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
[9019]	377	END DO
	378	END DO
[11234]	379	ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
	380	zztmp = grav * r1_2
	381	DO jj = 2, jpjm1
[9019]	382	DO ji = fs_2, fs_jpim1 ! vector opt.
[11234]	383	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
	384	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
	385	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
	386	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
[9019]	387	END DO
	388	END DO
[11234]	389	ENDIF
[9019]	390	ENDIF
[11234]	391	!
	392	! != Add atmospheric pressure forcing =!
	393	! ! ---------------------------------- !
[9019]	394	IF( ln_bt_fw ) THEN ! Add wind forcing
	395	DO jj = 2, jpjm1
	396	DO ji = fs_2, fs_jpim1 ! vector opt.
	397	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)
	398	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)
	399	END DO
	400	END DO
[2724]	401	ELSE
[9043]	402	zztmp = r1_rau0 * r1_2
[9019]	403	DO jj = 2, jpjm1
	404	DO ji = fs_2, fs_jpim1 ! vector opt.
	405	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)
	406	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
	407	END DO
	408	END DO
[4292]	409	ENDIF
	410	!
[11234]	411	! !----------------!
	412	! !== sssh_frc ==! Right-Hand-Side of the barotropic ssh equation (over the FULL domain)
	413	! !----------------!
	414	! != Net water flux forcing applied to a water column =!
	415	! ! --------------------------------------------------- !
	416	IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
[11380]	417	zssh_frc(1:jpi,1:jpj) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
[11234]	418	ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
[9043]	419	zztmp = r1_rau0 * r1_2
[11380]	420	zssh_frc(1:jpi,1:jpj) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) + fwfisf(:,:) + fwfisf_b(:,:) )
[4292]	421	ENDIF
[11234]	422	! != Add Stokes drift divergence =! (if exist)
	423	IF( ln_sdw ) THEN ! ----------------------------- !
[11380]	424	zssh_frc(1:jpi,1:jpj) = zssh_frc(1:jpi,1:jpj) + div_sd(:,:)
[7646]	425	ENDIF
	426	!
[4292]	427	#if defined key_asminc
[11234]	428	! != Add the IAU weighted SSH increment =!
	429	! ! ------------------------------------ !
[4292]	430	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
[11380]	431	zssh_frc(1:jpi,1:jpj) = zssh_frc(1:jpi,1:jpj) - ssh_iau(:,:)
[4292]	432	ENDIF
	433	#endif
[11234]	434	! != Fill boundary data arrays for AGRIF
[5656]	435	! ! ------------------------------------
[4486]	436	#if defined key_agrif
	437	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
	438	#endif
[11380]	439
	440
[4292]	441	!
[358]	442	! -----------------------------------------------------------------------
[4292]	443	! Phase 2 : Integration of the barotropic equations
[358]	444	! -----------------------------------------------------------------------
[1502]	445	!
	446	! ! ==================== !
	447	! ! Initialisations !
[4292]	448	! ! ==================== !
[4370]	449	! Initialize barotropic variables:
[11380]	450	IF( ll_init ) THEN
[7753]	451	sshbb_e(:,:) = 0._wp
	452	ubb_e (:,:) = 0._wp
	453	vbb_e (:,:) = 0._wp
	454	sshb_e (:,:) = 0._wp
	455	ub_e (:,:) = 0._wp
	456	vb_e (:,:) = 0._wp
[4700]	457	ENDIF
	458	!
[11234]	459	IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0)
[11380]	460	zhtp2_e(1:jpi,1:jpj) = ht_n(1:jpi,1:jpj)
	461	zhup2_e(1:jpi,1:jpj) = hu_n(1:jpi,1:jpj)
	462	zhvp2_e(1:jpi,1:jpj) = hv_n(1:jpi,1:jpj)
[11234]	463	ENDIF
	464	!
[4370]	465	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
[11380]	466	sshn_e(1:jpi,1:jpj) = sshn(1:jpi,1:jpj)
	467	un_e (1:jpi,1:jpj) = un_b(1:jpi,1:jpj)
	468	vn_e (1:jpi,1:jpj) = vn_b(1:jpi,1:jpj)
[7753]	469	!
[11380]	470	hu_e (1:jpi,1:jpj) = hu_n(1:jpi,1:jpj)
	471	hv_e (1:jpi,1:jpj) = hv_n(1:jpi,1:jpj)
	472	hur_e (1:jpi,1:jpj) = r1_hu_n(1:jpi,1:jpj)
	473	hvr_e (1:jpi,1:jpj) = r1_hv_n(1:jpi,1:jpj)
[4370]	474	ELSE ! CENTRED integration: start from BEFORE fields
[11380]	475	sshn_e(1:jpi,1:jpj) = sshb(1:jpi,1:jpj)
	476	un_e (1:jpi,1:jpj) = ub_b(1:jpi,1:jpj)
	477	vn_e (1:jpi,1:jpj) = vb_b(1:jpi,1:jpj)
[7753]	478	!
[11380]	479	hu_e (1:jpi,1:jpj) = hu_b(1:jpi,1:jpj)
	480	hv_e (1:jpi,1:jpj) = hv_b(1:jpi,1:jpj)
	481	hur_e (1:jpi,1:jpj) = r1_hu_b(1:jpi,1:jpj)
	482	hvr_e (1:jpi,1:jpj) = r1_hv_b(1:jpi,1:jpj)
[4292]	483	ENDIF
	484	!
[11380]	485	hu_n_xtd(1:jpi,1:jpj) = hu_n(1:jpi,1:jpj)
	486	hv_n_xtd(1:jpi,1:jpj) = hv_n(1:jpi,1:jpj)
	487	!
	488	!
	489	! ! Extend arrays
	490	! ! --------------
	491	!
	492	IF( ln_linssh ) THEN
	493	CALL lbc_lnk_multi( 'dynspg_ts', hu_n_xtd, 'U', -1._wp, hv_n_xtd, 'V', -1._wp &
	494	& , zCdU_u , 'U', -1._wp, zCdU_v , 'V', -1._wp &
	495	& , zu_frc , 'U', -1._wp, zv_frc , 'V', -1._wp &
	496	& , un_e , 'U', -1._wp, vn_e , 'V', -1._wp &
	497	& , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp &
	498	& , hur_e , 'U', -1._wp, hvr_e , 'V', -1._wp &
	499	& , zhtp2_e , 'T', 1._wp, zhup2_e, 'U', -1._wp, zhvp2_e , 'V', -1._wp &
	500	& , zssh_frc, 'T', 1._wp, sshn_e , 'T', 1._wp, khlcom = nn_hls+nn_hlts )
	501	ELSE
	502	CALL lbc_lnk_multi( 'dynspg_ts', hu_n_xtd, 'U', -1._wp, hv_n_xtd, 'V', -1._wp &
	503	& , zCdU_u , 'U', -1._wp, zCdU_v , 'V', -1._wp &
	504	& , zu_frc , 'U', -1._wp, zv_frc , 'V', -1._wp &
	505	& , un_e , 'U', -1._wp, vn_e , 'V', -1._wp &
	506	& , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp &
	507	& , hur_e , 'U', -1._wp, hvr_e , 'V', -1._wp &
	508	& , zssh_frc, 'T', 1._wp, sshn_e , 'T', 1._wp, khlcom = nn_hls+nn_hlts )
	509	END IF
	510	!
[4292]	511	! Initialize sums:
[7753]	512	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
	513	va_b (:,:) = 0._wp
	514	ssha (:,:) = 0._wp ! Sum for after averaged sea level
	515	un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop
	516	vn_adv(:,:) = 0._wp
[9528]	517	!
	518	IF( ln_wd_dl ) THEN
[9023]	519	zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary)
	520	zvwdmask(:,:) = 0._wp !
[9528]	521	zuwdav2 (:,:) = 0._wp
	522	zvwdav2 (:,:) = 0._wp
[11380]	523	END IF
[9023]	524
[11380]	525	ixtd = nn_hls + nn_hlts ! solution is now correct over the whole domain (interior + regular halos + time splitting halos)
	526	ibi = 1 - nn_hlts ; ibj = 1 - nn_hlts
	527	iei = jpi + nn_hlts ; iej = jpj + nn_hlts
[9528]	528	! ! ==================== !
[11380]	529	DO jm = 1, icycle ! sub-time-step loop !
[1502]	530	! ! ==================== !
[10425]	531	!
[11380]	532	l_full_nf_update = jm == icycle ! false: disable full North fold update (performances) for jm = 1 to icycle-1
[11234]	533	!
	534	! !== Update the forcing ==! (BDY and tides)
	535	!
[11380]	536	IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jm, kt_offset= noffset+1 )
	537	IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jm, kt_offset= noffset )
[4292]	538	!
[11380]	539	! !== extrapolation at mid-step ==! (jm+1/2)
[11234]	540	!
	541	! !* Set extrapolation coefficients for predictor step:
[11380]	542	IF( (jm<3) .AND. ll_init ) THEN ! Forward
[4292]	543	za1 = 1._wp
	544	za2 = 0._wp
	545	za3 = 0._wp
	546	ELSE ! AB3-AM4 Coefficients: bet=0.281105
	547	za1 = 1.781105_wp ! za1 = 3/2 + bet
	548	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
	549	za3 = 0.281105_wp ! za3 = bet
	550	ENDIF
[11234]	551	!
[11380]	552	! !* Extrapolate barotropic velocities at mid-step (jm+1/2)
[11234]	553	!-- m+1/2 m m-1 m-2 --!
	554	!-- u = (3/2+beta) u -(1/2+2beta) u + beta u --!
	555	!-------------------------------------------------------------------------!
[11380]	556	ua_e(ibi:iei,ibj:iej) = za1 * un_e(ibi:iei,ibj:iej) + za2 * ub_e(ibi:iei,ibj:iej) + za3 * ubb_e(ibi:iei,ibj:iej)
	557	va_e(ibi:iei,ibj:iej) = za1 * vn_e(ibi:iei,ibj:iej) + za2 * vb_e(ibi:iei,ibj:iej) + za3 * vbb_e(ibi:iei,ibj:iej)
[4292]	558
[6140]	559	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
[4292]	560	! ! ------------------
[11380]	561	! !* Extrapolate Sea Level at mid-step (jm+1/2)
[11234]	562	!-- m+1/2 m m-1 m-2 --!
	563	!-- ssh = (3/2+beta) ssh -(1/2+2beta) ssh + beta ssh --!
	564	!--------------------------------------------------------------------------------!
[11380]	565	zsshp2_e(ibi:iei,ibj:iej) = za1 * sshn_e (ibi:iei,ibj:iej) + za2 * sshb_e(ibi:iei,ibj:iej) &
	566	& + za3 * sshbb_e(ibi:iei,ibj:iej)
[11234]	567	! set wetting & drying mask at tracer points for this barotropic mid-step
	568	IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask )
[9528]	569	!
[11240]	570	! ! ocean t-depth at mid-step
[11380]	571	zhtp2_e(ibi:iei,ibj:iej) = ht_0_xtd(ibi:iei,ibj:iej) + zsshp2_e(ibi:iei,ibj:iej)
[11240]	572	!
[11380]	573	! ! ocean u- and v-depth at mid-step
	574	DO jj = ibj, iej-1 ! not last column, not last row
	575	DO ji = ibi, iei-1
	576	zhup2_e(ji,jj) = hu_0_xtd(ji,jj) + r1_2 * r1_e1e2u_xtd(ji,jj) &
	577	& * ( e1e2t_xtd(ji ,jj) * zsshp2_e(ji ,jj) &
	578	& + e1e2t_xtd(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask_xtd(ji,jj)
	579	zhvp2_e(ji,jj) = hv_0_xtd(ji,jj) + r1_2 * r1_e1e2v_xtd(ji,jj) &
	580	& * ( e1e2t_xtd(ji,jj ) * zsshp2_e(ji,jj ) &
	581	& + e1e2t_xtd(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask_xtd(ji,jj)
[4292]	582	END DO
	583	END DO
	584	!
	585	ENDIF
	586	!
[11380]	587	! !== after SSH ==! (jm+1)
[11234]	588	!
	589	! ! update (ua_e,va_e) to enforce volume conservation at open boundaries
	590	! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
[11380]	591	IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jm, ua_e, va_e, zhup2_e, zhvp2_e )
[10481]	592	!
[11240]	593	! ! resulting flux at mid-step (not over the full domain)
[11380]	594	zhU(ibi:iei-1,ibj:iej-1) = e2u_xtd(ibi:iei-1,ibj:iej-1) * ua_e(ibi:iei-1,ibj:iej-1) * zhup2_e(ibi:iei-1,ibj:iej-1)
	595	zhV(ibi:iei-1,ibj:iej-1) = e1v_xtd(ibi:iei-1,ibj:iej-1) * va_e(ibi:iei-1,ibj:iej-1) * zhvp2_e(ibi:iei-1,ibj:iej-1)
[4486]	596	!
	597	#if defined key_agrif
[6140]	598	! Set fluxes during predictor step to ensure volume conservation
	599	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
[4486]	600	IF((nbondi == -1).OR.(nbondi == 2)) THEN
[9019]	601	DO jj = 1, jpj
[11234]	602	zhU(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)
	603	zhV(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)
[4486]	604	END DO
	605	ENDIF
	606	IF((nbondi == 1).OR.(nbondi == 2)) THEN
	607	DO jj=1,jpj
[11234]	608	zhU(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
	609	zhV(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj)
[4486]	610	END DO
	611	ENDIF
	612	IF((nbondj == -1).OR.(nbondj == 2)) THEN
	613	DO ji=1,jpi
[11234]	614	zhV(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)
	615	zhU(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)
[4486]	616	END DO
	617	ENDIF
	618	IF((nbondj == 1).OR.(nbondj == 2)) THEN
	619	DO ji=1,jpi
[11234]	620	zhV(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
	621	zhU(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1)
[4486]	622	END DO
	623	ENDIF
	624	ENDIF
	625	#endif
[11240]	626	IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rdtbt) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV
[9023]	627
[11234]	628	IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where
	629	! ! the direction of the flow is from dry cells
[11240]	630	CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V
[9528]	631	!
	632	ENDIF
[11380]	633	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
	634	IF ( ln_wd_dl_bc ) THEN
	635	zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column
	636	zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row
	637	END IF
	638	!
	639	! Sum over sub-time-steps to compute advective velocities (only correct on interior domain)
	640	za2 = wgtbtp2(jm) ! zhU, zhV hold fluxes extrapolated at jm+1/2
	641	un_adv(1:jpi,1:jpj) = un_adv(1:jpi,1:jpj) + za2 * zhU(1:jpi,1:jpj) * r1_e2u(1:jpi,1:jpj)
	642	vn_adv(1:jpi,1:jpj) = vn_adv(1:jpi,1:jpj) + za2 * zhV(1:jpi,1:jpj) * r1_e1v(1:jpi,1:jpj)
[11240]	643	!
	644	!
[11234]	645	! Compute Sea Level at step jit+1
	646	!-- m+1 m m+1/2 --!
	647	!-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --!
	648	!-------------------------------------------------------------------------!
[11380]	649	! correct domain reduction
	650	ixtd = ixtd - 1
	651	ibi = ibi + 1 ; ibj = ibj + 1
	652	iei = iei - 1 ; iej = iej - 1
	653	DO jj = ibj, iej
	654	DO ji = ibi, iei
	655	zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t_xtd(ji,jj)
	656	ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask_xtd(ji,jj)
[11240]	657	END DO
	658	END DO
	659	!
[11380]	660	IF( nn_hlts == 0 ) THEN
	661	CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )
	662	ixtd = nn_hls + nn_hlts ! solution is now correct over the whole domain
	663	ibi = 1 - nn_hlts ; ibj = 1 - nn_hlts
	664	iei = jpi + nn_hlts ; iej = jpj + nn_hlts
[11372]	665	END IF
[11380]	666
[11240]	667	!
[6140]	668	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
[11380]	669	IF( ln_bdy ) THEN
	670	CALL swap_bdyptr ! bdy treatment is now done on extended domain
	671	CALL bdy_ssh( ssha_e, idbi, idei, idbj, idej, ldcomall=.true., pmask=ssmask_xtd, khlcom=nn_hls+nn_hlts )
	672	CALL swap_bdyptr ! bdy treatment is now done on regular domain
	673	END IF
	674
[4292]	675	#if defined key_agrif
[11380]	676	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jm )
[4292]	677	#endif
[11380]	678	!
	679	! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
	680	!-- m+1/2 m+1 m m-1 m-2 --!
	681	!-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --!
	682	!-----------------------------------------------------------------------------------------!
	683	CALL ts_bck_interp( jm, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation
	684	zsshp2_e(ibi:iei,ibj:iej) = za0 * ssha_e(ibi:iei,ibj:iej) + za1 * sshn_e (ibi:iei,ibj:iej) &
	685	& + za2 * sshb_e(ibi:iei,ibj:iej) + za3 * sshbb_e(ibi:iei,ibj:iej)
	686	!
	687	!
[4292]	688	! Sea Surface Height at u-,v-points (vvl case only)
[11380]	689	IF( .NOT.ln_linssh ) THEN
	690	DO jj = ibj, iej-1
	691	DO ji = ibi, iei-1
	692	zsshu_a(ji,jj) = r1_2 * ssumask_xtd(ji,jj) * r1_e1e2u_xtd(ji,jj) &
	693	& * ( e1e2t_xtd(ji ,jj ) * ssha_e(ji ,jj ) &
	694	& + e1e2t_xtd(ji+1,jj ) * ssha_e(ji+1,jj ) )
	695	zsshv_a(ji,jj) = r1_2 * ssvmask_xtd(ji,jj) * r1_e1e2v_xtd(ji,jj) &
	696	& * ( e1e2t_xtd(ji ,jj ) * ssha_e(ji ,jj ) &
	697	& + e1e2t_xtd(ji ,jj+1) * ssha_e(ji ,jj+1) )
[4292]	698	END DO
[358]	699	END DO
[11380]	700	ENDIF
[1502]	701	!
[11234]	702	! ! Surface pressure gradient
	703	zldg = ( 1._wp - rn_scal_load ) * grav ! local factor
[11380]	704	DO jj = ibj, iej-1
	705	DO ji = ibi, iei-1
	706	zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u_xtd(ji,jj)
	707	zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v_xtd(ji,jj)
[4292]	708	END DO
[11234]	709	END DO
	710	IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient
	711	CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters
	712	zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
	713	zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
[4292]	714	ENDIF
	715	!
	716	! Add Coriolis trend:
[11380]	717	! - zwz array used in dyn_cor_2d or triads normally depend on sea level with ln_linssh=F and should be updated
[4292]	718	! at each time step. We however keep them constant here for optimization.
[11380]	719	! - Recall that zhU and zhV hold fluxes at jm+1/2 (extrapolated not backward interpolated)
	720	! - zu_trd_xtd and zv_trd_xtd are only correct on (ibi+1:iei-1,ibj+1:iej-1)
	721	! NOTE : input flux arguments have to be correct (ibi:iei,ibj:iej) -> a lbc call between input arguments computation
	722	! and this call without fluxes (typically after ssh at step m+1 computation) would not yield correct results
	723	CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV , idbi, idei, idbj, idej &
	724	& , zu_trd, zv_trd , idbi, idei, idbj, idej )
[4292]	725	!
	726	! Add tidal astronomical forcing if defined
[11380]	727	! pot_astro is correct on 1:jpi,1:jpj
[7646]	728	IF ( ln_tide .AND. ln_tide_pot ) THEN
[4292]	729	DO jj = 2, jpjm1
	730	DO ji = fs_2, fs_jpim1 ! vector opt.
[11234]	731	zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
	732	zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
[4292]	733	END DO
	734	END DO
	735	ENDIF
	736	!
[9023]	737	! Add bottom stresses:
	738	!jth do implicitly instead
	739	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
[11380]	740	DO jj = ibj, iej
	741	DO ji = ibi, iei
[9045]	742	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
	743	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
	744	END DO
	745	END DO
[11234]	746	ENDIF
[4292]	747	!
	748	! Set next velocities:
[11234]	749	! Compute barotropic speeds at step jit+1 (h : total height of the water colomn)
	750	!-- VECTOR FORM
[11261]	751	!-- m+1 m / m+1/2 \ --!
	752	!-- u = u + delta_t' * \ (1-r)g grad_x( ssh') - f * k vect u + frc / --!
	753	!-- --!
[11380]	754	!-- FLUX FORM --!
[11261]	755	!-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --!
	756	!-- u = m+1 \| h * u + delta_t' * \ h * (1-r)g grad_x( ssh') - h * f * k vect u + h * frc / \| --!
	757	!-- h \ / --!
	758	!------------------------------------------------------------------------------------------------------------------------!
[11380]	759	! correct domain reduction
	760	ixtd = ixtd - 1
	761	ibi = ibi + 1 ; ibj = ibj + 1
	762	iei = iei - 1 ; iej = iej - 1
	763	!
[9023]	764	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
[11380]	765	DO jj = ibj, iej
	766	DO ji = ibi, iei
[5930]	767	ua_e(ji,jj) = ( un_e(ji,jj) &
[11234]	768	& + rdtbt * ( zu_spg(ji,jj) &
[4292]	769	& + zu_trd(ji,jj) &
	770	& + zu_frc(ji,jj) ) &
[11380]	771	& ) * ssumask_xtd(ji,jj)
[358]	772
[5930]	773	va_e(ji,jj) = ( vn_e(ji,jj) &
[11234]	774	& + rdtbt * ( zv_spg(ji,jj) &
[4292]	775	& + zv_trd(ji,jj) &
	776	& + zv_frc(ji,jj) ) &
[11380]	777	& ) * ssvmask_xtd(ji,jj)
[4292]	778	END DO
	779	END DO
[6140]	780	!
[9023]	781	ELSE !* Flux form
[11380]	782	DO jj = ibj, iej
	783	DO ji = ibi, iei
	784	! ! hu_e, hv_e hold depth at jm, zhup2_e, zhvp2_e hold extrapolated depth at jm+1/2
	785	! ! backward interpolated depth used in spg terms at jm+1/2
	786	zhu_bck = hu_0_xtd(ji,jj) + r1_2r1_e1e2u_xtd(ji,jj) &
	787	& ( e1e2t_xtd(ji ,jj) * zsshp2_e(ji ,jj) &
	788	& + e1e2t_xtd(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask_xtd(ji,jj)
	789	zhv_bck = hv_0_xtd(ji,jj) + r1_2r1_e1e2v_xtd(ji,jj) &
	790	& ( e1e2t_xtd(ji,jj ) * zsshp2_e(ji,jj ) &
	791	& + e1e2t_xtd(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask_xtd(ji,jj)
	792	! ! inverse depth at jm+1
	793	z1_hu = ssumask_xtd(ji,jj) / ( hu_0_xtd(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask_xtd(ji,jj) )
	794	z1_hv = ssvmask_xtd(ji,jj) / ( hv_0_xtd(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask_xtd(ji,jj) )
[11234]	795	!
[11380]	796	ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) &
	797	& + rdtbt * ( zhu_bck * zu_spg(ji,jj) & !
	798	& + zhup2_e(ji,jj) * zu_trd(ji,jj) & !
	799	& + hu_n_xtd (ji,jj) * zu_frc(ji,jj) ) ) * z1_hu
[11234]	800	!
[11380]	801	va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) &
	802	& + rdtbt * ( zhv_bck * zv_spg(ji,jj) & !
	803	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) & !
	804	& + hv_n_xtd (ji,jj) * zv_frc(ji,jj) ) ) * z1_hv
[592]	805	END DO
	806	END DO
[4292]	807	ENDIF
[10272]	808	!jth implicit bottom friction:
	809	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
	810	DO jj = 2, jpjm1
	811	DO ji = fs_2, fs_jpim1 ! vector opt.
	812	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
	813	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
	814	END DO
	815	END DO
	816	ENDIF
[11241]	817
[11380]	818	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) and inverse depth
	819	hu_e (ibi:iei,ibj:iej) = hu_0_xtd(ibi:iei,ibj:iej) + zsshu_a(ibi:iei,ibj:iej)
	820	hv_e (ibi:iei,ibj:iej) = hv_0_xtd(ibi:iei,ibj:iej) + zsshv_a(ibi:iei,ibj:iej)
	821	hur_e(ibi:iei,ibj:iej) = ssumask_xtd(ibi:iei,ibj:iej) / ( hu_e(ibi:iei,ibj:iej) + 1._wp - ssumask_xtd(ibi:iei,ibj:iej) )
	822	hvr_e(ibi:iei,ibj:iej) = ssvmask_xtd(ibi:iei,ibj:iej) / ( hv_e(ibi:iei,ibj:iej) + 1._wp - ssvmask_xtd(ibi:iei,ibj:iej) )
[1438]	823	ENDIF
[11380]	824
	825	IF( ixtd == 0 ) THEN
	826	IF( .NOT. ln_linssh ) THEN
	827	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & ! after
	828	& , un_e , 'U', -1._wp, vn_e , 'V', -1._wp & ! now
	829	& , ub_e , 'U', -1._wp, vb_e , 'V', -1._wp & ! before
	830	& , ubb_e , 'U', -1._wp, vbb_e , 'V', -1._wp & ! before before
	831	& , ssha_e, 'T', 1._wp, sshn_e , 'T', 1._wp & ! after, now
	832	& , sshb_e, 'T', 1._wp, sshbb_e, 'T', 1._wp & ! before, before before
	833	& , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp &
	834	& , hur_e , 'U', -1._wp, hvr_e , 'V', -1._wp &
	835	& , khlcom = nn_hls+nn_hlts )
	836	ELSE
	837	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & ! after
	838	& , un_e , 'U', -1._wp, vn_e , 'V', -1._wp & ! now
	839	& , ub_e , 'U', -1._wp, vb_e , 'V', -1._wp & ! before
	840	& , ubb_e , 'U', -1._wp, vbb_e , 'V', -1._wp & ! before before
	841	& , ssha_e, 'T', 1._wp, sshn_e , 'T', 1._wp & ! after, now
	842	& , sshb_e, 'T', 1._wp, sshbb_e, 'T', 1._wp & ! before, before before
	843	& , khlcom = nn_hls+nn_hlts )
	844	END IF
	845	ixtd = nn_hls + nn_hlts ! solution is now correct over the whole domain
	846	ibi = 1 - nn_hlts ; ibj = 1 - nn_hlts
	847	iei = jpi + nn_hlts ; iej = jpj + nn_hlts
	848	END IF
[4292]	849	!
[11261]	850	!
[11380]	851	! ! open boundaries
	852	! ! bdy treatment is here done on regular domain (nn_hlts forced to 1 if ln_bdy or ln_tides)
	853	IF( ln_bdy ) CALL bdy_dyn2d( jm, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e, idbi, idei, idbj, idej &
	854	& , ldcomall=.true., pumask=ssumask_xtd, pvmask=ssvmask_xtd, khlcom=nn_hls+nn_hlts )
[4486]	855	#if defined key_agrif
[11380]	856	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jm ) ! Agrif
[4292]	857	#endif
	858	! !* Swap
	859	! ! ----
[7753]	860	ubb_e (:,:) = ub_e (:,:)
	861	ub_e (:,:) = un_e (:,:)
	862	un_e (:,:) = ua_e (:,:)
	863	!
	864	vbb_e (:,:) = vb_e (:,:)
	865	vb_e (:,:) = vn_e (:,:)
	866	vn_e (:,:) = va_e (:,:)
	867	!
	868	sshbb_e(:,:) = sshb_e(:,:)
	869	sshb_e (:,:) = sshn_e(:,:)
	870	sshn_e (:,:) = ssha_e(:,:)
[4292]	871
	872	! !* Sum over whole bt loop
	873	! ! ----------------------
[11380]	874	za1 = wgtbtp1(jm)
[6140]	875	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
[11380]	876	ua_b(1:jpi,1:jpj) = ua_b(1:jpi,1:jpj) + za1 * ua_e(1:jpi,1:jpj)
	877	va_b(1:jpi,1:jpj) = va_b(1:jpi,1:jpj) + za1 * va_e(1:jpi,1:jpj)
[9023]	878	ELSE ! Sum transports
	879	IF ( .NOT.ln_wd_dl ) THEN
[11380]	880	ua_b(1:jpi,1:jpj) = ua_b(1:jpi,1:jpj) + za1 * ua_e(1:jpi,1:jpj) * hu_e(1:jpi,1:jpj)
	881	va_b(1:jpi,1:jpj) = va_b(1:jpi,1:jpj) + za1 * va_e(1:jpi,1:jpj) * hv_e(1:jpi,1:jpj)
[9023]	882	ELSE
[11380]	883	ua_b(1:jpi,1:jpj) = ua_b(1:jpi,1:jpj) + za1 * ua_e(1:jpi,1:jpj) * hu_e(1:jpi,1:jpj) * zuwdmask(1:jpi,1:jpj)
	884	va_b(1:jpi,1:jpj) = va_b(1:jpi,1:jpj) + za1 * va_e(1:jpi,1:jpj) * hv_e(1:jpi,1:jpj) * zvwdmask(1:jpi,1:jpj)
[9023]	885	END IF
[4292]	886	ENDIF
[9023]	887	! ! Sum sea level
[11380]	888	ssha(1:jpi,1:jpj) = ssha(1:jpi,1:jpj) + za1 * ssha_e(1:jpi,1:jpj)
[9023]	889
[358]	890	! ! ==================== !
	891	END DO ! end loop !
	892	! ! ==================== !
[1438]	893	! -----------------------------------------------------------------------------
[1502]	894	! Phase 3. update the general trend with the barotropic trend
[1438]	895	! -----------------------------------------------------------------------------
[11380]	896	! Correction on regular halos
	897	CALL lbc_lnk_multi( 'dynspg_ts', un_adv, 'U', -1._wp, vn_adv, 'V', -1._wp &
	898	& , ua_b , 'U', -1._wp, va_b , 'V', -1._wp &
	899	& , ssha , 'T', -1._wp )
[1502]	900	!
[4292]	901	! Set advection velocity correction:
[9023]	902	IF (ln_bt_fw) THEN
	903	IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
[11234]	904	DO jj = 1, jpj
	905	DO ji = 1, jpi
	906	zun_save = un_adv(ji,jj)
	907	zvn_save = vn_adv(ji,jj)
	908	! ! apply the previously computed correction
	909	un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) )
	910	vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) )
	911	! ! Update corrective fluxes for next time step
	912	un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
	913	vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
	914	! ! Save integrated transport for next computation
	915	ub2_b(ji,jj) = zun_save
	916	vb2_b(ji,jj) = zvn_save
	917	END DO
	918	END DO
[9023]	919	ELSE
[11234]	920	un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero
	921	vn_bf(:,:) = 0._wp
	922	ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation
	923	vb2_b(:,:) = vn_adv(:,:)
	924	END IF
[4292]	925	ENDIF
[9023]	926
	927
[4292]	928	!
	929	! Update barotropic trend:
[6140]	930	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
[4292]	931	DO jk=1,jpkm1
[9043]	932	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
	933	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
[4292]	934	END DO
	935	ELSE
[5930]	936	! At this stage, ssha has been corrected: compute new depths at velocity points
	937	DO jj = 1, jpjm1
	938	DO ji = 1, jpim1 ! NO Vector Opt.
[9554]	939	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
	940	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
[5930]	941	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
[9554]	942	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
	943	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
[5930]	944	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
	945	END DO
[11380]	946	END DO ! Boundary conditions
	947	CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp, khlcom=nn_hls+nn_hlts ) ! change array used?
[5930]	948	!
[4292]	949	DO jk=1,jpkm1
[9043]	950	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
	951	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
[4292]	952	END DO
	953	! Save barotropic velocities not transport:
[11380]	954	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(1:jpi,1:jpj) + 1._wp - ssumask(:,:) )
	955	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(1:jpi,1:jpj) + 1._wp - ssvmask(:,:) )
[4292]	956	ENDIF
[9023]	957
	958
	959	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
[4292]	960	DO jk = 1, jpkm1
[9023]	961	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)r1_hu_n(:,:) - un_b(:,:) ) umask(:,:,jk)
	962	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)r1_hv_n(:,:) - vn_b(:,:) ) vmask(:,:,jk)
[358]	963	END DO
[9023]	964
	965	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
	966	DO jk = 1, jpkm1
[9109]	967	un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
	968	& + zuwdav2(:,:)(un(:,:,jk) - un_adv(:,:)r1_hu_n(:,:)) ) * umask(:,:,jk)
	969	vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) &
	970	& + zvwdav2(:,:)(vn(:,:,jk) - vn_adv(:,:)r1_hv_n(:,:)) ) * vmask(:,:,jk)
[9023]	971	END DO
	972	END IF
	973
	974
	975	CALL iom_put( "ubar", un_adv(:,:)*r1_hu_n(:,:) ) ! barotropic i-current
	976	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv_n(:,:) ) ! barotropic i-current
[1502]	977	!
[4486]	978	#if defined key_agrif
	979	! Save time integrated fluxes during child grid integration
[5656]	980	! (used to update coarse grid transports at next time step)
[4486]	981	!
[6140]	982	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
	983	IF( Agrif_NbStepint() == 0 ) THEN
[7753]	984	ub2_i_b(:,:) = 0._wp
	985	vb2_i_b(:,:) = 0._wp
[4486]	986	END IF
	987	!
	988	za1 = 1._wp / REAL(Agrif_rhot(), wp)
[7753]	989	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
	990	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
[4486]	991	ENDIF
	992	#endif
[1502]	993	! !* write time-spliting arrays in the restart
[6140]	994	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
[508]	995	!
[9023]	996	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
	997	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
[1662]	998	!
[9019]	999	IF( ln_diatmb ) THEN
[9554]	1000	CALL iom_put( "baro_u" , un_bssumask(:,:)+zmdi(1.-ssumask(:,:) ) ) ! Barotropic U Velocity
	1001	CALL iom_put( "baro_v" , vn_bssvmask(:,:)+zmdi(1.-ssvmask(:,:) ) ) ! Barotropic V Velocity
[6140]	1002	ENDIF
[2715]	1003	!
[11380]	1004
	1005	! deallocate temporary arrays
	1006	DEALLOCATE( zu_trd , zv_trd &
	1007	& , zu_frc , zv_frc &
	1008	& , zu_spg , zv_spg &
	1009	& , zsshu_a , zsshv_a &
	1010	& , zhup2_e , zhvp2_e &
	1011	& , zCdU_u , zCdU_v &
	1012	& , zhU , zhV &
	1013	& , zssh_frc, zsshp2_e &
	1014	& , zhtp2_e &
	1015	& , hu_n_xtd, hv_n_xtd )
	1016	!
[508]	1017	END SUBROUTINE dyn_spg_ts
	1018
[6140]	1019
[4292]	1020	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
	1021	!!---------------------------------------------------------------------
	1022	!! * ROUTINE ts_wgt *
	1023	!!
	1024	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
	1025	!!----------------------------------------------------------------------
[11380]	1026	LOGICAL, INTENT(in ) :: ll_av ! temporal averaging=.true.
	1027	LOGICAL, INTENT(in ) :: ll_fw ! forward time splitting =.true.
	1028	INTEGER, INTENT(inout) :: jpit ! cycle length
[4292]	1029	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
	1030	zwgt2 ! Secondary weights
	1031
[11380]	1032	INTEGER :: jic, jm, ji ! temporary integers
[4292]	1033	REAL(wp) :: za1, za2
	1034	!!----------------------------------------------------------------------
[508]	1035
[4292]	1036	zwgt1(:) = 0._wp
	1037	zwgt2(:) = 0._wp
	1038
	1039	! Set time index when averaged value is requested
	1040	IF (ll_fw) THEN
	1041	jic = nn_baro
	1042	ELSE
	1043	jic = 2 * nn_baro
	1044	ENDIF
	1045
	1046	! Set primary weights:
	1047	IF (ll_av) THEN
	1048	! Define simple boxcar window for primary weights
	1049	! (width = nn_baro, centered around jic)
	1050	SELECT CASE ( nn_bt_flt )
	1051	CASE( 0 ) ! No averaging
	1052	zwgt1(jic) = 1._wp
	1053	jpit = jic
	1054
	1055	CASE( 1 ) ! Boxcar, width = nn_baro
[11380]	1056	DO jm = 1, 3*nn_baro
	1057	za1 = ABS(float(jm-jic))/float(nn_baro)
[4292]	1058	IF (za1 < 0.5_wp) THEN
[11380]	1059	zwgt1(jm) = 1._wp
	1060	jpit = jm
[4292]	1061	ENDIF
	1062	ENDDO
	1063
	1064	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
[11380]	1065	DO jm = 1, 3*nn_baro
	1066	za1 = ABS(float(jm-jic))/float(nn_baro)
[4292]	1067	IF (za1 < 1._wp) THEN
[11380]	1068	zwgt1(jm) = 1._wp
	1069	jpit = jm
[4292]	1070	ENDIF
	1071	ENDDO
	1072	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
	1073	END SELECT
	1074
	1075	ELSE ! No time averaging
	1076	zwgt1(jic) = 1._wp
	1077	jpit = jic
	1078	ENDIF
	1079
	1080	! Set secondary weights
[11380]	1081	DO jm = 1, jpit
	1082	DO ji = jm, jpit
	1083	zwgt2(jm) = zwgt2(jm) + zwgt1(ji)
[4292]	1084	END DO
	1085	END DO
	1086
	1087	! Normalize weigths:
	1088	za1 = 1._wp / SUM(zwgt1(1:jpit))
	1089	za2 = 1._wp / SUM(zwgt2(1:jpit))
[11380]	1090	DO jm = 1, jpit
	1091	zwgt1(jm) = zwgt1(jm) * za1
	1092	zwgt2(jm) = zwgt2(jm) * za2
[4292]	1093	END DO
	1094	!
	1095	END SUBROUTINE ts_wgt
	1096
[6140]	1097
[508]	1098	SUBROUTINE ts_rst( kt, cdrw )
	1099	!!---------------------------------------------------------------------
	1100	!! * ROUTINE ts_rst *
	1101	!!
	1102	!! ** Purpose : Read or write time-splitting arrays in restart file
	1103	!!----------------------------------------------------------------------
[9528]	1104	INTEGER , INTENT(in) :: kt ! ocean time-step
	1105	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
[508]	1106	!!----------------------------------------------------------------------
	1107	!
[9506]	1108	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
	1109	! ! ---------------
[10256]	1110	IF( ln_rstart .AND. ln_bt_fw .AND. (neuler/=0) ) THEN !* Read the restart file
[9506]	1111	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios )
	1112	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios )
	1113	CALL iom_get( numror, jpdom_autoglo, 'un_bf' , un_bf (:,:), ldxios = lrxios )
	1114	CALL iom_get( numror, jpdom_autoglo, 'vn_bf' , vn_bf (:,:), ldxios = lrxios )
	1115	IF( .NOT.ln_bt_av ) THEN
	1116	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:), ldxios = lrxios )
	1117	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:), ldxios = lrxios )
	1118	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:), ldxios = lrxios )
	1119	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:), ldxios = lrxios )
	1120	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:), ldxios = lrxios )
	1121	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:), ldxios = lrxios )
	1122	ENDIF
[4486]	1123	#if defined key_agrif
[9506]	1124	! Read time integrated fluxes
	1125	IF ( .NOT.Agrif_Root() ) THEN
	1126	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lrxios )
	1127	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lrxios )
	1128	ENDIF
	1129	#endif
	1130	ELSE !* Start from rest
	1131	IF(lwp) WRITE(numout,*)
	1132	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
	1133	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
	1134	un_adv(:,:) = 0._wp ; vn_adv(:,:) = 0._wp ! used in the 1st interpol of agrif
	1135	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
	1136	#if defined key_agrif
	1137	IF ( .NOT.Agrif_Root() ) THEN
	1138	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
	1139	ENDIF
	1140	#endif
[4486]	1141	ENDIF
[9506]	1142	!
	1143	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
	1144	! ! -------------------
	1145	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
[9367]	1146	IF( lwxios ) CALL iom_swap( cwxios_context )
	1147	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios )
	1148	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios )
	1149	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios )
	1150	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios )
[4292]	1151	!
	1152	IF (.NOT.ln_bt_av) THEN
[9367]	1153	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:), ldxios = lwxios )
	1154	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:), ldxios = lwxios )
	1155	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:), ldxios = lwxios )
	1156	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:), ldxios = lwxios )
	1157	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:), ldxios = lwxios )
	1158	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:), ldxios = lwxios )
[4292]	1159	ENDIF
[4486]	1160	#if defined key_agrif
	1161	! Save time integrated fluxes
	1162	IF ( .NOT.Agrif_Root() ) THEN
[9367]	1163	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lwxios )
	1164	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lwxios )
[4486]	1165	ENDIF
	1166	#endif
[9367]	1167	IF( lwxios ) CALL iom_swap( cxios_context )
[4292]	1168	ENDIF
	1169	!
	1170	END SUBROUTINE ts_rst
[2528]	1171
[6140]	1172
	1173	SUBROUTINE dyn_spg_ts_init
[4292]	1174	!!---------------------------------------------------------------------
	1175	!! * ROUTINE dyn_spg_ts_init *
	1176	!!
	1177	!! ** Purpose : Set time splitting options
	1178	!!----------------------------------------------------------------------
[6140]	1179	INTEGER :: ji ,jj ! dummy loop indices
	1180	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
[9019]	1181	REAL(wp), DIMENSION(jpi,jpj) :: zcu
[11234]	1182	INTEGER :: inum
[4292]	1183	!!----------------------------------------------------------------------
[4370]	1184	!
[5930]	1185	! Max courant number for ext. grav. waves
[4370]	1186	!
[5930]	1187	DO jj = 1, jpj
	1188	DO ji =1, jpi
	1189	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
	1190	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
[7646]	1191	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
[4370]	1192	END DO
[5930]	1193	END DO
	1194	!
[5836]	1195	zcmax = MAXVAL( zcu(:,:) )
[10425]	1196	CALL mpp_max( 'dynspg_ts', zcmax )
[2528]	1197
[4370]	1198	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
[6140]	1199	IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
[4292]	1200
[5836]	1201	rdtbt = rdt / REAL( nn_baro , wp )
[4292]	1202	zcmax = zcmax * rdtbt
[9023]	1203	! Print results
[4292]	1204	IF(lwp) WRITE(numout,*)
[9169]	1205	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
	1206	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
[5930]	1207	IF( ln_bt_auto ) THEN
[9169]	1208	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_baro '
[4370]	1209	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
[4292]	1210	ELSE
[9169]	1211	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist nn_baro = ', nn_baro
[358]	1212	ENDIF
[4292]	1213
	1214	IF(ln_bt_av) THEN
[9169]	1215	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_baro time steps is on '
[4292]	1216	ELSE
[9169]	1217	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
[4292]	1218	ENDIF
[508]	1219	!
[4292]	1220	!
	1221	IF(ln_bt_fw) THEN
[4370]	1222	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
[4292]	1223	ELSE
[4370]	1224	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
[4292]	1225	ENDIF
	1226	!
[4486]	1227	#if defined key_agrif
	1228	! Restrict the use of Agrif to the forward case only
[9023]	1229	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
[4486]	1230	#endif
	1231	!
[4370]	1232	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
[4292]	1233	SELECT CASE ( nn_bt_flt )
[6140]	1234	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
	1235	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
	1236	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
[9169]	1237	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
[4292]	1238	END SELECT
	1239	!
[4370]	1240	IF(lwp) WRITE(numout,*) ' '
	1241	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
	1242	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
	1243	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
	1244	!
[9023]	1245	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
	1246	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
	1247	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
	1248	ENDIF
	1249	!
[6140]	1250	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
[4292]	1251	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
	1252	ENDIF
[6140]	1253	IF( zcmax>0.9_wp ) THEN
[4292]	1254	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
	1255	ENDIF
	1256	!
[9124]	1257	! ! Allocate time-splitting arrays
	1258	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
	1259	!
	1260	! ! read restart when needed
[9506]	1261	CALL ts_rst( nit000, 'READ' )
[9124]	1262	!
[9367]	1263	IF( lwxios ) THEN
	1264	! define variables in restart file when writing with XIOS
	1265	CALL iom_set_rstw_var_active('ub2_b')
	1266	CALL iom_set_rstw_var_active('vb2_b')
	1267	CALL iom_set_rstw_var_active('un_bf')
	1268	CALL iom_set_rstw_var_active('vn_bf')
	1269	!
	1270	IF (.NOT.ln_bt_av) THEN
	1271	CALL iom_set_rstw_var_active('sshbb_e')
	1272	CALL iom_set_rstw_var_active('ubb_e')
	1273	CALL iom_set_rstw_var_active('vbb_e')
	1274	CALL iom_set_rstw_var_active('sshb_e')
	1275	CALL iom_set_rstw_var_active('ub_e')
	1276	CALL iom_set_rstw_var_active('vb_e')
	1277	ENDIF
	1278	#if defined key_agrif
	1279	! Save time integrated fluxes
	1280	IF ( .NOT.Agrif_Root() ) THEN
	1281	CALL iom_set_rstw_var_active('ub2_i_b')
	1282	CALL iom_set_rstw_var_active('vb2_i_b')
	1283	ENDIF
	1284	#endif
	1285	ENDIF
	1286	!
[11380]	1287	! initialize extended scale factors
	1288	ht_0_xtd (1:jpi,1:jpj) = ht_0 (1:jpi,1:jpj)
	1289	hu_0_xtd (1:jpi,1:jpj) = hu_0 (1:jpi,1:jpj)
	1290	hv_0_xtd (1:jpi,1:jpj) = hv_0 (1:jpi,1:jpj)
	1291	r1_e1e2t_xtd(1:jpi,1:jpj) = r1_e1e2t(1:jpi,1:jpj)
	1292	r1_e1e2u_xtd(1:jpi,1:jpj) = r1_e1e2u(1:jpi,1:jpj)
	1293	r1_e1e2v_xtd(1:jpi,1:jpj) = r1_e1e2v(1:jpi,1:jpj)
	1294	e1e2t_xtd (1:jpi,1:jpj) = e1e2t (1:jpi,1:jpj)
	1295	ssmask_xtd (1:jpi,1:jpj) = ssmask (1:jpi,1:jpj)
	1296	ssumask_xtd (1:jpi,1:jpj) = ssumask (1:jpi,1:jpj)
	1297	ssvmask_xtd (1:jpi,1:jpj) = ssvmask (1:jpi,1:jpj)
	1298	e2u_xtd (1:jpi,1:jpj) = e2u (1:jpi,1:jpj)
	1299	e1v_xtd (1:jpi,1:jpj) = e1v (1:jpi,1:jpj)
	1300	r1_e1u_xtd (1:jpi,1:jpj) = r1_e1u (1:jpi,1:jpj)
	1301	r1_e2v_xtd (1:jpi,1:jpj) = r1_e2v (1:jpi,1:jpj)
	1302	!
	1303	CALL lbc_lnk_multi( 'dynspg_ts', ht_0_xtd , 'T', 1._wp, hu_0_xtd , 'U', -1._wp, hv_0_xtd , 'V', -1._wp &
	1304	& , r1_e1e2t_xtd, 'T', 1._wp, r1_e1e2u_xtd, 'U', -1._wp, r1_e1e2v_xtd, 'V', -1._wp &
	1305	& , ssmask_xtd , 'T', 1._wp, ssumask_xtd , 'U', -1._wp, ssvmask_xtd , 'V', -1._wp &
	1306	& , e1e2t_xtd , 'T', 1._wp, e2u_xtd , 'U', -1._wp, e1v_xtd , 'V', -1._wp &
	1307	& , r1_e1u_xtd , 'U', -1._wp, r1_e2v_xtd , 'V', -1._wp &
	1308	& , khlcom = nn_hls+nn_hlts )
	1309	IF( ln_dynvor_enT ) THEN
	1310	ff_t_xtd (1:jpi,1:jpj) = ff_t (1:jpi,1:jpj)
	1311	CALL lbc_lnk_multi( 'dynspg_ts', ff_t_xtd , 'F', -1._wp, khlcom = nn_hls+nn_hlts )
	1312	END IF
	1313
	1314	!
[4292]	1315	END SUBROUTINE dyn_spg_ts_init
[508]	1316
[11234]	1317
[11380]	1318	SUBROUTINE dyn_cor_2d_init( kdbi, kdei, kdbj, kdej )
[11234]	1319	!!---------------------------------------------------------------------
	1320	!! * ROUTINE dyn_cor_2d_init *
	1321	!!
	1322	!! ** Purpose : Set time splitting options
	1323	!! Set arrays to remove/compute coriolis trend.
	1324	!! Do it once during initialization if volume is fixed, else at each long time step.
	1325	!! Note that these arrays are also used during barotropic loop. These are however frozen
	1326	!! although they should be updated in the variable volume case. Not a big approximation.
	1327	!! To remove this approximation, copy lines below inside barotropic loop
	1328	!! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
	1329	!!
	1330	!! Compute zwz = f / ( height of the water colomn )
	1331	!!----------------------------------------------------------------------
[11380]	1332	INTEGER , INTENT(in ) :: kdbi, kdei, kdbj, kdej
[11234]	1333	INTEGER :: ji ,jj, jk ! dummy loop indices
	1334	REAL(wp) :: z1_ht
	1335	REAL(wp), DIMENSION(jpi,jpj) :: zhf
	1336	!!----------------------------------------------------------------------
	1337	!
	1338	SELECT CASE( nvor_scheme )
	1339	CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)
	1340	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
	1341	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
[11380]	1342	DO jj = 1, jpj-1
	1343	DO ji = 1, jpi-1
	1344	zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &
	1345	& ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp
[11234]	1346	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
	1347	END DO
	1348	END DO
	1349	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
[11380]	1350	DO jj = 1, jpj-1
	1351	DO ji = 1, jpi-1
	1352	zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) &
[11234]	1353	& + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) &
	1354	& / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) &
	1355	& + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) )
	1356	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
	1357	END DO
	1358	END DO
	1359	END SELECT
	1360	!
[11380]	1361	DO jj = 2, jpj-1
	1362	DO ji = 2, jpi-1
[11234]	1363	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	1364	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	1365	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	1366	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	1367	END DO
	1368	END DO
[11380]	1369	CALL lbc_lnk_multi( 'dynspg_ts', ftne, 'F', 1._wp, ftnw, 'F', 1._wp &
	1370	& , ftse, 'F', 1._wp, ftsw, 'F', 1._wp, khlcom = nn_hls+nn_hlts )
[11234]	1371	!
	1372	CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme)
	1373	DO jj = 2, jpj
	1374	DO ji = 2, jpi
	1375	z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )
	1376	ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht
	1377	ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht
	1378	ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
	1379	ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht
	1380	END DO
	1381	END DO
[11380]	1382	CALL lbc_lnk_multi( 'dynspg_ts', ftne, 'F', 1._wp, ftnw, 'F', 1._wp &
	1383	& , ftse, 'F', 1._wp, ftsw, 'F', 1._wp, khlcom = nn_hls+nn_hlts )
[11234]	1384	!
	1385	CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !
	1386	!
	1387	zwz(:,:) = 0._wp
	1388	zhf(:,:) = 0._wp
	1389
	1390	!!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
	1391	!!gm A priori a better value should be something like :
	1392	!!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
	1393	!!gm divided by the sum of the corresponding mask
	1394	!!gm
	1395	!!
	1396	IF( .NOT.ln_sco ) THEN
	1397
	1398	!!gm agree the JC comment : this should be done in a much clear way
	1399
	1400	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
	1401	! Set it to zero for the time being
	1402	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
	1403	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
	1404	! ENDIF
	1405	! zhf(:,:) = gdepw_0(:,:,jk+1)
	1406	!
	1407	ELSE
	1408	!
	1409	!zhf(:,:) = hbatf(:,:)
	1410	DO jj = 1, jpjm1
	1411	DO ji = 1, jpim1
	1412	zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &
	1413	& + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &
	1414	& / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &
	1415	& + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )
	1416	END DO
	1417	END DO
	1418	ENDIF
	1419	!
	1420	DO jj = 1, jpjm1
	1421	zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
	1422	END DO
	1423	!
	1424	DO jk = 1, jpkm1
	1425	DO jj = 1, jpjm1
	1426	zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
	1427	END DO
	1428	END DO
[11380]	1429	! JC: TBC. hf should be greater than 0
	1430	DO jj = 2, jpjm1
	1431	DO ji = 2, jpim1
[11234]	1432	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj)
	1433	END DO
	1434	END DO
[11380]	1435	zwz(2:jpim1,2:jpjm1) = ff_f(2:jpim1,2:jpjm1) * zwz(2:jpim1,2:jpjm1)
	1436	CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp, khlcom = nn_hls+nn_hlts )
[11234]	1437	END SELECT
	1438
	1439	END SUBROUTINE dyn_cor_2d_init
	1440
	1441
	1442
[11380]	1443	SUBROUTINE dyn_cor_2d( phgtt, phgtu, phgtv, pun, pvn, phU, phV, kdbi , kdei , kdbj , kdej , pu_trd, pv_trd &
	1444	& , kdbi2, kdei2, kdbj2, kdej2 )
[11234]	1445	!!---------------------------------------------------------------------
	1446	!! * ROUTINE dyn_cor_2d *
	1447	!!
	1448	!! ** Purpose : Compute u and v coriolis trends
[11380]	1449	!!
	1450	!! kdXX2 are useful in the initialisation where some arrays are not over the whole domain
	1451	!! and some are
[11234]	1452	!!----------------------------------------------------------------------
[11380]	1453	REAL(wp), DIMENSION(kdbi :kdei ,kdbj :kdej ), INTENT(in ) :: phgtt, phgtu, phgtv, pun, pvn ! height, speed
	1454	INTEGER , INTENT(in ) :: kdbi , kdei , kdbj , kdej ! arrays size
	1455	REAL(wp), DIMENSION(kdbi2:kdei2,kdbj2:kdej2), INTENT(in ) :: phU, phV ! flux
	1456	REAL(wp), DIMENSION(kdbi2:kdei2,kdbj2:kdej2), INTENT( out) :: pu_trd, pv_trd
	1457	INTEGER , INTENT(in ) :: kdbi2, kdei2, kdbj2, kdej2 ! arrays size
	1458	INTEGER :: ji, jj ! dummy loop indices
	1459	REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! local integer
[11234]	1460	!!----------------------------------------------------------------------
	1461	SELECT CASE( nvor_scheme )
	1462	CASE( np_ENT ) ! enstrophy conserving scheme (f-point)
[11380]	1463	DO jj = kdbj+1, kdej-1
	1464	DO ji = kdbi+1, kdei-1
	1465	z1_hu = ssumask_xtd(ji,jj) / ( phgtu(ji,jj) + 1._wp - ssumask_xtd(ji,jj) )
	1466	z1_hv = ssvmask_xtd(ji,jj) / ( phgtv(ji,jj) + 1._wp - ssvmask_xtd(ji,jj) )
	1467	pu_trd(ji,jj) = + r1_4 * r1_e1e2u_xtd(ji,jj) * z1_hu &
	1468	& * ( e1e2t_xtd(ji+1,jj)phgtt(ji+1,jj)ff_t_xtd(ji+1,jj) * ( pvn(ji+1,jj) + pvn(ji+1,jj-1) ) &
	1469	& + e1e2t_xtd(ji ,jj)phgtt(ji ,jj)ff_t_xtd(ji ,jj) * ( pvn(ji ,jj) + pvn(ji ,jj-1) ) )
[11234]	1470	!
[11380]	1471	pv_trd(ji,jj) = - r1_4 * r1_e1e2v_xtd(ji,jj) * z1_hv &
	1472	& * ( e1e2t_xtd(ji,jj+1)phgtt(ji,jj+1)ff_t_xtd(ji,jj+1) * ( pun(ji,jj+1) + pun(ji-1,jj+1) ) &
	1473	& + e1e2t_xtd(ji,jj )phgtt(ji,jj )ff_t_xtd(ji,jj ) * ( pun(ji,jj ) + pun(ji-1,jj ) ) )
[11234]	1474	END DO
[11380]	1475	END DO
[11234]	1476	!
	1477	CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX
[11380]	1478	DO jj = kdbj+1, kdej-1
	1479	DO ji = kdbi+1, kdei-1
	1480	zy1 = ( phV(ji,jj-1) + phV(ji+1,jj-1) ) * r1_e1u_xtd(ji,jj)
	1481	zy2 = ( phV(ji,jj ) + phV(ji+1,jj ) ) * r1_e1u_xtd(ji,jj)
	1482	zx1 = ( phU(ji-1,jj) + phU(ji-1,jj+1) ) * r1_e2v_xtd(ji,jj)
	1483	zx2 = ( phU(ji ,jj) + phU(ji ,jj+1) ) * r1_e2v_xtd(ji,jj)
[11234]	1484	! energy conserving formulation for planetary vorticity term
[11380]	1485	pu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	1486	pv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
[11234]	1487	END DO
	1488	END DO
	1489	!
	1490	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
[11380]	1491	DO jj = kdbj+1, kdej-1
	1492	DO ji = kdbi+1, kdei-1
	1493	zy1 = r1_8 * ( phV(ji ,jj-1) + phV(ji+1,jj-1) &
	1494	& + phV(ji ,jj ) + phV(ji+1,jj ) ) * r1_e1u_xtd(ji,jj)
	1495	zx1 = - r1_8 * ( phU(ji-1,jj ) + phU(ji-1,jj+1) &
	1496	& + phU(ji ,jj ) + phU(ji ,jj+1) ) * r1_e2v_xtd(ji,jj)
	1497	pu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	1498	pv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
[11234]	1499	END DO
	1500	END DO
	1501	!
[11380]	1502	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
	1503	DO jj = kdbj+1, kdej-1
	1504	DO ji = kdbi+1, kdei-1
	1505	pu_trd(ji,jj) = + r1_12 * r1_e1u_xtd(ji,jj) * ( ftne(ji,jj ) * phV(ji ,jj ) &
	1506	& + ftnw(ji+1,jj) * phV(ji+1,jj ) &
	1507	& + ftse(ji,jj ) * phV(ji ,jj-1) &
	1508	& + ftsw(ji+1,jj) * phV(ji+1,jj-1) )
	1509	pv_trd(ji,jj) = - r1_12 * r1_e2v_xtd(ji,jj) * ( ftsw(ji,jj+1) * phU(ji-1,jj+1) &
	1510	& + ftse(ji,jj+1) * phU(ji ,jj+1) &
	1511	& + ftnw(ji,jj ) * phU(ji-1,jj ) &
	1512	& + ftne(ji,jj ) * phU(ji ,jj ) )
[11234]	1513	END DO
	1514	END DO
	1515	!
	1516	END SELECT
	1517	!
[11380]	1518	END SUBROUTINE dyn_cor_2d
[11234]	1519
	1520
	1521	SUBROUTINE wad_tmsk( pssh, ptmsk )
	1522	!!----------------------------------------------------------------------
	1523	!! * ROUTINE wad_lmt *
	1524	!!
	1525	!! ** Purpose : set wetting & drying mask at tracer points
	1526	!! for the current barotropic sub-step
	1527	!!
	1528	!! ** Method : ???
	1529	!!
	1530	!! ** Action : ptmsk : wetting & drying t-mask
	1531	!!----------------------------------------------------------------------
	1532	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh !
	1533	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk !
	1534	!
	1535	INTEGER :: ji, jj ! dummy loop indices
	1536	!!----------------------------------------------------------------------
	1537	!
	1538	IF( ln_wd_dl_rmp ) THEN
	1539	DO jj = 1, jpj
	1540	DO ji = 1, jpi
	1541	IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
	1542	! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
	1543	ptmsk(ji,jj) = 1._wp
	1544	ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
	1545	ptmsk(ji,jj) = TANH( 50._wp( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1) )
	1546	ELSE
	1547	ptmsk(ji,jj) = 0._wp
	1548	ENDIF
	1549	END DO
	1550	END DO
	1551	ELSE
	1552	DO jj = 1, jpj
	1553	DO ji = 1, jpi
	1554	IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp
	1555	ELSE ; ptmsk(ji,jj) = 0._wp
	1556	ENDIF
	1557	END DO
	1558	END DO
	1559	ENDIF
	1560	!
	1561	END SUBROUTINE wad_tmsk
	1562
	1563
	1564	SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
	1565	!!----------------------------------------------------------------------
	1566	!! * ROUTINE wad_lmt *
	1567	!!
	1568	!! ** Purpose : set wetting & drying mask at tracer points
	1569	!! for the current barotropic sub-step
	1570	!!
	1571	!! ** Method : ???
	1572	!!
	1573	!! ** Action : ptmsk : wetting & drying t-mask
	1574	!!----------------------------------------------------------------------
	1575	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
	1576	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
	1577	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
	1578	!
	1579	INTEGER :: ji, jj ! dummy loop indices
	1580	!!----------------------------------------------------------------------
	1581	!
	1582	DO jj = 1, jpj
	1583	DO ji = 1, jpim1 ! not jpi-column
	1584	IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj)
	1585	ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj)
	1586	ENDIF
	1587	phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
	1588	pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
	1589	END DO
	1590	END DO
	1591	!
	1592	DO jj = 1, jpjm1 ! not jpj-row
	1593	DO ji = 1, jpi
	1594	IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj )
	1595	ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1)
	1596	ENDIF
	1597	phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
	1598	pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
	1599	END DO
	1600	END DO
	1601	!
	1602	END SUBROUTINE wad_Umsk
	1603
	1604
	1605	SUBROUTINE wad_spg( sshn, zcpx, zcpy )
	1606	!!---------------------------------------------------------------------
	1607	!! * ROUTINE wad_sp *
	1608	!!
	1609	!! ** Purpose :
	1610	!!----------------------------------------------------------------------
	1611	INTEGER :: ji ,jj ! dummy loop indices
	1612	LOGICAL :: ll_tmp1, ll_tmp2
	1613	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: sshn
	1614	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
	1615	!!----------------------------------------------------------------------
	1616	DO jj = 2, jpjm1
	1617	DO ji = 2, jpim1
	1618	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > &
	1619	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
	1620	& MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) &
	1621	& > rn_wdmin1 + rn_wdmin2
	1622	ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( &
	1623	& MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > &
	1624	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
	1625	IF(ll_tmp1) THEN
	1626	zcpx(ji,jj) = 1.0_wp
	1627	ELSEIF(ll_tmp2) THEN
	1628	! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here
	1629	zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
	1630	& / (sshn(ji+1,jj) - sshn(ji ,jj)) )
	1631	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
	1632	ELSE
	1633	zcpx(ji,jj) = 0._wp
	1634	ENDIF
	1635	!
	1636	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > &
	1637	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
	1638	& MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) &
	1639	& > rn_wdmin1 + rn_wdmin2
	1640	ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( &
	1641	& MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > &
	1642	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
	1643
	1644	IF(ll_tmp1) THEN
	1645	zcpy(ji,jj) = 1.0_wp
	1646	ELSE IF(ll_tmp2) THEN
	1647	! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here
	1648	zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
	1649	& / (sshn(ji,jj+1) - sshn(ji,jj )) )
	1650	zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) )
	1651	ELSE
	1652	zcpy(ji,jj) = 0._wp
	1653	ENDIF
	1654	END DO
	1655	END DO
	1656
	1657	END SUBROUTINE wad_spg
	1658
	1659
	1660
[11380]	1661	SUBROUTINE drg_init( kdbi, kdei, kdbj, kdej, pu_frc, pv_frc, pCdU_u, pCdU_v )
[11234]	1662	!!----------------------------------------------------------------------
[11380]	1663	!! * ROUTINE drg_init *
[11234]	1664	!!
	1665	!! ** Purpose : - add the baroclinic top/bottom drag contribution to
	1666	!! the baroclinic part of the barotropic RHS
	1667	!! - compute the barotropic drag coefficients
	1668	!!
	1669	!! ** Method : computation done over the INNER domain only
	1670	!!----------------------------------------------------------------------
[11380]	1671	INTEGER , INTENT(in ) :: kdbi, kdei, kdbj, kdej ! arrays size
	1672	REAL(wp), DIMENSION(kdbi:kdei,kdbj:kdej), INTENT(inout) :: pu_frc, pv_frc ! baroclinic part of the barotropic RHS
	1673	REAL(wp), DIMENSION(kdbi:kdei,kdbj:kdej), INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients
[11234]	1674	!
	1675	INTEGER :: ji, jj ! dummy loop indices
	1676	INTEGER :: ikbu, ikbv, iktu, iktv
	1677	REAL(wp) :: zztmp
	1678	REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i
	1679	!!----------------------------------------------------------------------
	1680	!
	1681	! !== Set the barotropic drag coef. ==!
	1682	!
	1683	IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)
	1684
	1685	DO jj = 2, jpjm1
	1686	DO ji = 2, jpim1 ! INNER domain
	1687	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
	1688	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
	1689	END DO
	1690	END DO
	1691	ELSE ! bottom friction only
	1692	DO jj = 2, jpjm1
	1693	DO ji = 2, jpim1 ! INNER domain
	1694	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
	1695	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
	1696	END DO
	1697	END DO
	1698	ENDIF
	1699	!
	1700	! !== BOTTOM stress contribution from baroclinic velocities ==!
	1701	!
	1702	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
	1703
	1704	DO jj = 2, jpjm1
	1705	DO ji = 2, jpim1 ! INNER domain
	1706	ikbu = mbku(ji,jj)
	1707	ikbv = mbkv(ji,jj)
	1708	zu_i(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj)
	1709	zv_i(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
	1710	END DO
	1711	END DO
	1712	ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
	1713
	1714	DO jj = 2, jpjm1
	1715	DO ji = 2, jpim1 ! INNER domain
	1716	ikbu = mbku(ji,jj)
	1717	ikbv = mbkv(ji,jj)
	1718	zu_i(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj)
	1719	zv_i(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
	1720	END DO
	1721	END DO
	1722	ENDIF
	1723	!
	1724	IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please !
	1725	zztmp = -1._wp / rdtbt
	1726	DO jj = 2, jpjm1
	1727	DO ji = 2, jpim1 ! INNER domain
[11380]	1728	pu_frc(ji,jj) = pu_frc(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
[11234]	1729	& r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp )
[11380]	1730	pv_frc(ji,jj) = pv_frc(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &
[11234]	1731	& r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
	1732	END DO
	1733	END DO
	1734	ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
	1735
	1736	DO jj = 2, jpjm1
	1737	DO ji = 2, jpim1 ! INNER domain
[11380]	1738	pu_frc(ji,jj) = pu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) zu_i(ji,jj)
	1739	pv_frc(ji,jj) = pv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) zv_i(ji,jj)
[11234]	1740	END DO
	1741	END DO
	1742	END IF
	1743	!
	1744	! !== TOP stress contribution from baroclinic velocities ==! (no W/D case)
	1745	!
	1746	IF( ln_isfcav ) THEN
	1747	!
	1748	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
	1749
	1750	DO jj = 2, jpjm1
	1751	DO ji = 2, jpim1 ! INNER domain
	1752	iktu = miku(ji,jj)
	1753	iktv = mikv(ji,jj)
	1754	zu_i(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj)
	1755	zv_i(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
	1756	END DO
	1757	END DO
	1758	ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
	1759
	1760	DO jj = 2, jpjm1
	1761	DO ji = 2, jpim1 ! INNER domain
	1762	iktu = miku(ji,jj)
	1763	iktv = mikv(ji,jj)
	1764	zu_i(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj)
	1765	zv_i(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
	1766	END DO
	1767	END DO
	1768	ENDIF
	1769	!
	1770	! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
	1771
	1772	DO jj = 2, jpjm1
	1773	DO ji = 2, jpim1 ! INNER domain
[11380]	1774	pu_frc(ji,jj) = pu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) zu_i(ji,jj)
	1775	pv_frc(ji,jj) = pv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) zv_i(ji,jj)
[11234]	1776	END DO
	1777	END DO
	1778	!
	1779	ENDIF
	1780	!
[11380]	1781	END SUBROUTINE drg_init
[11234]	1782
[11380]	1783
	1784	SUBROUTINE ts_bck_interp( km, ld_init, & ! <== in
[11234]	1785	& za0, za1, za2, za3 ) ! ==> out
	1786	!!----------------------------------------------------------------------
[11380]	1787	INTEGER ,INTENT(in ) :: km ! index of sub time step
	1788	LOGICAL ,INTENT(in ) :: ld_init !
[11234]	1789	REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient
	1790	!
	1791	REAL(wp) :: zepsilon, zgamma ! - -
	1792	!!----------------------------------------------------------------------
	1793	! ! set Half-step back interpolation coefficient
[11380]	1794	IF ( km==1 .AND. ld_init ) THEN !* Forward-backward
[11234]	1795	za0 = 1._wp
	1796	za1 = 0._wp
	1797	za2 = 0._wp
	1798	za3 = 0._wp
[11380]	1799	ELSEIF( km==2 .AND. ld_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
[11234]	1800	za0 = 1.0833333333333_wp ! za0 = 1-gam-eps
	1801	za1 =-0.1666666666666_wp ! za1 = gam
	1802	za2 = 0.0833333333333_wp ! za2 = eps
	1803	za3 = 0._wp
	1804	ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
	1805	IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion
	1806	za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps
	1807	za1 = 0.285_wp ! za1 = 1/2 - 2gam - 3eps
	1808	za2 = 0.088_wp ! za2 = gam
	1809	za3 = 0.013_wp ! za3 = eps
	1810	ELSE ! no time diffusion
	1811	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
	1812	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
	1813	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
	1814	za1 = 1._wp - za0 - zgamma - zepsilon
	1815	za2 = zgamma
	1816	za3 = zepsilon
	1817	ENDIF
	1818	ENDIF
	1819	END SUBROUTINE ts_bck_interp
	1820
	1821
[358]	1822	!!======================================================================
	1823	END MODULE dynspg_ts

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

Download in other formats: