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/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN/dynspg_ts.F90 @ 13229

Last change on this file since 13229 was 13229, checked in by francesca, 4 years ago
dev_r12558_HPC-08_epico_Extra_Halo: merge with trunk@13218, see #2366
Property svn:keywords set to `Id`
File size: 78.7 KB

Rev	Line
[358]	1	MODULE dynspg_ts
[9023]	2
[12377]	3	!! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out !
[9023]	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
[12377]	33	USE isf_oce ! ice shelf variable (fwfisf)
[9019]	34	USE zdf_oce ! vertical physics: variables
	35	USE zdfdrg ! vertical physics: top/bottom drag coef.
[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
[12377]	46	USE tide_mod !
[7646]	47	USE sbcwave ! surface wave
[9019]	48	#if defined key_agrif
[9570]	49	USE agrif_oce_interp ! agrif
[9124]	50	USE agrif_oce
[9019]	51	#endif
	52	#if defined key_asminc
	53	USE asminc ! Assimilation increment
	54	#endif
[6140]	55	!
	56	USE in_out_manager ! I/O manager
[358]	57	USE lib_mpp ! distributed memory computing library
	58	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	59	USE prtctl ! Print control
[2715]	60	USE iom ! IOM library
[4292]	61	USE restart ! only for lrst_oce
[358]	62
[11536]	63	USE iom ! to remove
	64
[358]	65	IMPLICIT NONE
	66	PRIVATE
	67
[9124]	68	PUBLIC dyn_spg_ts ! called by dyn_spg
	69	PUBLIC dyn_spg_ts_init ! - - dyn_spg_init
[358]	70
[9019]	71	!! Time filtered arrays at baroclinic time step:
	72	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step
[9124]	73	!
[12489]	74	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_e <= 2.5 nn_e
	75	REAL(wp),SAVE :: rDt_e ! Barotropic time step
[9019]	76	!
	77	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields
[9124]	78	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points
	79	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
	80	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
[4292]	81
[9043]	82	REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios
	83	REAL(wp) :: r1_8 = 0.125_wp !
	84	REAL(wp) :: r1_4 = 0.25_wp !
	85	REAL(wp) :: r1_2 = 0.5_wp !
[508]	86
[358]	87	!! * Substitutions
[12377]	88	# include "do_loop_substitute.h90"
[2715]	89	!!----------------------------------------------------------------------
[9598]	90	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5217]	91	!! $Id$
[10068]	92	!! Software governed by the CeCILL license (see ./LICENSE)
[2715]	93	!!----------------------------------------------------------------------
[358]	94	CONTAINS
	95
[2715]	96	INTEGER FUNCTION dyn_spg_ts_alloc()
	97	!!----------------------------------------------------------------------
	98	!! * routine dyn_spg_ts_alloc *
	99	!!----------------------------------------------------------------------
[6140]	100	INTEGER :: ierr(3)
[4292]	101	!!----------------------------------------------------------------------
	102	ierr(:) = 0
[6140]	103	!
[12489]	104	ALLOCATE( wgtbtp1(3nn_e), wgtbtp2(3nn_e), zwz(jpi,jpj), STAT=ierr(1) )
[9528]	105	IF( ln_dynvor_een .OR. ln_dynvor_eeT ) &
[11536]	106	& ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2) )
[6140]	107	!
	108	ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) )
	109	!
	110	dyn_spg_ts_alloc = MAXVAL( ierr(:) )
	111	!
[10425]	112	CALL mpp_sum( 'dynspg_ts', dyn_spg_ts_alloc )
	113	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_stop( 'STOP', 'dyn_spg_ts_alloc: failed to allocate arrays' )
[2715]	114	!
	115	END FUNCTION dyn_spg_ts_alloc
	116
[5836]	117
[12377]	118	SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )
[358]	119	!!----------------------------------------------------------------------
	120	!!
[6140]	121	!! ** Purpose : - Compute the now trend due to the explicit time stepping
	122	!! of the quasi-linear barotropic system, and add it to the
	123	!! general momentum trend.
[358]	124	!!
[6140]	125	!! ** Method : - split-explicit schem (time splitting) :
[4374]	126	!! Barotropic variables are advanced from internal time steps
	127	!! "n" to "n+1" if ln_bt_fw=T
	128	!! or from
	129	!! "n-1" to "n+1" if ln_bt_fw=F
	130	!! thanks to a generalized forward-backward time stepping (see ref. below).
[358]	131	!!
[4374]	132	!! ** Action :
[12377]	133	!! -Update the filtered free surface at step "n+1" : pssh(:,:,Kaa)
	134	!! -Update filtered barotropic velocities at step "n+1" : puu_b(:,:,:,Kaa), vv_b(:,:,:,Kaa)
[9023]	135	!! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv
[4374]	136	!! These are used to advect tracers and are compliant with discrete
	137	!! continuity equation taken at the baroclinic time steps. This
	138	!! ensures tracers conservation.
[12377]	139	!! - (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) momentum trend updated with barotropic component.
[358]	140	!!
[6140]	141	!! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
[358]	142	!!---------------------------------------------------------------------
[12377]	143	INTEGER , INTENT( in ) :: kt ! ocean time-step index
	144	INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices
	145	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
	146	REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels
[2715]	147	!
[9554]	148	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[9019]	149	LOGICAL :: ll_fw_start ! =T : forward integration
[9554]	150	LOGICAL :: ll_init ! =T : special startup of 2d equations
[11536]	151	INTEGER :: noffset ! local integers : time offset for bdy update
[12489]	152	REAL(wp) :: r1_Dt_b, z1_hu, z1_hv ! local scalars
[9528]	153	REAL(wp) :: za0, za1, za2, za3 ! - -
[12206]	154	REAL(wp) :: zztmp, zldg ! - -
[11536]	155	REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - -
	156	REAL(wp) :: zun_save, zvn_save ! - -
	157	REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc
	158	REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg
	159	REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhtp2_e
	160	REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e
[9019]	161	REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points
[11536]	162	REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes
[3294]	163	!
[9023]	164	REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True.
	165
	166	INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point
	167
	168	REAL(wp) :: zepsilon, zgamma ! - -
[9019]	169	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef.
[9023]	170	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
	171	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask
[12377]	172	REAL(wp) :: zt0substep ! Time of day at the beginning of the time substep
[358]	173	!!----------------------------------------------------------------------
[3294]	174	!
[9023]	175	IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) )
	176	! !* Allocate temporary arrays
	177	IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
[3294]	178	!
[9023]	179	zwdramp = r_rn_wdmin1 ! simplest ramp
	180	! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
[11536]	181	! ! inverse of baroclinic time step
[12489]	182	r1_Dt_b = 1._wp / rDt
[4292]	183	!
[9023]	184	ll_init = ln_bt_av ! if no time averaging, then no specific restart
[4292]	185	ll_fw_start = .FALSE.
[9023]	186	! ! time offset in steps for bdy data update
[12489]	187	IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_e
[6140]	188	ELSE ; noffset = 0
	189	ENDIF
[4292]	190	!
[9023]	191	IF( kt == nit000 ) THEN !* initialisation
[508]	192	!
[358]	193	IF(lwp) WRITE(numout,*)
	194	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
	195	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
[4354]	196	IF(lwp) WRITE(numout,*)
[1502]	197	!
[12489]	198	IF( l_1st_euler ) ll_init=.TRUE.
[1502]	199	!
[12489]	200	IF( ln_bt_fw .OR. l_1st_euler ) THEN
[6140]	201	ll_fw_start =.TRUE.
	202	noffset = 0
[4292]	203	ELSE
[6140]	204	ll_fw_start =.FALSE.
[4292]	205	ENDIF
[11536]	206	! ! Set averaging weights and cycle length:
[6140]	207	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
[4292]	208	!
[12489]	209	ELSEIF( kt == nit000 + 1 ) THEN !* initialisation 2nd time-step
	210	!
	211	IF( .NOT.ln_bt_fw ) THEN
	212	! If we did an Euler timestep on the first timestep we need to reset ll_fw_start
	213	! and the averaging weights. We don't have an easy way of telling whether we did
	214	! an Euler timestep on the first timestep (because l_1st_euler is reset to .false.
	215	! at the end of the first timestep) so just do this in all cases.
	216	ll_fw_start = .FALSE.
	217	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
	218	ENDIF
	219	!
[4292]	220	ENDIF
	221	!
[358]	222	! -----------------------------------------------------------------------------
	223	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
	224	! -----------------------------------------------------------------------------
[1502]	225	!
[4292]	226	!
[11536]	227	! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends)
	228	! ! --------------------------- !
[12377]	229	zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm)
	230	zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm)
[1502]	231	!
[4292]	232	!
[12377]	233	! != U(Krhs) => baroclinic trend =! (remove its vertical mean)
	234	DO jk = 1, jpkm1 ! ----------------------------- !
	235	uu(:,:,jk,Krhs) = ( uu(:,:,jk,Krhs) - zu_frc(:,:) ) * umask(:,:,jk)
	236	vv(:,:,jk,Krhs) = ( vv(:,:,jk,Krhs) - zv_frc(:,:) ) * vmask(:,:,jk)
[1502]	237	END DO
[7646]	238
	239	!!gm Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
	240	!!gm Is it correct to do so ? I think so...
	241
[11536]	242	! != remove 2D Coriolis and pressure gradient trends =!
	243	! ! ------------------------------------------------- !
[9528]	244	!
[12377]	245	IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init( Kmm ) ! Set zwz, the barotropic Coriolis force coefficient
[11536]	246	! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes
[1502]	247	!
[11536]	248	! !* 2D Coriolis trends
[12377]	249	zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes
	250	zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column
[11536]	251	!
[12377]	252	CALL dyn_cor_2d( hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in
	253	& zu_trd, zv_trd ) ! ==>> out
[11536]	254	!
	255	IF( .NOT.ln_linssh ) THEN !* surface pressure gradient (variable volume only)
[508]	256	!
[11536]	257	IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg
[12377]	258	CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy
	259	DO_2D_00_00
	260	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) &
	261	& * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth
	262	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) &
	263	& * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth
	264	END_2D
[11536]	265	ELSE ! now suface pressure gradient
[12377]	266	DO_2D_00_00
	267	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e1u(ji,jj)
	268	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e2v(ji,jj)
	269	END_2D
[9019]	270	ENDIF
	271	!
[1502]	272	ENDIF
[9019]	273	!
[12377]	274	DO_2D_00_00
	275	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
	276	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
	277	END_2D
[4292]	278	!
[11536]	279	! != Add bottom stress contribution from baroclinic velocities =!
	280	! ! ----------------------------------------------------------- !
[12377]	281	CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients
[11536]	282	! != Add atmospheric pressure forcing =!
	283	! ! ---------------------------------- !
	284	IF( ln_apr_dyn ) THEN
	285	IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
[12377]	286	DO_2D_00_00
	287	zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
	288	zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
	289	END_2D
[11536]	290	ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
	291	zztmp = grav * r1_2
[12377]	292	DO_2D_00_00
	293	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
	294	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
	295	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
	296	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
	297	END_2D
	298	ENDIF
[9019]	299	ENDIF
[11536]	300	!
	301	! != Add atmospheric pressure forcing =!
	302	! ! ---------------------------------- !
[9019]	303	IF( ln_bt_fw ) THEN ! Add wind forcing
[12377]	304	DO_2D_00_00
[12489]	305	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)
	306	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)
[12377]	307	END_2D
[2724]	308	ELSE
[12489]	309	zztmp = r1_rho0 * r1_2
[12377]	310	DO_2D_00_00
	311	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm)
	312	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm)
	313	END_2D
[4292]	314	ENDIF
	315	!
[11536]	316	! !----------------!
	317	! !== sssh_frc ==! Right-Hand-Side of the barotropic ssh equation (over the FULL domain)
	318	! !----------------!
	319	! != Net water flux forcing applied to a water column =!
	320	! ! --------------------------------------------------- !
	321	IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
[12489]	322	zssh_frc(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) )
[11536]	323	ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
[12489]	324	zztmp = r1_rho0 * r1_2
[12377]	325	zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) &
	326	& - rnf(:,:) - rnf_b(:,:) &
	327	& + fwfisf_cav(:,:) + fwfisf_cav_b(:,:) &
	328	& + fwfisf_par(:,:) + fwfisf_par_b(:,:) )
[4292]	329	ENDIF
[11536]	330	! != Add Stokes drift divergence =! (if exist)
	331	IF( ln_sdw ) THEN ! ----------------------------- !
[7753]	332	zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
[7646]	333	ENDIF
	334	!
[12377]	335	! ! ice sheet coupling
	336	IF ( ln_isf .AND. ln_isfcpl ) THEN
	337	!
	338	! ice sheet coupling
	339	IF( ln_rstart .AND. kt == nit000 ) THEN
	340	zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_ssh(:,:)
	341	END IF
	342	!
	343	! conservation option
	344	IF( ln_isfcpl_cons ) THEN
	345	zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_cons_ssh(:,:)
	346	END IF
	347	!
	348	END IF
	349	!
[4292]	350	#if defined key_asminc
[11536]	351	! != Add the IAU weighted SSH increment =!
	352	! ! ------------------------------------ !
[4292]	353	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
[7753]	354	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
[4292]	355	ENDIF
	356	#endif
[11536]	357	! != Fill boundary data arrays for AGRIF
[5656]	358	! ! ------------------------------------
[4486]	359	#if defined key_agrif
	360	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
	361	#endif
[4292]	362	!
[358]	363	! -----------------------------------------------------------------------
[4292]	364	! Phase 2 : Integration of the barotropic equations
[358]	365	! -----------------------------------------------------------------------
[1502]	366	!
	367	! ! ==================== !
	368	! ! Initialisations !
[4292]	369	! ! ==================== !
[4370]	370	! Initialize barotropic variables:
[4770]	371	IF( ll_init )THEN
[7753]	372	sshbb_e(:,:) = 0._wp
	373	ubb_e (:,:) = 0._wp
	374	vbb_e (:,:) = 0._wp
	375	sshb_e (:,:) = 0._wp
	376	ub_e (:,:) = 0._wp
	377	vb_e (:,:) = 0._wp
[4700]	378	ENDIF
	379	!
[11536]	380	IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0)
[12377]	381	zhup2_e(:,:) = hu(:,:,Kmm)
	382	zhvp2_e(:,:) = hv(:,:,Kmm)
	383	zhtp2_e(:,:) = ht(:,:)
[11536]	384	ENDIF
	385	!
[4370]	386	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
[12377]	387	sshn_e(:,:) = pssh(:,:,Kmm)
	388	un_e (:,:) = puu_b(:,:,Kmm)
	389	vn_e (:,:) = pvv_b(:,:,Kmm)
[7753]	390	!
[12377]	391	hu_e (:,:) = hu(:,:,Kmm)
	392	hv_e (:,:) = hv(:,:,Kmm)
	393	hur_e (:,:) = r1_hu(:,:,Kmm)
	394	hvr_e (:,:) = r1_hv(:,:,Kmm)
[4370]	395	ELSE ! CENTRED integration: start from BEFORE fields
[12377]	396	sshn_e(:,:) = pssh(:,:,Kbb)
	397	un_e (:,:) = puu_b(:,:,Kbb)
	398	vn_e (:,:) = pvv_b(:,:,Kbb)
[7753]	399	!
[12377]	400	hu_e (:,:) = hu(:,:,Kbb)
	401	hv_e (:,:) = hv(:,:,Kbb)
	402	hur_e (:,:) = r1_hu(:,:,Kbb)
	403	hvr_e (:,:) = r1_hv(:,:,Kbb)
[4292]	404	ENDIF
	405	!
	406	! Initialize sums:
[12377]	407	puu_b (:,:,Kaa) = 0._wp ! After barotropic velocities (or transport if flux form)
	408	pvv_b (:,:,Kaa) = 0._wp
	409	pssh (:,:,Kaa) = 0._wp ! Sum for after averaged sea level
[7753]	410	un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop
	411	vn_adv(:,:) = 0._wp
[9528]	412	!
	413	IF( ln_wd_dl ) THEN
[9023]	414	zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary)
	415	zvwdmask(:,:) = 0._wp !
[9528]	416	zuwdav2 (:,:) = 0._wp
	417	zvwdav2 (:,:) = 0._wp
[9023]	418	END IF
	419
[9528]	420	! ! ==================== !
[4292]	421	DO jn = 1, icycle ! sub-time-step loop !
[1502]	422	! ! ==================== !
[10425]	423	!
	424	l_full_nf_update = jn == icycle ! false: disable full North fold update (performances) for jn = 1 to icycle-1
[4292]	425	!
[11536]	426	! !== Update the forcing ==! (BDY and tides)
	427	!
[12377]	428	IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, pt_offset= REAL(noffset+1,wp) )
	429	! Update tide potential at the beginning of current time substep
	430	IF( ln_tide_pot .AND. ln_tide ) THEN
[12489]	431	zt0substep = REAL(nsec_day, wp) - 0.5_wprn_Dt + (jn + noffset - 1) rn_Dt / REAL(nn_e, wp)
[12377]	432	CALL upd_tide(zt0substep, Kmm)
	433	END IF
[11536]	434	!
	435	! !== extrapolation at mid-step ==! (jn+1/2)
	436	!
	437	! !* Set extrapolation coefficients for predictor step:
[4292]	438	IF ((jn<3).AND.ll_init) THEN ! Forward
	439	za1 = 1._wp
	440	za2 = 0._wp
	441	za3 = 0._wp
	442	ELSE ! AB3-AM4 Coefficients: bet=0.281105
	443	za1 = 1.781105_wp ! za1 = 3/2 + bet
	444	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
	445	za3 = 0.281105_wp ! za3 = bet
	446	ENDIF
[11536]	447	!
	448	! !* Extrapolate barotropic velocities at mid-step (jn+1/2)
	449	!-- m+1/2 m m-1 m-2 --!
	450	!-- u = (3/2+beta) u -(1/2+2beta) u + beta u --!
	451	!-------------------------------------------------------------------------!
[7753]	452	ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
	453	va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
[4292]	454
[6140]	455	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
[4292]	456	! ! ------------------
	457	! Extrapolate Sea Level at step jit+0.5:
[11536]	458	!-- m+1/2 m m-1 m-2 --!
	459	!-- ssh = (3/2+beta) ssh -(1/2+2beta) ssh + beta ssh --!
	460	!--------------------------------------------------------------------------------!
[7753]	461	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
[9023]	462
[11536]	463	! set wetting & drying mask at tracer points for this barotropic mid-step
	464	IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask )
[9528]	465	!
[11536]	466	! ! ocean t-depth at mid-step
	467	zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:)
	468	!
	469	! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk)
[12377]	470	DO_2D_11_10
	471	zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) &
	472	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
	473	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
	474	END_2D
	475	DO_2D_10_11
	476	zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) &
	477	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
	478	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
	479	END_2D
[4292]	480	!
	481	ENDIF
	482	!
[11536]	483	! !== after SSH ==! (jn+1)
	484	!
	485	! ! update (ua_e,va_e) to enforce volume conservation at open boundaries
	486	! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
[10481]	487	IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
[12377]	488	!
[11536]	489	! ! resulting flux at mid-step (not over the full domain)
	490	zhU(1:jpim1,1:jpj ) = e2u(1:jpim1,1:jpj ) * ua_e(1:jpim1,1:jpj ) * zhup2_e(1:jpim1,1:jpj ) ! not jpi-column
	491	zhV(1:jpi ,1:jpjm1) = e1v(1:jpi ,1:jpjm1) * va_e(1:jpi ,1:jpjm1) * zhvp2_e(1:jpi ,1:jpjm1) ! not jpj-row
[4486]	492	!
	493	#if defined key_agrif
[6140]	494	! Set fluxes during predictor step to ensure volume conservation
[12377]	495	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV )
[4486]	496	#endif
[12489]	497	IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rDt_e) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV
[9023]	498
[11536]	499	IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where
	500	! ! the direction of the flow is from dry cells
	501	CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V
[9528]	502	!
	503	ENDIF
[11536]	504	!
	505	!
	506	! Compute Sea Level at step jit+1
	507	!-- m+1 m m+1/2 --!
	508	!-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --!
	509	!-------------------------------------------------------------------------!
[12377]	510	DO_2D_00_00
	511	zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj)
[12489]	512	ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)
[12377]	513	END_2D
[11536]	514	!
[13229]	515	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
	516	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
	517	#if defined key_agrif
	518	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
	519	#endif
	520
[11536]	521	CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )
	522	!
	523	! ! Sum over sub-time-steps to compute advective velocities
	524	za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5
	525	un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
	526	vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
	527	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
	528	IF ( ln_wd_dl_bc ) THEN
	529	zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column
	530	zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row
	531	END IF
	532	!
[4292]	533	!
	534	! Sea Surface Height at u-,v-points (vvl case only)
[6140]	535	IF( .NOT.ln_linssh ) THEN
[12377]	536	DO_2D_00_00
	537	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
	538	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
	539	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
	540	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
	541	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
	542	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
	543	END_2D
[4292]	544	ENDIF
[11536]	545	!
	546	! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
	547	!-- m+1/2 m+1 m m-1 m-2 --!
	548	!-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --!
	549	!------------------------------------------------------------------------------------------!
	550	CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation
[9528]	551	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
	552	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
[1502]	553	!
[11536]	554	! ! Surface pressure gradient
	555	zldg = ( 1._wp - rn_scal_load ) * grav ! local factor
[12377]	556	DO_2D_00_00
	557	zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
	558	zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
	559	END_2D
[11536]	560	IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient
	561	CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters
	562	zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
	563	zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
[4292]	564	ENDIF
	565	!
	566	! Add Coriolis trend:
[6140]	567	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
[4292]	568	! at each time step. We however keep them constant here for optimization.
[11536]	569	! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
	570	CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )
[4292]	571	!
	572	! Add tidal astronomical forcing if defined
[7646]	573	IF ( ln_tide .AND. ln_tide_pot ) THEN
[12377]	574	DO_2D_00_00
	575	zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
	576	zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
	577	END_2D
[4292]	578	ENDIF
	579	!
[9023]	580	! Add bottom stresses:
	581	!jth do implicitly instead
	582	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
[12377]	583	DO_2D_00_00
	584	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
	585	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
	586	END_2D
[11536]	587	ENDIF
[4292]	588	!
	589	! Set next velocities:
[11536]	590	! Compute barotropic speeds at step jit+1 (h : total height of the water colomn)
	591	!-- VECTOR FORM
	592	!-- m+1 m / m+1/2 \ --!
	593	!-- u = u + delta_t' * \ (1-r)g grad_x( ssh') - f * k vect u + frc / --!
	594	!-- --!
	595	!-- FLUX FORM --!
	596	!-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --!
	597	!-- u = m+1 \| h * u + delta_t' * \ h * (1-r)g grad_x( ssh') - h * f * k vect u + h * frc / \| --!
	598	!-- h \ / --!
	599	!------------------------------------------------------------------------------------------------------------------------!
[9023]	600	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
[12377]	601	DO_2D_00_00
	602	ua_e(ji,jj) = ( un_e(ji,jj) &
[12489]	603	& + rDt_e * ( zu_spg(ji,jj) &
[12377]	604	& + zu_trd(ji,jj) &
	605	& + zu_frc(ji,jj) ) &
	606	& ) * ssumask(ji,jj)
[358]	607
[12377]	608	va_e(ji,jj) = ( vn_e(ji,jj) &
[12489]	609	& + rDt_e * ( zv_spg(ji,jj) &
[12377]	610	& + zv_trd(ji,jj) &
	611	& + zv_frc(ji,jj) ) &
	612	& ) * ssvmask(ji,jj)
	613	END_2D
[6140]	614	!
[9023]	615	ELSE !* Flux form
[12377]	616	DO_2D_00_00
	617	! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
	618	! ! backward interpolated depth used in spg terms at jn+1/2
	619	zhu_bck = hu_0(ji,jj) + r1_2r1_e1e2u(ji,jj) ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
	620	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
	621	zhv_bck = hv_0(ji,jj) + r1_2r1_e1e2v(ji,jj) ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
	622	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
	623	! ! inverse depth at jn+1
	624	z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
	625	z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
	626	!
	627	ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) &
[12489]	628	& + rDt_e * ( zhu_bck * zu_spg (ji,jj) & !
[12377]	629	& + zhup2_e(ji,jj) * zu_trd (ji,jj) & !
	630	& + hu(ji,jj,Kmm) * zu_frc (ji,jj) ) ) * z1_hu
	631	!
	632	va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) &
[12489]	633	& + rDt_e * ( zhv_bck * zv_spg (ji,jj) & !
[12377]	634	& + zhvp2_e(ji,jj) * zv_trd (ji,jj) & !
	635	& + hv(ji,jj,Kmm) * zv_frc (ji,jj) ) ) * z1_hv
	636	END_2D
[4292]	637	ENDIF
[10272]	638	!jth implicit bottom friction:
	639	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
[12377]	640	DO_2D_00_00
[12489]	641	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj))
	642	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj))
[12377]	643	END_2D
[10272]	644	ENDIF
[11536]	645
[13229]	646	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
[11536]	647	hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
	648	hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
	649	hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
	650	hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
[13229]	651	ENDIF
	652	! ! open boundaries
	653	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
	654	#if defined key_agrif
	655	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
	656	#endif
	657	!
	658	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
[11536]	659	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp &
[12072]	660	& , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp &
	661	& , hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp )
[11536]	662	ELSE
	663	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
[1438]	664	ENDIF
[4292]	665	!
	666	! !* Swap
	667	! ! ----
[7753]	668	ubb_e (:,:) = ub_e (:,:)
	669	ub_e (:,:) = un_e (:,:)
	670	un_e (:,:) = ua_e (:,:)
	671	!
	672	vbb_e (:,:) = vb_e (:,:)
	673	vb_e (:,:) = vn_e (:,:)
	674	vn_e (:,:) = va_e (:,:)
	675	!
	676	sshbb_e(:,:) = sshb_e(:,:)
	677	sshb_e (:,:) = sshn_e(:,:)
	678	sshn_e (:,:) = ssha_e(:,:)
[4292]	679
	680	! !* Sum over whole bt loop
	681	! ! ----------------------
	682	za1 = wgtbtp1(jn)
[6140]	683	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
[12377]	684	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:)
	685	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:)
[9023]	686	ELSE ! Sum transports
	687	IF ( .NOT.ln_wd_dl ) THEN
[12377]	688	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:)
	689	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:)
[9023]	690	ELSE
[12377]	691	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
	692	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
[9023]	693	END IF
[4292]	694	ENDIF
[9023]	695	! ! Sum sea level
[12377]	696	pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
[9023]	697
[358]	698	! ! ==================== !
	699	END DO ! end loop !
	700	! ! ==================== !
[1438]	701	! -----------------------------------------------------------------------------
[1502]	702	! Phase 3. update the general trend with the barotropic trend
[1438]	703	! -----------------------------------------------------------------------------
[1502]	704	!
[4292]	705	! Set advection velocity correction:
[9023]	706	IF (ln_bt_fw) THEN
[12489]	707	IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
[12377]	708	DO_2D_11_11
	709	zun_save = un_adv(ji,jj)
	710	zvn_save = vn_adv(ji,jj)
	711	! ! apply the previously computed correction
[12489]	712	un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
	713	vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
[12377]	714	! ! Update corrective fluxes for next time step
[12489]	715	un_bf(ji,jj) = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
	716	vn_bf(ji,jj) = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
[12377]	717	! ! Save integrated transport for next computation
	718	ub2_b(ji,jj) = zun_save
	719	vb2_b(ji,jj) = zvn_save
	720	END_2D
[9023]	721	ELSE
[11536]	722	un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero
	723	vn_bf(:,:) = 0._wp
	724	ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation
	725	vb2_b(:,:) = vn_adv(:,:)
	726	END IF
[4292]	727	ENDIF
[9023]	728
	729
[4292]	730	!
	731	! Update barotropic trend:
[6140]	732	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
[4292]	733	DO jk=1,jpkm1
[12489]	734	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b
	735	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b
[4292]	736	END DO
	737	ELSE
[12377]	738	! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
	739	DO_2D_10_10
	740	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
	741	& * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) &
	742	& + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) )
	743	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
	744	& * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) &
	745	& + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) )
	746	END_2D
[10425]	747	CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
[5930]	748	!
[4292]	749	DO jk=1,jpkm1
[12489]	750	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b
	751	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b
[4292]	752	END DO
	753	! Save barotropic velocities not transport:
[12377]	754	puu_b(:,:,Kaa) = puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
	755	pvv_b(:,:,Kaa) = pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
[4292]	756	ENDIF
[9023]	757
	758
	759	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
[4292]	760	DO jk = 1, jpkm1
[12377]	761	puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) umask(:,:,jk)
	762	pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) vmask(:,:,jk)
[358]	763	END DO
[9023]	764
	765	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
	766	DO jk = 1, jpkm1
[12377]	767	puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
	768	& + zuwdav2(:,:)(puu(:,:,jk,Kmm) - un_adv(:,:)r1_hu(:,:,Kmm)) ) * umask(:,:,jk)
	769	pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) &
	770	& + zvwdav2(:,:)(pvv(:,:,jk,Kmm) - vn_adv(:,:)r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)
[9023]	771	END DO
	772	END IF
	773
	774
[12377]	775	CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current
	776	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current
[1502]	777	!
[4486]	778	#if defined key_agrif
	779	! Save time integrated fluxes during child grid integration
[5656]	780	! (used to update coarse grid transports at next time step)
[4486]	781	!
[6140]	782	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
	783	IF( Agrif_NbStepint() == 0 ) THEN
[7753]	784	ub2_i_b(:,:) = 0._wp
	785	vb2_i_b(:,:) = 0._wp
[4486]	786	END IF
	787	!
	788	za1 = 1._wp / REAL(Agrif_rhot(), wp)
[7753]	789	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
	790	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
[4486]	791	ENDIF
	792	#endif
[1502]	793	! !* write time-spliting arrays in the restart
[6140]	794	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
[508]	795	!
[9023]	796	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
	797	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
[1662]	798	!
[12377]	799	CALL iom_put( "baro_u" , puu_b(:,:,Kmm) ) ! Barotropic U Velocity
	800	CALL iom_put( "baro_v" , pvv_b(:,:,Kmm) ) ! Barotropic V Velocity
[2715]	801	!
[508]	802	END SUBROUTINE dyn_spg_ts
	803
[6140]	804
[4292]	805	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
	806	!!---------------------------------------------------------------------
	807	!! * ROUTINE ts_wgt *
	808	!!
	809	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
	810	!!----------------------------------------------------------------------
	811	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
	812	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
	813	INTEGER, INTENT(inout) :: jpit ! cycle length
[12489]	814	REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, & ! Primary weights
[4292]	815	zwgt2 ! Secondary weights
	816
	817	INTEGER :: jic, jn, ji ! temporary integers
	818	REAL(wp) :: za1, za2
	819	!!----------------------------------------------------------------------
[508]	820
[4292]	821	zwgt1(:) = 0._wp
	822	zwgt2(:) = 0._wp
	823
	824	! Set time index when averaged value is requested
	825	IF (ll_fw) THEN
[12489]	826	jic = nn_e
[4292]	827	ELSE
[12489]	828	jic = 2 * nn_e
[4292]	829	ENDIF
	830
	831	! Set primary weights:
	832	IF (ll_av) THEN
	833	! Define simple boxcar window for primary weights
[12489]	834	! (width = nn_e, centered around jic)
[4292]	835	SELECT CASE ( nn_bt_flt )
	836	CASE( 0 ) ! No averaging
	837	zwgt1(jic) = 1._wp
	838	jpit = jic
	839
[12489]	840	CASE( 1 ) ! Boxcar, width = nn_e
	841	DO jn = 1, 3*nn_e
	842	za1 = ABS(float(jn-jic))/float(nn_e)
[4292]	843	IF (za1 < 0.5_wp) THEN
	844	zwgt1(jn) = 1._wp
	845	jpit = jn
	846	ENDIF
	847	ENDDO
	848
[12489]	849	CASE( 2 ) ! Boxcar, width = 2 * nn_e
	850	DO jn = 1, 3*nn_e
	851	za1 = ABS(float(jn-jic))/float(nn_e)
[4292]	852	IF (za1 < 1._wp) THEN
	853	zwgt1(jn) = 1._wp
	854	jpit = jn
	855	ENDIF
	856	ENDDO
	857	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
	858	END SELECT
	859
	860	ELSE ! No time averaging
	861	zwgt1(jic) = 1._wp
	862	jpit = jic
	863	ENDIF
	864
	865	! Set secondary weights
	866	DO jn = 1, jpit
	867	DO ji = jn, jpit
	868	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
	869	END DO
	870	END DO
	871
	872	! Normalize weigths:
	873	za1 = 1._wp / SUM(zwgt1(1:jpit))
	874	za2 = 1._wp / SUM(zwgt2(1:jpit))
	875	DO jn = 1, jpit
	876	zwgt1(jn) = zwgt1(jn) * za1
	877	zwgt2(jn) = zwgt2(jn) * za2
	878	END DO
	879	!
	880	END SUBROUTINE ts_wgt
	881
[6140]	882
[508]	883	SUBROUTINE ts_rst( kt, cdrw )
	884	!!---------------------------------------------------------------------
	885	!! * ROUTINE ts_rst *
	886	!!
	887	!! ** Purpose : Read or write time-splitting arrays in restart file
	888	!!----------------------------------------------------------------------
[9528]	889	INTEGER , INTENT(in) :: kt ! ocean time-step
	890	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
[508]	891	!!----------------------------------------------------------------------
	892	!
[9506]	893	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
	894	! ! ---------------
[12489]	895	IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN !* Read the restart file
[12738]	896	CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	897	CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
	898	CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	899	CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
[9506]	900	IF( .NOT.ln_bt_av ) THEN
[12738]	901	CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp, ldxios = lrxios )
	902	CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	903	CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
	904	CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp, ldxios = lrxios )
	905	CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	906	CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
[9506]	907	ENDIF
[4486]	908	#if defined key_agrif
[9506]	909	! Read time integrated fluxes
	910	IF ( .NOT.Agrif_Root() ) THEN
[12738]	911	CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	912	CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
[9506]	913	ENDIF
	914	#endif
	915	ELSE !* Start from rest
	916	IF(lwp) WRITE(numout,*)
	917	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
	918	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
	919	un_adv(:,:) = 0._wp ; vn_adv(:,:) = 0._wp ! used in the 1st interpol of agrif
	920	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
	921	#if defined key_agrif
	922	IF ( .NOT.Agrif_Root() ) THEN
	923	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
	924	ENDIF
	925	#endif
[4486]	926	ENDIF
[9506]	927	!
	928	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
	929	! ! -------------------
	930	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
[9367]	931	IF( lwxios ) CALL iom_swap( cwxios_context )
	932	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios )
	933	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios )
	934	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios )
	935	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios )
[4292]	936	!
	937	IF (.NOT.ln_bt_av) THEN
[9367]	938	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:), ldxios = lwxios )
	939	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:), ldxios = lwxios )
	940	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:), ldxios = lwxios )
	941	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:), ldxios = lwxios )
	942	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:), ldxios = lwxios )
	943	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:), ldxios = lwxios )
[4292]	944	ENDIF
[4486]	945	#if defined key_agrif
	946	! Save time integrated fluxes
	947	IF ( .NOT.Agrif_Root() ) THEN
[9367]	948	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lwxios )
	949	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lwxios )
[4486]	950	ENDIF
	951	#endif
[9367]	952	IF( lwxios ) CALL iom_swap( cxios_context )
[4292]	953	ENDIF
	954	!
	955	END SUBROUTINE ts_rst
[2528]	956
[6140]	957
	958	SUBROUTINE dyn_spg_ts_init
[4292]	959	!!---------------------------------------------------------------------
	960	!! * ROUTINE dyn_spg_ts_init *
	961	!!
	962	!! ** Purpose : Set time splitting options
	963	!!----------------------------------------------------------------------
[6140]	964	INTEGER :: ji ,jj ! dummy loop indices
	965	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
[9019]	966	REAL(wp), DIMENSION(jpi,jpj) :: zcu
[4292]	967	!!----------------------------------------------------------------------
[4370]	968	!
[5930]	969	! Max courant number for ext. grav. waves
[4370]	970	!
[13065]	971	DO_2D_00_00
[12377]	972	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
	973	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
	974	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
	975	END_2D
[5930]	976	!
[13065]	977	zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) )
[10425]	978	CALL mpp_max( 'dynspg_ts', zcmax )
[2528]	979
[4370]	980	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
[12489]	981	IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
[4292]	982
[12489]	983	rDt_e = rn_Dt / REAL( nn_e , wp )
	984	zcmax = zcmax * rDt_e
[9023]	985	! Print results
[4292]	986	IF(lwp) WRITE(numout,*)
[9169]	987	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
	988	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
[5930]	989	IF( ln_bt_auto ) THEN
[12489]	990	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e '
[4370]	991	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
[4292]	992	ELSE
[12489]	993	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e
[358]	994	ENDIF
[4292]	995
	996	IF(ln_bt_av) THEN
[12489]	997	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on '
[4292]	998	ELSE
[9169]	999	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
[4292]	1000	ENDIF
[508]	1001	!
[4292]	1002	!
	1003	IF(ln_bt_fw) THEN
[4370]	1004	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
[4292]	1005	ELSE
[4370]	1006	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
[4292]	1007	ENDIF
	1008	!
[4486]	1009	#if defined key_agrif
	1010	! Restrict the use of Agrif to the forward case only
[9023]	1011	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
[4486]	1012	#endif
	1013	!
[4370]	1014	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
[4292]	1015	SELECT CASE ( nn_bt_flt )
[6140]	1016	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
[12489]	1017	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e'
	1018	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_e'
[9169]	1019	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
[4292]	1020	END SELECT
	1021	!
[4370]	1022	IF(lwp) WRITE(numout,*) ' '
[12489]	1023	IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e
	1024	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rDt_e
[4370]	1025	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
	1026	!
[9023]	1027	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
	1028	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
	1029	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
	1030	ENDIF
	1031	!
[6140]	1032	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
[4292]	1033	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
	1034	ENDIF
[6140]	1035	IF( zcmax>0.9_wp ) THEN
[12489]	1036	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )
[4292]	1037	ENDIF
	1038	!
[9124]	1039	! ! Allocate time-splitting arrays
	1040	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
	1041	!
	1042	! ! read restart when needed
[9506]	1043	CALL ts_rst( nit000, 'READ' )
[9124]	1044	!
[9367]	1045	IF( lwxios ) THEN
	1046	! define variables in restart file when writing with XIOS
	1047	CALL iom_set_rstw_var_active('ub2_b')
	1048	CALL iom_set_rstw_var_active('vb2_b')
	1049	CALL iom_set_rstw_var_active('un_bf')
	1050	CALL iom_set_rstw_var_active('vn_bf')
	1051	!
	1052	IF (.NOT.ln_bt_av) THEN
	1053	CALL iom_set_rstw_var_active('sshbb_e')
	1054	CALL iom_set_rstw_var_active('ubb_e')
	1055	CALL iom_set_rstw_var_active('vbb_e')
	1056	CALL iom_set_rstw_var_active('sshb_e')
	1057	CALL iom_set_rstw_var_active('ub_e')
	1058	CALL iom_set_rstw_var_active('vb_e')
	1059	ENDIF
	1060	#if defined key_agrif
	1061	! Save time integrated fluxes
	1062	IF ( .NOT.Agrif_Root() ) THEN
	1063	CALL iom_set_rstw_var_active('ub2_i_b')
	1064	CALL iom_set_rstw_var_active('vb2_i_b')
	1065	ENDIF
	1066	#endif
	1067	ENDIF
	1068	!
[4292]	1069	END SUBROUTINE dyn_spg_ts_init
[508]	1070
[11536]	1071
[12377]	1072	SUBROUTINE dyn_cor_2D_init( Kmm )
[11536]	1073	!!---------------------------------------------------------------------
[12377]	1074	!! * ROUTINE dyn_cor_2D_init *
[11536]	1075	!!
	1076	!! ** Purpose : Set time splitting options
	1077	!! Set arrays to remove/compute coriolis trend.
	1078	!! Do it once during initialization if volume is fixed, else at each long time step.
	1079	!! Note that these arrays are also used during barotropic loop. These are however frozen
	1080	!! although they should be updated in the variable volume case. Not a big approximation.
	1081	!! To remove this approximation, copy lines below inside barotropic loop
	1082	!! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
	1083	!!
	1084	!! Compute zwz = f / ( height of the water colomn )
	1085	!!----------------------------------------------------------------------
[12377]	1086	INTEGER, INTENT(in) :: Kmm ! Time index
[11536]	1087	INTEGER :: ji ,jj, jk ! dummy loop indices
	1088	REAL(wp) :: z1_ht
	1089	REAL(wp), DIMENSION(jpi,jpj) :: zhf
	1090	!!----------------------------------------------------------------------
	1091	!
	1092	SELECT CASE( nvor_scheme )
	1093	CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)
	1094	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
	1095	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
[12377]	1096	DO_2D_10_10
	1097	zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + &
	1098	& ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp
	1099	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
	1100	END_2D
[11536]	1101	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
[12377]	1102	DO_2D_10_10
	1103	zwz(ji,jj) = ( ht (ji ,jj+1) + ht (ji+1,jj+1) &
	1104	& + ht (ji ,jj ) + ht (ji+1,jj ) ) &
	1105	& / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) &
	1106	& + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) )
	1107	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
	1108	END_2D
[11536]	1109	END SELECT
	1110	CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
	1111	!
	1112	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[12377]	1113	DO_2D_01_01
	1114	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	1115	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	1116	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	1117	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	1118	END_2D
[11536]	1119	!
	1120	CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme)
	1121	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[12377]	1122	DO_2D_01_01
	1123	z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
	1124	ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht
	1125	ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht
	1126	ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
	1127	ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht
	1128	END_2D
[11536]	1129	!
	1130	CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !
	1131	!
	1132	zwz(:,:) = 0._wp
	1133	zhf(:,:) = 0._wp
	1134
	1135	!!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
	1136	!!gm A priori a better value should be something like :
	1137	!!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
	1138	!!gm divided by the sum of the corresponding mask
	1139	!!gm
	1140	!!
	1141	IF( .NOT.ln_sco ) THEN
	1142
	1143	!!gm agree the JC comment : this should be done in a much clear way
	1144
	1145	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
	1146	! Set it to zero for the time being
	1147	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
	1148	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
	1149	! ENDIF
	1150	! zhf(:,:) = gdepw_0(:,:,jk+1)
	1151	!
	1152	ELSE
	1153	!
	1154	!zhf(:,:) = hbatf(:,:)
[12377]	1155	DO_2D_10_10
	1156	zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &
	1157	& + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &
	1158	& / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &
	1159	& + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )
	1160	END_2D
[11536]	1161	ENDIF
	1162	!
	1163	DO jj = 1, jpjm1
	1164	zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
	1165	END DO
	1166	!
	1167	DO jk = 1, jpkm1
	1168	DO jj = 1, jpjm1
[12377]	1169	zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
[11536]	1170	END DO
	1171	END DO
	1172	CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
	1173	! JC: TBC. hf should be greater than 0
[12377]	1174	DO_2D_11_11
	1175	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj)
	1176	END_2D
[11536]	1177	zwz(:,:) = ff_f(:,:) * zwz(:,:)
	1178	END SELECT
	1179
	1180	END SUBROUTINE dyn_cor_2d_init
	1181
	1182
	1183
[12377]	1184	SUBROUTINE dyn_cor_2d( phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )
[11536]	1185	!!---------------------------------------------------------------------
	1186	!! * ROUTINE dyn_cor_2d *
	1187	!!
	1188	!! ** Purpose : Compute u and v coriolis trends
	1189	!!----------------------------------------------------------------------
	1190	INTEGER :: ji ,jj ! dummy loop indices
	1191	REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - -
[12377]	1192	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: phu, phv, punb, pvnb, zhU, zhV
[11536]	1193	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd
	1194	!!----------------------------------------------------------------------
	1195	SELECT CASE( nvor_scheme )
	1196	CASE( np_ENT ) ! enstrophy conserving scheme (f-point)
[12377]	1197	DO_2D_00_00
	1198	z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
	1199	z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
	1200	zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu &
	1201	& * ( e1e2t(ji+1,jj)ht(ji+1,jj)ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) &
	1202	& + e1e2t(ji ,jj)ht(ji ,jj)ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) )
	1203	!
	1204	zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv &
	1205	& * ( e1e2t(ji,jj+1)ht(ji,jj+1)ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) &
	1206	& + e1e2t(ji,jj )ht(ji,jj )ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) )
	1207	END_2D
[11536]	1208	!
	1209	CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX
[12377]	1210	DO_2D_00_00
	1211	zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
	1212	zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
	1213	zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
	1214	zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
	1215	! energy conserving formulation for planetary vorticity term
	1216	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	1217	zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
	1218	END_2D
[11536]	1219	!
	1220	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
[12377]	1221	DO_2D_00_00
	1222	zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) &
	1223	& + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
	1224	zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) &
	1225	& + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
	1226	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	1227	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
	1228	END_2D
[11536]	1229	!
	1230	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
[12377]	1231	DO_2D_00_00
	1232	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) &
	1233	& + ftnw(ji+1,jj) * zhV(ji+1,jj ) &
	1234	& + ftse(ji,jj ) * zhV(ji ,jj-1) &
	1235	& + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
	1236	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
	1237	& + ftse(ji,jj+1) * zhU(ji ,jj+1) &
	1238	& + ftnw(ji,jj ) * zhU(ji-1,jj ) &
	1239	& + ftne(ji,jj ) * zhU(ji ,jj ) )
	1240	END_2D
[11536]	1241	!
	1242	END SELECT
	1243	!
	1244	END SUBROUTINE dyn_cor_2D
	1245
	1246
	1247	SUBROUTINE wad_tmsk( pssh, ptmsk )
	1248	!!----------------------------------------------------------------------
	1249	!! * ROUTINE wad_lmt *
	1250	!!
	1251	!! ** Purpose : set wetting & drying mask at tracer points
	1252	!! for the current barotropic sub-step
	1253	!!
	1254	!! ** Method : ???
	1255	!!
	1256	!! ** Action : ptmsk : wetting & drying t-mask
	1257	!!----------------------------------------------------------------------
	1258	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh !
	1259	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk !
	1260	!
	1261	INTEGER :: ji, jj ! dummy loop indices
	1262	!!----------------------------------------------------------------------
	1263	!
	1264	IF( ln_wd_dl_rmp ) THEN
[12377]	1265	DO_2D_11_11
	1266	IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
	1267	! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
	1268	ptmsk(ji,jj) = 1._wp
	1269	ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
	1270	ptmsk(ji,jj) = TANH( 50._wp( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1) )
	1271	ELSE
	1272	ptmsk(ji,jj) = 0._wp
	1273	ENDIF
	1274	END_2D
[11536]	1275	ELSE
[12377]	1276	DO_2D_11_11
	1277	IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp
	1278	ELSE ; ptmsk(ji,jj) = 0._wp
	1279	ENDIF
	1280	END_2D
[11536]	1281	ENDIF
	1282	!
	1283	END SUBROUTINE wad_tmsk
	1284
	1285
	1286	SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
	1287	!!----------------------------------------------------------------------
	1288	!! * ROUTINE wad_lmt *
	1289	!!
	1290	!! ** Purpose : set wetting & drying mask at tracer points
	1291	!! for the current barotropic sub-step
	1292	!!
	1293	!! ** Method : ???
	1294	!!
	1295	!! ** Action : ptmsk : wetting & drying t-mask
	1296	!!----------------------------------------------------------------------
	1297	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
	1298	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
	1299	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
	1300	!
	1301	INTEGER :: ji, jj ! dummy loop indices
	1302	!!----------------------------------------------------------------------
	1303	!
[12377]	1304	DO_2D_11_10
	1305	IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj)
	1306	ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj)
	1307	ENDIF
	1308	phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
	1309	pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
	1310	END_2D
[11536]	1311	!
[12377]	1312	DO_2D_10_11
	1313	IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj )
	1314	ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1)
	1315	ENDIF
	1316	phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
	1317	pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
	1318	END_2D
[11536]	1319	!
	1320	END SUBROUTINE wad_Umsk
	1321
	1322
[12377]	1323	SUBROUTINE wad_spg( pshn, zcpx, zcpy )
[11536]	1324	!!---------------------------------------------------------------------
	1325	!! * ROUTINE wad_sp *
	1326	!!
	1327	!! ** Purpose :
	1328	!!----------------------------------------------------------------------
	1329	INTEGER :: ji ,jj ! dummy loop indices
	1330	LOGICAL :: ll_tmp1, ll_tmp2
[12377]	1331	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn
[11536]	1332	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
	1333	!!----------------------------------------------------------------------
[12377]	1334	DO_2D_00_00
	1335	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > &
	1336	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
	1337	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji+1,jj) + ht_0(ji+1,jj) ) &
	1338	& > rn_wdmin1 + rn_wdmin2
	1339	ll_tmp2 = ( ABS( pshn(ji+1,jj) - pshn(ji ,jj)) > 1.E-12 ).AND.( &
	1340	& MAX( pshn(ji,jj) , pshn(ji+1,jj) ) > &
	1341	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
	1342	IF(ll_tmp1) THEN
	1343	zcpx(ji,jj) = 1.0_wp
	1344	ELSEIF(ll_tmp2) THEN
	1345	! no worries about pshn(ji+1,jj) - pshn(ji ,jj) = 0, it won't happen ! here
	1346	zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
	1347	& / (pshn(ji+1,jj) - pshn(ji ,jj)) )
	1348	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
	1349	ELSE
	1350	zcpx(ji,jj) = 0._wp
	1351	ENDIF
	1352	!
	1353	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji,jj+1) ) > &
	1354	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
	1355	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji,jj+1) + ht_0(ji,jj+1) ) &
	1356	& > rn_wdmin1 + rn_wdmin2
	1357	ll_tmp2 = ( ABS( pshn(ji,jj) - pshn(ji,jj+1)) > 1.E-12 ).AND.( &
	1358	& MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > &
	1359	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
	1360
	1361	IF(ll_tmp1) THEN
	1362	zcpy(ji,jj) = 1.0_wp
	1363	ELSE IF(ll_tmp2) THEN
	1364	! no worries about pshn(ji,jj+1) - pshn(ji,jj ) = 0, it won't happen ! here
	1365	zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
	1366	& / (pshn(ji,jj+1) - pshn(ji,jj )) )
	1367	zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) )
	1368	ELSE
	1369	zcpy(ji,jj) = 0._wp
	1370	ENDIF
	1371	END_2D
[11536]	1372
	1373	END SUBROUTINE wad_spg
	1374
	1375
	1376
[12377]	1377	SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
[11536]	1378	!!----------------------------------------------------------------------
	1379	!! * ROUTINE dyn_drg_init *
	1380	!!
	1381	!! ** Purpose : - add the baroclinic top/bottom drag contribution to
	1382	!! the baroclinic part of the barotropic RHS
	1383	!! - compute the barotropic drag coefficients
	1384	!!
	1385	!! ** Method : computation done over the INNER domain only
	1386	!!----------------------------------------------------------------------
[12377]	1387	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
	1388	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation
	1389	REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels
	1390	REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS
	1391	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients
[11536]	1392	!
	1393	INTEGER :: ji, jj ! dummy loop indices
	1394	INTEGER :: ikbu, ikbv, iktu, iktv
	1395	REAL(wp) :: zztmp
	1396	REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i
	1397	!!----------------------------------------------------------------------
	1398	!
	1399	! !== Set the barotropic drag coef. ==!
	1400	!
	1401	IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)
	1402
[12377]	1403	DO_2D_00_00
	1404	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
	1405	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
	1406	END_2D
[11536]	1407	ELSE ! bottom friction only
[12377]	1408	DO_2D_00_00
	1409	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
	1410	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
	1411	END_2D
[11536]	1412	ENDIF
	1413	!
	1414	! !== BOTTOM stress contribution from baroclinic velocities ==!
	1415	!
	1416	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
	1417
[12377]	1418	DO_2D_00_00
	1419	ikbu = mbku(ji,jj)
	1420	ikbv = mbkv(ji,jj)
	1421	zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
	1422	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
	1423	END_2D
[11536]	1424	ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
	1425
[12377]	1426	DO_2D_00_00
	1427	ikbu = mbku(ji,jj)
	1428	ikbv = mbkv(ji,jj)
	1429	zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
	1430	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
	1431	END_2D
[11536]	1432	ENDIF
	1433	!
	1434	IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please !
[12489]	1435	zztmp = -1._wp / rDt_e
[12377]	1436	DO_2D_00_00
	1437	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
	1438	& r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp )
	1439	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &
	1440	& r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
	1441	END_2D
[11536]	1442	ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
	1443
[12377]	1444	DO_2D_00_00
	1445	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) zu_i(ji,jj)
	1446	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) zv_i(ji,jj)
	1447	END_2D
[11536]	1448	END IF
	1449	!
	1450	! !== TOP stress contribution from baroclinic velocities ==! (no W/D case)
	1451	!
	1452	IF( ln_isfcav ) THEN
	1453	!
	1454	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
	1455
[12377]	1456	DO_2D_00_00
	1457	iktu = miku(ji,jj)
	1458	iktv = mikv(ji,jj)
	1459	zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
	1460	zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
	1461	END_2D
[11536]	1462	ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
	1463
[12377]	1464	DO_2D_00_00
	1465	iktu = miku(ji,jj)
	1466	iktv = mikv(ji,jj)
	1467	zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
	1468	zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
	1469	END_2D
[11536]	1470	ENDIF
	1471	!
	1472	! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
	1473
[12377]	1474	DO_2D_00_00
	1475	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) zu_i(ji,jj)
	1476	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) zv_i(ji,jj)
	1477	END_2D
[11536]	1478	!
	1479	ENDIF
	1480	!
	1481	END SUBROUTINE dyn_drg_init
	1482
	1483	SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in
	1484	& za0, za1, za2, za3 ) ! ==> out
	1485	!!----------------------------------------------------------------------
	1486	INTEGER ,INTENT(in ) :: jn ! index of sub time step
	1487	LOGICAL ,INTENT(in ) :: ll_init !
	1488	REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient
	1489	!
	1490	REAL(wp) :: zepsilon, zgamma ! - -
	1491	!!----------------------------------------------------------------------
	1492	! ! set Half-step back interpolation coefficient
	1493	IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward
	1494	za0 = 1._wp
	1495	za1 = 0._wp
	1496	za2 = 0._wp
	1497	za3 = 0._wp
	1498	ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
	1499	za0 = 1.0833333333333_wp ! za0 = 1-gam-eps
	1500	za1 =-0.1666666666666_wp ! za1 = gam
	1501	za2 = 0.0833333333333_wp ! za2 = eps
	1502	za3 = 0._wp
	1503	ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
	1504	IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion
	1505	za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps
	1506	za1 = 0.285_wp ! za1 = 1/2 - 2gam - 3eps
	1507	za2 = 0.088_wp ! za2 = gam
	1508	za3 = 0.013_wp ! za3 = eps
	1509	ELSE ! no time diffusion
	1510	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
	1511	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
	1512	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
	1513	za1 = 1._wp - za0 - zgamma - zepsilon
	1514	za2 = zgamma
	1515	za3 = zepsilon
	1516	ENDIF
	1517	ENDIF
	1518	END SUBROUTINE ts_bck_interp
	1519
	1520
[358]	1521	!!======================================================================
	1522	END MODULE dynspg_ts

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

Download in other formats: