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_r12377_KERNEL-06_techene_e3/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynspg_ts.F90 @ 13193

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

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

Download in other formats: