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_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/DYN/dynspg_ts.F90 @ 13886

Last change on this file since 13886 was 13886, checked in by rlod, 3 years ago
phasing with trunk at revision r13787
Property svn:keywords set to `Id`
File size: 79.6 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"
[13237]	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	!
[12489]	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
[12489]	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
[13237]	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
[12489]	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
[12489]	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	!
[12489]	200	IF( l_1st_euler ) ll_init=.TRUE.
[1502]	201	!
[12489]	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	!
[12489]	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	! ! --------------------------- !
[13237]	231	DO jk = 1 , jpk
	232	ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
	233	ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
	234	END DO
[1502]	235	!
[13237]	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	!
[13237]	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	!
[13237]	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
[13472]	266	DO_2D( 0, 0, 0, 0 ) ! SPG with the application of W/D gravity filters
[12377]	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
[13295]	273	DO_2D( 0, 0, 0, 0 )
[12377]	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	!
[13497]	281	DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend
[12377]	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)
[13295]	293	DO_2D( 0, 0, 0, 0 )
[12377]	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
[13295]	299	DO_2D( 0, 0, 0, 0 )
[12377]	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
[13295]	311	DO_2D( 0, 0, 0, 0 )
[12489]	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
[12489]	316	zztmp = r1_rho0 * r1_2
[13295]	317	DO_2D( 0, 0, 0, 0 )
[12377]	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)
[12489]	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)
[12489]	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)
	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
[12489]	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)
[13497]	477	DO_2D( 1, 1, 1, 0 ) ! not jpi-column
[12377]	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
[13497]	482	DO_2D( 1, 0, 1, 1 ) ! not jpj-row
[12377]	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
[12489]	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	!-------------------------------------------------------------------------!
[13295]	517	DO_2D( 0, 0, 0, 0 )
[12377]	518	zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj)
[12489]	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	!
[13289]	522	CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )
	523	!
[13216]	524	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
	525	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
	526	#if defined key_agrif
	527	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
	528	#endif
[11536]	529	!
	530	! ! Sum over sub-time-steps to compute advective velocities
	531	za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5
	532	un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
	533	vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
	534	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
	535	IF ( ln_wd_dl_bc ) THEN
	536	zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column
	537	zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row
	538	END IF
	539	!
[4292]	540	!
	541	! Sea Surface Height at u-,v-points (vvl case only)
[6140]	542	IF( .NOT.ln_linssh ) THEN
[13295]	543	DO_2D( 0, 0, 0, 0 )
[12377]	544	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
	545	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
	546	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
	547	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
	548	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
	549	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
	550	END_2D
[4292]	551	ENDIF
[11536]	552	!
	553	! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
	554	!-- m+1/2 m+1 m m-1 m-2 --!
	555	!-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --!
	556	!------------------------------------------------------------------------------------------!
	557	CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation
[9528]	558	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
	559	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
[1502]	560	!
[11536]	561	! ! Surface pressure gradient
	562	zldg = ( 1._wp - rn_scal_load ) * grav ! local factor
[13295]	563	DO_2D( 0, 0, 0, 0 )
[12377]	564	zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
	565	zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
	566	END_2D
[11536]	567	IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient
	568	CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters
	569	zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
	570	zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
[4292]	571	ENDIF
	572	!
	573	! Add Coriolis trend:
[6140]	574	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
[4292]	575	! at each time step. We however keep them constant here for optimization.
[11536]	576	! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
[13237]	577	CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )
[4292]	578	!
	579	! Add tidal astronomical forcing if defined
[7646]	580	IF ( ln_tide .AND. ln_tide_pot ) THEN
[13295]	581	DO_2D( 0, 0, 0, 0 )
[12377]	582	zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
	583	zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
	584	END_2D
[4292]	585	ENDIF
	586	!
[9023]	587	! Add bottom stresses:
	588	!jth do implicitly instead
	589	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
[13295]	590	DO_2D( 0, 0, 0, 0 )
[12377]	591	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
	592	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
	593	END_2D
[11536]	594	ENDIF
[4292]	595	!
	596	! Set next velocities:
[11536]	597	! Compute barotropic speeds at step jit+1 (h : total height of the water colomn)
	598	!-- VECTOR FORM
	599	!-- m+1 m / m+1/2 \ --!
	600	!-- u = u + delta_t' * \ (1-r)g grad_x( ssh') - f * k vect u + frc / --!
	601	!-- --!
	602	!-- FLUX FORM --!
	603	!-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --!
	604	!-- u = m+1 \| h * u + delta_t' * \ h * (1-r)g grad_x( ssh') - h * f * k vect u + h * frc / \| --!
	605	!-- h \ / --!
	606	!------------------------------------------------------------------------------------------------------------------------!
[9023]	607	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
[13295]	608	DO_2D( 0, 0, 0, 0 )
[12377]	609	ua_e(ji,jj) = ( un_e(ji,jj) &
[12489]	610	& + rDt_e * ( zu_spg(ji,jj) &
[12377]	611	& + zu_trd(ji,jj) &
	612	& + zu_frc(ji,jj) ) &
	613	& ) * ssumask(ji,jj)
[358]	614
[12377]	615	va_e(ji,jj) = ( vn_e(ji,jj) &
[12489]	616	& + rDt_e * ( zv_spg(ji,jj) &
[12377]	617	& + zv_trd(ji,jj) &
	618	& + zv_frc(ji,jj) ) &
	619	& ) * ssvmask(ji,jj)
	620	END_2D
[6140]	621	!
[9023]	622	ELSE !* Flux form
[13295]	623	DO_2D( 0, 0, 0, 0 )
[12377]	624	! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
	625	! ! backward interpolated depth used in spg terms at jn+1/2
	626	zhu_bck = hu_0(ji,jj) + r1_2r1_e1e2u(ji,jj) ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
	627	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
	628	zhv_bck = hv_0(ji,jj) + r1_2r1_e1e2v(ji,jj) ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
	629	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
	630	! ! inverse depth at jn+1
	631	z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
	632	z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
	633	!
	634	ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) &
[12489]	635	& + rDt_e * ( zhu_bck * zu_spg (ji,jj) & !
[12377]	636	& + zhup2_e(ji,jj) * zu_trd (ji,jj) & !
	637	& + hu(ji,jj,Kmm) * zu_frc (ji,jj) ) ) * z1_hu
	638	!
	639	va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) &
[12489]	640	& + rDt_e * ( zhv_bck * zv_spg (ji,jj) & !
[12377]	641	& + zhvp2_e(ji,jj) * zv_trd (ji,jj) & !
	642	& + hv(ji,jj,Kmm) * zv_frc (ji,jj) ) ) * z1_hv
	643	END_2D
[4292]	644	ENDIF
[10272]	645	!jth implicit bottom friction:
	646	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
[13295]	647	DO_2D( 0, 0, 0, 0 )
[12489]	648	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj))
	649	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj))
[12377]	650	END_2D
[10272]	651	ENDIF
[11536]	652
[13216]	653	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
[11536]	654	hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
	655	hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
	656	hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
	657	hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) )
[13216]	658	ENDIF
	659	!
	660	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
[11536]	661	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp &
[12072]	662	& , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp &
	663	& , hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp )
[11536]	664	ELSE
	665	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
[1438]	666	ENDIF
[13289]	667	! ! open boundaries
	668	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
	669	#if defined key_agrif
	670	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
	671	#endif
[4292]	672	! !* Swap
	673	! ! ----
[7753]	674	ubb_e (:,:) = ub_e (:,:)
	675	ub_e (:,:) = un_e (:,:)
	676	un_e (:,:) = ua_e (:,:)
	677	!
	678	vbb_e (:,:) = vb_e (:,:)
	679	vb_e (:,:) = vn_e (:,:)
	680	vn_e (:,:) = va_e (:,:)
	681	!
	682	sshbb_e(:,:) = sshb_e(:,:)
	683	sshb_e (:,:) = sshn_e(:,:)
	684	sshn_e (:,:) = ssha_e(:,:)
[4292]	685
	686	! !* Sum over whole bt loop
	687	! ! ----------------------
	688	za1 = wgtbtp1(jn)
[6140]	689	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
[12377]	690	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:)
	691	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:)
[9023]	692	ELSE ! Sum transports
	693	IF ( .NOT.ln_wd_dl ) THEN
[12377]	694	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:)
	695	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:)
[9023]	696	ELSE
[12377]	697	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
	698	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
[9023]	699	END IF
[4292]	700	ENDIF
[9023]	701	! ! Sum sea level
[12377]	702	pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
[9023]	703
[358]	704	! ! ==================== !
	705	END DO ! end loop !
	706	! ! ==================== !
[1438]	707	! -----------------------------------------------------------------------------
[1502]	708	! Phase 3. update the general trend with the barotropic trend
[1438]	709	! -----------------------------------------------------------------------------
[1502]	710	!
[4292]	711	! Set advection velocity correction:
[9023]	712	IF (ln_bt_fw) THEN
[12489]	713	IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
[13295]	714	DO_2D( 1, 1, 1, 1 )
[12377]	715	zun_save = un_adv(ji,jj)
	716	zvn_save = vn_adv(ji,jj)
	717	! ! apply the previously computed correction
[12489]	718	un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
	719	vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
[12377]	720	! ! Update corrective fluxes for next time step
[12489]	721	un_bf(ji,jj) = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
	722	vn_bf(ji,jj) = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
[12377]	723	! ! Save integrated transport for next computation
	724	ub2_b(ji,jj) = zun_save
	725	vb2_b(ji,jj) = zvn_save
	726	END_2D
[9023]	727	ELSE
[11536]	728	un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero
	729	vn_bf(:,:) = 0._wp
	730	ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation
	731	vb2_b(:,:) = vn_adv(:,:)
	732	END IF
[4292]	733	ENDIF
[9023]	734
	735
[4292]	736	!
	737	! Update barotropic trend:
[6140]	738	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
[4292]	739	DO jk=1,jpkm1
[12489]	740	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b
	741	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b
[4292]	742	END DO
	743	ELSE
[12377]	744	! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
[13295]	745	DO_2D( 1, 0, 1, 0 )
[12377]	746	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
	747	& * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) &
	748	& + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) )
	749	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
	750	& * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) &
	751	& + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) )
	752	END_2D
[10425]	753	CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
[5930]	754	!
[4292]	755	DO jk=1,jpkm1
[12489]	756	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b
	757	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b
[4292]	758	END DO
	759	! Save barotropic velocities not transport:
[12377]	760	puu_b(:,:,Kaa) = puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
	761	pvv_b(:,:,Kaa) = pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
[4292]	762	ENDIF
[9023]	763
	764
	765	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
[4292]	766	DO jk = 1, jpkm1
[12377]	767	puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) umask(:,:,jk)
	768	pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) vmask(:,:,jk)
[358]	769	END DO
[9023]	770
	771	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
	772	DO jk = 1, jpkm1
[12377]	773	puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
	774	& + zuwdav2(:,:)(puu(:,:,jk,Kmm) - un_adv(:,:)r1_hu(:,:,Kmm)) ) * umask(:,:,jk)
	775	pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) &
	776	& + zvwdav2(:,:)(pvv(:,:,jk,Kmm) - vn_adv(:,:)r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)
[9023]	777	END DO
	778	END IF
	779
	780
[12377]	781	CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current
	782	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current
[1502]	783	!
[4486]	784	#if defined key_agrif
	785	! Save time integrated fluxes during child grid integration
[5656]	786	! (used to update coarse grid transports at next time step)
[4486]	787	!
[6140]	788	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
	789	IF( Agrif_NbStepint() == 0 ) THEN
[7753]	790	ub2_i_b(:,:) = 0._wp
	791	vb2_i_b(:,:) = 0._wp
[4486]	792	END IF
	793	!
	794	za1 = 1._wp / REAL(Agrif_rhot(), wp)
[7753]	795	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
	796	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
[4486]	797	ENDIF
	798	#endif
[1502]	799	! !* write time-spliting arrays in the restart
[6140]	800	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
[508]	801	!
[9023]	802	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
	803	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
[1662]	804	!
[12377]	805	CALL iom_put( "baro_u" , puu_b(:,:,Kmm) ) ! Barotropic U Velocity
	806	CALL iom_put( "baro_v" , pvv_b(:,:,Kmm) ) ! Barotropic V Velocity
[2715]	807	!
[508]	808	END SUBROUTINE dyn_spg_ts
	809
[6140]	810
[4292]	811	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
	812	!!---------------------------------------------------------------------
	813	!! * ROUTINE ts_wgt *
	814	!!
	815	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
	816	!!----------------------------------------------------------------------
	817	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
	818	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
	819	INTEGER, INTENT(inout) :: jpit ! cycle length
[12489]	820	REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, & ! Primary weights
[4292]	821	zwgt2 ! Secondary weights
	822
	823	INTEGER :: jic, jn, ji ! temporary integers
	824	REAL(wp) :: za1, za2
	825	!!----------------------------------------------------------------------
[508]	826
[4292]	827	zwgt1(:) = 0._wp
	828	zwgt2(:) = 0._wp
	829
	830	! Set time index when averaged value is requested
	831	IF (ll_fw) THEN
[12489]	832	jic = nn_e
[4292]	833	ELSE
[12489]	834	jic = 2 * nn_e
[4292]	835	ENDIF
	836
	837	! Set primary weights:
	838	IF (ll_av) THEN
	839	! Define simple boxcar window for primary weights
[12489]	840	! (width = nn_e, centered around jic)
[4292]	841	SELECT CASE ( nn_bt_flt )
	842	CASE( 0 ) ! No averaging
	843	zwgt1(jic) = 1._wp
	844	jpit = jic
	845
[12489]	846	CASE( 1 ) ! Boxcar, width = nn_e
	847	DO jn = 1, 3*nn_e
	848	za1 = ABS(float(jn-jic))/float(nn_e)
[4292]	849	IF (za1 < 0.5_wp) THEN
	850	zwgt1(jn) = 1._wp
	851	jpit = jn
	852	ENDIF
	853	ENDDO
	854
[12489]	855	CASE( 2 ) ! Boxcar, width = 2 * nn_e
	856	DO jn = 1, 3*nn_e
	857	za1 = ABS(float(jn-jic))/float(nn_e)
[4292]	858	IF (za1 < 1._wp) THEN
	859	zwgt1(jn) = 1._wp
	860	jpit = jn
	861	ENDIF
	862	ENDDO
	863	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
	864	END SELECT
	865
	866	ELSE ! No time averaging
	867	zwgt1(jic) = 1._wp
	868	jpit = jic
	869	ENDIF
	870
	871	! Set secondary weights
	872	DO jn = 1, jpit
	873	DO ji = jn, jpit
	874	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
	875	END DO
	876	END DO
	877
	878	! Normalize weigths:
	879	za1 = 1._wp / SUM(zwgt1(1:jpit))
	880	za2 = 1._wp / SUM(zwgt2(1:jpit))
	881	DO jn = 1, jpit
	882	zwgt1(jn) = zwgt1(jn) * za1
	883	zwgt2(jn) = zwgt2(jn) * za2
	884	END DO
	885	!
	886	END SUBROUTINE ts_wgt
	887
[6140]	888
[508]	889	SUBROUTINE ts_rst( kt, cdrw )
	890	!!---------------------------------------------------------------------
	891	!! * ROUTINE ts_rst *
	892	!!
	893	!! ** Purpose : Read or write time-splitting arrays in restart file
	894	!!----------------------------------------------------------------------
[9528]	895	INTEGER , INTENT(in) :: kt ! ocean time-step
	896	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
[508]	897	!!----------------------------------------------------------------------
	898	!
[9506]	899	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
	900	! ! ---------------
[12489]	901	IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN !* Read the restart file
[13286]	902	CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	903	CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
	904	CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	905	CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
[9506]	906	IF( .NOT.ln_bt_av ) THEN
[13286]	907	CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp, ldxios = lrxios )
	908	CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	909	CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
	910	CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp, ldxios = lrxios )
	911	CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	912	CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
[9506]	913	ENDIF
[4486]	914	#if defined key_agrif
[9506]	915	! Read time integrated fluxes
	916	IF ( .NOT.Agrif_Root() ) THEN
[13286]	917	CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	918	CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
[13886]	919	ELSE
	920	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
[9506]	921	ENDIF
	922	#endif
	923	ELSE !* Start from rest
	924	IF(lwp) WRITE(numout,*)
	925	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
[13886]	926	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
	927	un_adv (:,:) = 0._wp ; vn_adv (:,:) = 0._wp ! used in the 1st interpol of agrif
	928	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
[9506]	929	#if defined key_agrif
[13886]	930	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
[9506]	931	#endif
[4486]	932	ENDIF
[9506]	933	!
	934	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
	935	! ! -------------------
	936	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
[9367]	937	IF( lwxios ) CALL iom_swap( cwxios_context )
	938	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios )
	939	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios )
	940	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios )
	941	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios )
[4292]	942	!
	943	IF (.NOT.ln_bt_av) THEN
[9367]	944	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:), ldxios = lwxios )
	945	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:), ldxios = lwxios )
	946	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:), ldxios = lwxios )
	947	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:), ldxios = lwxios )
	948	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:), ldxios = lwxios )
	949	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:), ldxios = lwxios )
[4292]	950	ENDIF
[4486]	951	#if defined key_agrif
	952	! Save time integrated fluxes
	953	IF ( .NOT.Agrif_Root() ) THEN
[9367]	954	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lwxios )
	955	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lwxios )
[4486]	956	ENDIF
	957	#endif
[9367]	958	IF( lwxios ) CALL iom_swap( cxios_context )
[4292]	959	ENDIF
	960	!
	961	END SUBROUTINE ts_rst
[2528]	962
[6140]	963
	964	SUBROUTINE dyn_spg_ts_init
[4292]	965	!!---------------------------------------------------------------------
	966	!! * ROUTINE dyn_spg_ts_init *
	967	!!
	968	!! ** Purpose : Set time splitting options
	969	!!----------------------------------------------------------------------
[6140]	970	INTEGER :: ji ,jj ! dummy loop indices
	971	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
[9019]	972	REAL(wp), DIMENSION(jpi,jpj) :: zcu
[4292]	973	!!----------------------------------------------------------------------
[4370]	974	!
[5930]	975	! Max courant number for ext. grav. waves
[4370]	976	!
[13295]	977	DO_2D( 0, 0, 0, 0 )
[12377]	978	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
	979	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
	980	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
	981	END_2D
[5930]	982	!
[13286]	983	zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) )
[10425]	984	CALL mpp_max( 'dynspg_ts', zcmax )
[2528]	985
[4370]	986	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
[12489]	987	IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
[4292]	988
[12489]	989	rDt_e = rn_Dt / REAL( nn_e , wp )
	990	zcmax = zcmax * rDt_e
[9023]	991	! Print results
[4292]	992	IF(lwp) WRITE(numout,*)
[9169]	993	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
	994	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
[5930]	995	IF( ln_bt_auto ) THEN
[12489]	996	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e '
[4370]	997	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
[4292]	998	ELSE
[12489]	999	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e
[358]	1000	ENDIF
[4292]	1001
	1002	IF(ln_bt_av) THEN
[12489]	1003	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on '
[4292]	1004	ELSE
[9169]	1005	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
[4292]	1006	ENDIF
[508]	1007	!
[4292]	1008	!
	1009	IF(ln_bt_fw) THEN
[4370]	1010	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
[4292]	1011	ELSE
[4370]	1012	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
[4292]	1013	ENDIF
	1014	!
[4486]	1015	#if defined key_agrif
	1016	! Restrict the use of Agrif to the forward case only
[9023]	1017	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
[4486]	1018	#endif
	1019	!
[4370]	1020	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
[4292]	1021	SELECT CASE ( nn_bt_flt )
[6140]	1022	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
[12489]	1023	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e'
	1024	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_e'
[9169]	1025	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
[4292]	1026	END SELECT
	1027	!
[4370]	1028	IF(lwp) WRITE(numout,*) ' '
[12489]	1029	IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e
	1030	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rDt_e
[4370]	1031	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
	1032	!
[9023]	1033	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
	1034	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
	1035	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
	1036	ENDIF
	1037	!
[6140]	1038	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
[4292]	1039	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
	1040	ENDIF
[6140]	1041	IF( zcmax>0.9_wp ) THEN
[12489]	1042	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )
[4292]	1043	ENDIF
	1044	!
[9124]	1045	! ! Allocate time-splitting arrays
	1046	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
	1047	!
	1048	! ! read restart when needed
[9506]	1049	CALL ts_rst( nit000, 'READ' )
[9124]	1050	!
[9367]	1051	IF( lwxios ) THEN
	1052	! define variables in restart file when writing with XIOS
	1053	CALL iom_set_rstw_var_active('ub2_b')
	1054	CALL iom_set_rstw_var_active('vb2_b')
	1055	CALL iom_set_rstw_var_active('un_bf')
	1056	CALL iom_set_rstw_var_active('vn_bf')
	1057	!
	1058	IF (.NOT.ln_bt_av) THEN
	1059	CALL iom_set_rstw_var_active('sshbb_e')
	1060	CALL iom_set_rstw_var_active('ubb_e')
	1061	CALL iom_set_rstw_var_active('vbb_e')
	1062	CALL iom_set_rstw_var_active('sshb_e')
	1063	CALL iom_set_rstw_var_active('ub_e')
	1064	CALL iom_set_rstw_var_active('vb_e')
	1065	ENDIF
	1066	#if defined key_agrif
	1067	! Save time integrated fluxes
	1068	IF ( .NOT.Agrif_Root() ) THEN
	1069	CALL iom_set_rstw_var_active('ub2_i_b')
	1070	CALL iom_set_rstw_var_active('vb2_i_b')
	1071	ENDIF
	1072	#endif
	1073	ENDIF
	1074	!
[4292]	1075	END SUBROUTINE dyn_spg_ts_init
[508]	1076
[11536]	1077
[12377]	1078	SUBROUTINE dyn_cor_2D_init( Kmm )
[11536]	1079	!!---------------------------------------------------------------------
[12377]	1080	!! * ROUTINE dyn_cor_2D_init *
[11536]	1081	!!
	1082	!! ** Purpose : Set time splitting options
	1083	!! Set arrays to remove/compute coriolis trend.
	1084	!! Do it once during initialization if volume is fixed, else at each long time step.
	1085	!! Note that these arrays are also used during barotropic loop. These are however frozen
	1086	!! although they should be updated in the variable volume case. Not a big approximation.
	1087	!! To remove this approximation, copy lines below inside barotropic loop
	1088	!! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
	1089	!!
	1090	!! Compute zwz = f / ( height of the water colomn )
	1091	!!----------------------------------------------------------------------
[12377]	1092	INTEGER, INTENT(in) :: Kmm ! Time index
[11536]	1093	INTEGER :: ji ,jj, jk ! dummy loop indices
	1094	REAL(wp) :: z1_ht
	1095	REAL(wp), DIMENSION(jpi,jpj) :: zhf
	1096	!!----------------------------------------------------------------------
	1097	!
	1098	SELECT CASE( nvor_scheme )
[13237]	1099	CASE( np_EEN ) != EEN scheme using e3f energy & enstrophy scheme
[11536]	1100	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
	1101	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
[13295]	1102	DO_2D( 1, 0, 1, 0 )
[12377]	1103	zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + &
	1104	& ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp
	1105	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
	1106	END_2D
[11536]	1107	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
[13295]	1108	DO_2D( 1, 0, 1, 0 )
[12377]	1109	zwz(ji,jj) = ( ht (ji ,jj+1) + ht (ji+1,jj+1) &
	1110	& + ht (ji ,jj ) + ht (ji+1,jj ) ) &
	1111	& / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) &
	1112	& + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) )
	1113	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
	1114	END_2D
[11536]	1115	END SELECT
	1116	CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
	1117	!
	1118	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[13295]	1119	DO_2D( 0, 1, 0, 1 )
[12377]	1120	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	1121	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	1122	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	1123	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	1124	END_2D
[11536]	1125	!
[13237]	1126	CASE( np_EET ) != EEN scheme using e3t energy conserving scheme
[11536]	1127	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[13295]	1128	DO_2D( 0, 1, 0, 1 )
[12377]	1129	z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
	1130	ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht
	1131	ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht
	1132	ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
	1133	ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht
	1134	END_2D
[11536]	1135	!
	1136	CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !
	1137	!
	1138	zwz(:,:) = 0._wp
	1139	zhf(:,:) = 0._wp
	1140
	1141	!!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
	1142	!!gm A priori a better value should be something like :
	1143	!!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
	1144	!!gm divided by the sum of the corresponding mask
	1145	!!gm
	1146	!!
	1147	IF( .NOT.ln_sco ) THEN
	1148
	1149	!!gm agree the JC comment : this should be done in a much clear way
	1150
	1151	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
	1152	! Set it to zero for the time being
	1153	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
	1154	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
	1155	! ENDIF
	1156	! zhf(:,:) = gdepw_0(:,:,jk+1)
	1157	!
	1158	ELSE
	1159	!
	1160	!zhf(:,:) = hbatf(:,:)
[13295]	1161	DO_2D( 1, 0, 1, 0 )
[12377]	1162	zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &
	1163	& + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &
	1164	& / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &
	1165	& + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )
	1166	END_2D
[11536]	1167	ENDIF
	1168	!
	1169	DO jj = 1, jpjm1
	1170	zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
	1171	END DO
	1172	!
	1173	DO jk = 1, jpkm1
	1174	DO jj = 1, jpjm1
[12377]	1175	zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
[11536]	1176	END DO
	1177	END DO
	1178	CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
	1179	! JC: TBC. hf should be greater than 0
[13295]	1180	DO_2D( 1, 1, 1, 1 )
[12377]	1181	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj)
	1182	END_2D
[11536]	1183	zwz(:,:) = ff_f(:,:) * zwz(:,:)
	1184	END SELECT
	1185
	1186	END SUBROUTINE dyn_cor_2d_init
	1187
	1188
	1189
[13237]	1190	SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )
[11536]	1191	!!---------------------------------------------------------------------
	1192	!! * ROUTINE dyn_cor_2d *
	1193	!!
	1194	!! ** Purpose : Compute u and v coriolis trends
	1195	!!----------------------------------------------------------------------
	1196	INTEGER :: ji ,jj ! dummy loop indices
	1197	REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - -
[13237]	1198	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhU, zhV
[11536]	1199	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd
	1200	!!----------------------------------------------------------------------
	1201	SELECT CASE( nvor_scheme )
	1202	CASE( np_ENT ) ! enstrophy conserving scheme (f-point)
[13295]	1203	DO_2D( 0, 0, 0, 0 )
[12377]	1204	z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
	1205	z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
	1206	zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu &
[13237]	1207	& * ( e1e2t(ji+1,jj)pht(ji+1,jj)ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) &
	1208	& + e1e2t(ji ,jj)pht(ji ,jj)ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) )
[12377]	1209	!
	1210	zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv &
[13237]	1211	& * ( e1e2t(ji,jj+1)pht(ji,jj+1)ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) &
	1212	& + e1e2t(ji,jj )pht(ji,jj )ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) )
[12377]	1213	END_2D
[11536]	1214	!
	1215	CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX
[13295]	1216	DO_2D( 0, 0, 0, 0 )
[12377]	1217	zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
	1218	zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
	1219	zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
	1220	zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
	1221	! energy conserving formulation for planetary vorticity term
	1222	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	1223	zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
	1224	END_2D
[11536]	1225	!
	1226	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
[13295]	1227	DO_2D( 0, 0, 0, 0 )
[12377]	1228	zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) &
	1229	& + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
	1230	zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) &
	1231	& + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
	1232	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	1233	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
	1234	END_2D
[11536]	1235	!
	1236	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
[13295]	1237	DO_2D( 0, 0, 0, 0 )
[12377]	1238	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) &
	1239	& + ftnw(ji+1,jj) * zhV(ji+1,jj ) &
	1240	& + ftse(ji,jj ) * zhV(ji ,jj-1) &
	1241	& + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
	1242	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
	1243	& + ftse(ji,jj+1) * zhU(ji ,jj+1) &
	1244	& + ftnw(ji,jj ) * zhU(ji-1,jj ) &
	1245	& + ftne(ji,jj ) * zhU(ji ,jj ) )
	1246	END_2D
[11536]	1247	!
	1248	END SELECT
	1249	!
	1250	END SUBROUTINE dyn_cor_2D
	1251
	1252
	1253	SUBROUTINE wad_tmsk( pssh, ptmsk )
	1254	!!----------------------------------------------------------------------
	1255	!! * ROUTINE wad_lmt *
	1256	!!
	1257	!! ** Purpose : set wetting & drying mask at tracer points
	1258	!! for the current barotropic sub-step
	1259	!!
	1260	!! ** Method : ???
	1261	!!
	1262	!! ** Action : ptmsk : wetting & drying t-mask
	1263	!!----------------------------------------------------------------------
	1264	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh !
	1265	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk !
	1266	!
	1267	INTEGER :: ji, jj ! dummy loop indices
	1268	!!----------------------------------------------------------------------
	1269	!
	1270	IF( ln_wd_dl_rmp ) THEN
[13295]	1271	DO_2D( 1, 1, 1, 1 )
[12377]	1272	IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
	1273	! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
	1274	ptmsk(ji,jj) = 1._wp
	1275	ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
	1276	ptmsk(ji,jj) = TANH( 50._wp( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1) )
	1277	ELSE
	1278	ptmsk(ji,jj) = 0._wp
	1279	ENDIF
	1280	END_2D
[11536]	1281	ELSE
[13295]	1282	DO_2D( 1, 1, 1, 1 )
[12377]	1283	IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp
	1284	ELSE ; ptmsk(ji,jj) = 0._wp
	1285	ENDIF
	1286	END_2D
[11536]	1287	ENDIF
	1288	!
	1289	END SUBROUTINE wad_tmsk
	1290
	1291
	1292	SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
	1293	!!----------------------------------------------------------------------
	1294	!! * ROUTINE wad_lmt *
	1295	!!
	1296	!! ** Purpose : set wetting & drying mask at tracer points
	1297	!! for the current barotropic sub-step
	1298	!!
	1299	!! ** Method : ???
	1300	!!
	1301	!! ** Action : ptmsk : wetting & drying t-mask
	1302	!!----------------------------------------------------------------------
	1303	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
	1304	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
	1305	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
	1306	!
	1307	INTEGER :: ji, jj ! dummy loop indices
	1308	!!----------------------------------------------------------------------
	1309	!
[13497]	1310	DO_2D( 1, 1, 1, 0 ) ! not jpi-column
[12377]	1311	IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj)
	1312	ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj)
	1313	ENDIF
	1314	phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
	1315	pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
	1316	END_2D
[11536]	1317	!
[13497]	1318	DO_2D( 1, 0, 1, 1 ) ! not jpj-row
[12377]	1319	IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj )
	1320	ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1)
	1321	ENDIF
	1322	phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
	1323	pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
	1324	END_2D
[11536]	1325	!
	1326	END SUBROUTINE wad_Umsk
	1327
	1328
[12377]	1329	SUBROUTINE wad_spg( pshn, zcpx, zcpy )
[11536]	1330	!!---------------------------------------------------------------------
	1331	!! * ROUTINE wad_sp *
	1332	!!
	1333	!! ** Purpose :
	1334	!!----------------------------------------------------------------------
	1335	INTEGER :: ji ,jj ! dummy loop indices
	1336	LOGICAL :: ll_tmp1, ll_tmp2
[12377]	1337	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn
[11536]	1338	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
	1339	!!----------------------------------------------------------------------
[13295]	1340	DO_2D( 0, 0, 0, 0 )
[12377]	1341	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > &
	1342	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
	1343	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji+1,jj) + ht_0(ji+1,jj) ) &
	1344	& > rn_wdmin1 + rn_wdmin2
	1345	ll_tmp2 = ( ABS( pshn(ji+1,jj) - pshn(ji ,jj)) > 1.E-12 ).AND.( &
	1346	& MAX( pshn(ji,jj) , pshn(ji+1,jj) ) > &
	1347	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
	1348	IF(ll_tmp1) THEN
	1349	zcpx(ji,jj) = 1.0_wp
	1350	ELSEIF(ll_tmp2) THEN
	1351	! no worries about pshn(ji+1,jj) - pshn(ji ,jj) = 0, it won't happen ! here
	1352	zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
	1353	& / (pshn(ji+1,jj) - pshn(ji ,jj)) )
	1354	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
	1355	ELSE
	1356	zcpx(ji,jj) = 0._wp
	1357	ENDIF
	1358	!
	1359	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji,jj+1) ) > &
	1360	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
	1361	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji,jj+1) + ht_0(ji,jj+1) ) &
	1362	& > rn_wdmin1 + rn_wdmin2
	1363	ll_tmp2 = ( ABS( pshn(ji,jj) - pshn(ji,jj+1)) > 1.E-12 ).AND.( &
	1364	& MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > &
	1365	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
	1366
	1367	IF(ll_tmp1) THEN
	1368	zcpy(ji,jj) = 1.0_wp
	1369	ELSE IF(ll_tmp2) THEN
	1370	! no worries about pshn(ji,jj+1) - pshn(ji,jj ) = 0, it won't happen ! here
	1371	zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
	1372	& / (pshn(ji,jj+1) - pshn(ji,jj )) )
	1373	zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) )
	1374	ELSE
	1375	zcpy(ji,jj) = 0._wp
	1376	ENDIF
	1377	END_2D
[11536]	1378
	1379	END SUBROUTINE wad_spg
	1380
	1381
	1382
[12377]	1383	SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
[11536]	1384	!!----------------------------------------------------------------------
	1385	!! * ROUTINE dyn_drg_init *
	1386	!!
	1387	!! ** Purpose : - add the baroclinic top/bottom drag contribution to
	1388	!! the baroclinic part of the barotropic RHS
	1389	!! - compute the barotropic drag coefficients
	1390	!!
	1391	!! ** Method : computation done over the INNER domain only
	1392	!!----------------------------------------------------------------------
[12377]	1393	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
	1394	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation
	1395	REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels
	1396	REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS
	1397	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients
[11536]	1398	!
	1399	INTEGER :: ji, jj ! dummy loop indices
	1400	INTEGER :: ikbu, ikbv, iktu, iktv
	1401	REAL(wp) :: zztmp
	1402	REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i
	1403	!!----------------------------------------------------------------------
	1404	!
	1405	! !== Set the barotropic drag coef. ==!
	1406	!
[13472]	1407	IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities)
[11536]	1408
[13295]	1409	DO_2D( 0, 0, 0, 0 )
[12377]	1410	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
	1411	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
	1412	END_2D
[11536]	1413	ELSE ! bottom friction only
[13295]	1414	DO_2D( 0, 0, 0, 0 )
[12377]	1415	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
	1416	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
	1417	END_2D
[11536]	1418	ENDIF
	1419	!
	1420	! !== BOTTOM stress contribution from baroclinic velocities ==!
	1421	!
	1422	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
	1423
[13295]	1424	DO_2D( 0, 0, 0, 0 )
[12377]	1425	ikbu = mbku(ji,jj)
	1426	ikbv = mbkv(ji,jj)
	1427	zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
	1428	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
	1429	END_2D
[11536]	1430	ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
	1431
[13295]	1432	DO_2D( 0, 0, 0, 0 )
[12377]	1433	ikbu = mbku(ji,jj)
	1434	ikbv = mbkv(ji,jj)
	1435	zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
	1436	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
	1437	END_2D
[11536]	1438	ENDIF
	1439	!
	1440	IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please !
[12489]	1441	zztmp = -1._wp / rDt_e
[13295]	1442	DO_2D( 0, 0, 0, 0 )
[12377]	1443	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
	1444	& r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp )
	1445	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &
	1446	& r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
	1447	END_2D
[11536]	1448	ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
	1449
[13295]	1450	DO_2D( 0, 0, 0, 0 )
[12377]	1451	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)
	1452	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)
	1453	END_2D
[11536]	1454	END IF
	1455	!
	1456	! !== TOP stress contribution from baroclinic velocities ==! (no W/D case)
	1457	!
[13472]	1458	IF( ln_isfcav.OR.ln_drgice_imp ) THEN
[11536]	1459	!
	1460	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
	1461
[13295]	1462	DO_2D( 0, 0, 0, 0 )
[12377]	1463	iktu = miku(ji,jj)
	1464	iktv = mikv(ji,jj)
	1465	zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
	1466	zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
	1467	END_2D
[11536]	1468	ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
	1469
[13295]	1470	DO_2D( 0, 0, 0, 0 )
[12377]	1471	iktu = miku(ji,jj)
	1472	iktv = mikv(ji,jj)
	1473	zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
	1474	zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
	1475	END_2D
[11536]	1476	ENDIF
	1477	!
	1478	! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
	1479
[13295]	1480	DO_2D( 0, 0, 0, 0 )
[12377]	1481	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)
	1482	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)
	1483	END_2D
[11536]	1484	!
	1485	ENDIF
	1486	!
	1487	END SUBROUTINE dyn_drg_init
	1488
	1489	SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in
	1490	& za0, za1, za2, za3 ) ! ==> out
	1491	!!----------------------------------------------------------------------
	1492	INTEGER ,INTENT(in ) :: jn ! index of sub time step
	1493	LOGICAL ,INTENT(in ) :: ll_init !
	1494	REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient
	1495	!
	1496	REAL(wp) :: zepsilon, zgamma ! - -
	1497	!!----------------------------------------------------------------------
	1498	! ! set Half-step back interpolation coefficient
	1499	IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward
	1500	za0 = 1._wp
	1501	za1 = 0._wp
	1502	za2 = 0._wp
	1503	za3 = 0._wp
	1504	ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
	1505	za0 = 1.0833333333333_wp ! za0 = 1-gam-eps
	1506	za1 =-0.1666666666666_wp ! za1 = gam
	1507	za2 = 0.0833333333333_wp ! za2 = eps
	1508	za3 = 0._wp
	1509	ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
	1510	IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion
	1511	za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps
	1512	za1 = 0.285_wp ! za1 = 1/2 - 2gam - 3eps
	1513	za2 = 0.088_wp ! za2 = gam
	1514	za3 = 0.013_wp ! za3 = eps
	1515	ELSE ! no time diffusion
	1516	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
	1517	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
	1518	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
	1519	za1 = 1._wp - za0 - zgamma - zepsilon
	1520	za2 = zgamma
	1521	za3 = zepsilon
	1522	ENDIF
	1523	ENDIF
	1524	END SUBROUTINE ts_bck_interp
	1525
	1526
[358]	1527	!!======================================================================
	1528	END MODULE dynspg_ts

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

Download in other formats: