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_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynspg_ts.F90 @ 14018

Last change on this file since 14018 was 14018, checked in by techene, 3 years ago
#2385 branch updated with trunk 13970
Property svn:keywords set to `Id`
File size: 77.7 KB

Rev	Line
[358]	1	MODULE dynspg_ts
[9023]	2
[12377]	3	!! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out !
[9023]	4
[358]	5	!!======================================================================
[6140]	6	!! * MODULE dynspg_ts *
	7	!! Ocean dynamics: surface pressure gradient trend, split-explicit scheme
	8	!!======================================================================
[1502]	9	!! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code
	10	!! - ! 2005-11 (V. Garnier, G. Madec) optimization
	11	!! - ! 2006-08 (S. Masson) distributed restart using iom
	12	!! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines
	13	!! - ! 2008-01 (R. Benshila) change averaging method
	14	!! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
[2528]	15	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
[2724]	16	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
[4292]	17	!! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping
	18	!! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility
[5930]	19	!! 3.7 ! 2015-11 (J. Chanut) free surface simplification
[7646]	20	!! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence
[9019]	21	!! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90)
[2724]	22	!!---------------------------------------------------------------------
[6140]	23
[358]	24	!!----------------------------------------------------------------------
[6140]	25	!! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme
	26	!! dyn_spg_ts_init: initialisation of the time-splitting scheme
	27	!! ts_wgt : set time-splitting weights for temporal averaging (or not)
	28	!! ts_rst : read/write time-splitting fields in restart file
[358]	29	!!----------------------------------------------------------------------
	30	USE oce ! ocean dynamics and tracers
	31	USE dom_oce ! ocean space and time domain
[888]	32	USE sbc_oce ! surface boundary condition: ocean
[12377]	33	USE isf_oce ! ice shelf variable (fwfisf)
[9019]	34	USE zdf_oce ! vertical physics: variables
	35	USE zdfdrg ! vertical physics: top/bottom drag coef.
[6140]	36	USE sbcapr ! surface boundary condition: atmospheric pressure
	37	USE dynadv , ONLY: ln_dynadv_vec
[9528]	38	USE dynvor ! vortivity scheme indicators
[358]	39	USE phycst ! physical constants
	40	USE dynvor ! vorticity term
[6152]	41	USE wet_dry ! wetting/drying flux limter
[7646]	42	USE bdy_oce ! open boundary
[10481]	43	USE bdyvol ! open boundary volume conservation
[5930]	44	USE bdytides ! open boundary condition data
[3294]	45	USE bdydyn2d ! open boundary conditions on barotropic variables
[12377]	46	USE tide_mod !
[7646]	47	USE sbcwave ! surface wave
[9019]	48	#if defined key_agrif
[9570]	49	USE agrif_oce_interp ! agrif
[9124]	50	USE agrif_oce
[9019]	51	#endif
	52	#if defined key_asminc
	53	USE asminc ! Assimilation increment
	54	#endif
[6140]	55	!
	56	USE in_out_manager ! I/O manager
[358]	57	USE lib_mpp ! distributed memory computing library
	58	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	59	USE prtctl ! Print control
[2715]	60	USE iom ! IOM library
[4292]	61	USE restart ! only for lrst_oce
[358]	62
[11536]	63	USE iom ! to remove
	64
[358]	65	IMPLICIT NONE
	66	PRIVATE
	67
[9124]	68	PUBLIC dyn_spg_ts ! called by dyn_spg
	69	PUBLIC dyn_spg_ts_init ! - - dyn_spg_init
[358]	70
[9019]	71	!! Time filtered arrays at baroclinic time step:
	72	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step
[9124]	73	!
[12489]	74	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_e <= 2.5 nn_e
	75	REAL(wp),SAVE :: rDt_e ! Barotropic time step
[9019]	76	!
	77	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields
[9124]	78	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points
	79	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
	80	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
[4292]	81
[9043]	82	REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios
	83	REAL(wp) :: r1_8 = 0.125_wp !
	84	REAL(wp) :: r1_4 = 0.25_wp !
	85	REAL(wp) :: r1_2 = 0.5_wp !
[508]	86
[358]	87	!! * Substitutions
[12377]	88	# include "do_loop_substitute.h90"
[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
[13998]	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	!
[13998]	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	!
[13427]	308	! != Add wind forcing =!
	309	! ! ------------------ !
	310	IF( ln_bt_fw ) THEN
[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)
[13427]	388	zhup2_e(:,:) = hu_0(:,:)
	389	zhvp2_e(:,:) = hv_0(:,:)
	390	zhtp2_e(:,:) = ht_0(:,:)
[11536]	391	ENDIF
	392	!
[13427]	393	IF( ln_bt_fw ) THEN ! FORWARD integration: start from NOW fields
	394	sshn_e(:,:) = pssh (:,:,Kmm)
[12377]	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
[13427]	403	sshn_e(:,:) = pssh (:,:,Kbb)
[12377]	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:
[13427]	414	puu_b (:,:,Kaa) = 0._wp ! After barotropic velocities (or transport if flux form)
	415	pvv_b (:,:,Kaa) = 0._wp
[12377]	416	pssh (:,:,Kaa) = 0._wp ! Sum for after averaged sea level
[13427]	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)
[13427]	477	#if defined key_qcoTest_FluxForm
	478	! ! 'key_qcoTest_FluxForm' : simple ssh average
[13998]	479	DO_2D( 1, 1, 1, 0 ) ! not jpi-column
[13427]	480	zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj ) ) * ssumask(ji,jj)
	481	END_2D
	482	DO_2D( 1, 0, 1, 1 )
	483	zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji ,jj+1) ) * ssvmask(ji,jj)
	484	END_2D
	485	#else
	486	! ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
[13998]	487	DO_2D( 1, 1, 1, 0 ) ! not jpi-column
[12377]	488	zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) &
	489	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
	490	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
	491	END_2D
[13998]	492	DO_2D( 1, 0, 1, 1 ) ! not jpj-row
[12377]	493	zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) &
	494	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
	495	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
	496	END_2D
[13427]	497	#endif
[4292]	498	!
	499	ENDIF
	500	!
[11536]	501	! !== after SSH ==! (jn+1)
	502	!
	503	! ! update (ua_e,va_e) to enforce volume conservation at open boundaries
	504	! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
[10481]	505	IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
[12377]	506	!
[11536]	507	! ! resulting flux at mid-step (not over the full domain)
	508	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
	509	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]	510	!
	511	#if defined key_agrif
[6140]	512	! Set fluxes during predictor step to ensure volume conservation
[12377]	513	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV )
[4486]	514	#endif
[12489]	515	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]	516
[11536]	517	IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where
	518	! ! the direction of the flow is from dry cells
	519	CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V
[9528]	520	!
	521	ENDIF
[11536]	522	!
	523	!
	524	! Compute Sea Level at step jit+1
	525	!-- m+1 m m+1/2 --!
	526	!-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --!
	527	!-------------------------------------------------------------------------!
[13295]	528	DO_2D( 0, 0, 0, 0 )
[12377]	529	zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj)
[12489]	530	ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)
[12377]	531	END_2D
[11536]	532	!
[13289]	533	CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )
	534	!
[13216]	535	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
	536	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
	537	#if defined key_agrif
	538	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
	539	#endif
[11536]	540	!
	541	! ! Sum over sub-time-steps to compute advective velocities
	542	za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5
	543	un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
	544	vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
	545	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
	546	IF ( ln_wd_dl_bc ) THEN
	547	zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column
	548	zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row
	549	END IF
	550	!
[4292]	551	!
	552	! Sea Surface Height at u-,v-points (vvl case only)
[13427]	553	IF( .NOT.ln_linssh ) THEN
	554	#if defined key_qcoTest_FluxForm
	555	! ! 'key_qcoTest_FluxForm' : simple ssh average
	556	DO_2D( 1, 1, 1, 0 )
	557	zsshu_a(ji,jj) = r1_2 * ( ssha_e(ji,jj) + ssha_e(ji+1,jj ) ) * ssumask(ji,jj)
	558	END_2D
	559	DO_2D( 1, 0, 1, 1 )
	560	zsshv_a(ji,jj) = r1_2 * ( ssha_e(ji,jj) + ssha_e(ji ,jj+1) ) * ssvmask(ji,jj)
	561	END_2D
	562	#else
[13295]	563	DO_2D( 0, 0, 0, 0 )
[13427]	564	zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
	565	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) * ssumask(ji,jj)
	566	zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
	567	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) * ssvmask(ji,jj)
[12377]	568	END_2D
[13427]	569	#endif
	570	ENDIF
[11536]	571	!
	572	! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
	573	!-- m+1/2 m+1 m m-1 m-2 --!
	574	!-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --!
	575	!------------------------------------------------------------------------------------------!
	576	CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation
[9528]	577	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
	578	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
[1502]	579	!
[11536]	580	! ! Surface pressure gradient
	581	zldg = ( 1._wp - rn_scal_load ) * grav ! local factor
[13295]	582	DO_2D( 0, 0, 0, 0 )
[12377]	583	zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
	584	zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
	585	END_2D
[11536]	586	IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient
	587	CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters
	588	zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
	589	zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
[4292]	590	ENDIF
	591	!
	592	! Add Coriolis trend:
[6140]	593	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
[4292]	594	! at each time step. We however keep them constant here for optimization.
[11536]	595	! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
[13237]	596	CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )
[4292]	597	!
	598	! Add tidal astronomical forcing if defined
[7646]	599	IF ( ln_tide .AND. ln_tide_pot ) THEN
[13295]	600	DO_2D( 0, 0, 0, 0 )
[12377]	601	zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
	602	zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
	603	END_2D
[4292]	604	ENDIF
	605	!
[9023]	606	! Add bottom stresses:
	607	!jth do implicitly instead
	608	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
[13295]	609	DO_2D( 0, 0, 0, 0 )
[12377]	610	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
	611	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
	612	END_2D
[11536]	613	ENDIF
[4292]	614	!
	615	! Set next velocities:
[11536]	616	! Compute barotropic speeds at step jit+1 (h : total height of the water colomn)
	617	!-- VECTOR FORM
	618	!-- m+1 m / m+1/2 \ --!
	619	!-- u = u + delta_t' * \ (1-r)g grad_x( ssh') - f * k vect u + frc / --!
	620	!-- --!
	621	!-- FLUX FORM --!
	622	!-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --!
	623	!-- u = m+1 \| h * u + delta_t' * \ h * (1-r)g grad_x( ssh') - h * f * k vect u + h * frc / \| --!
	624	!-- h \ / --!
	625	!------------------------------------------------------------------------------------------------------------------------!
[9023]	626	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
[13295]	627	DO_2D( 0, 0, 0, 0 )
[12377]	628	ua_e(ji,jj) = ( un_e(ji,jj) &
[12489]	629	& + rDt_e * ( zu_spg(ji,jj) &
[12377]	630	& + zu_trd(ji,jj) &
	631	& + zu_frc(ji,jj) ) &
	632	& ) * ssumask(ji,jj)
[358]	633
[12377]	634	va_e(ji,jj) = ( vn_e(ji,jj) &
[12489]	635	& + rDt_e * ( zv_spg(ji,jj) &
[12377]	636	& + zv_trd(ji,jj) &
	637	& + zv_frc(ji,jj) ) &
	638	& ) * ssvmask(ji,jj)
	639	END_2D
[6140]	640	!
[9023]	641	ELSE !* Flux form
[13295]	642	DO_2D( 0, 0, 0, 0 )
[12377]	643	! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2
	644	! ! backward interpolated depth used in spg terms at jn+1/2
[13427]	645	#if defined key_qcoTest_FluxForm
	646	! ! 'key_qcoTest_FluxForm' : simple ssh average
	647	zhu_bck = hu_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj ) ) * ssumask(ji,jj)
	648	zhv_bck = hv_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji ,jj+1) ) * ssvmask(ji,jj)
	649	#else
[12377]	650	zhu_bck = hu_0(ji,jj) + r1_2r1_e1e2u(ji,jj) ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
	651	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj)
	652	zhv_bck = hv_0(ji,jj) + r1_2r1_e1e2v(ji,jj) ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
	653	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj)
[13427]	654	#endif
[12377]	655	! ! inverse depth at jn+1
	656	z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
	657	z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
	658	!
	659	ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) &
[12489]	660	& + rDt_e * ( zhu_bck * zu_spg (ji,jj) & !
[12377]	661	& + zhup2_e(ji,jj) * zu_trd (ji,jj) & !
	662	& + hu(ji,jj,Kmm) * zu_frc (ji,jj) ) ) * z1_hu
	663	!
	664	va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) &
[12489]	665	& + rDt_e * ( zhv_bck * zv_spg (ji,jj) & !
[12377]	666	& + zhvp2_e(ji,jj) * zv_trd (ji,jj) & !
	667	& + hv(ji,jj,Kmm) * zv_frc (ji,jj) ) ) * z1_hv
	668	END_2D
[4292]	669	ENDIF
[10272]	670	!jth implicit bottom friction:
	671	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
[13295]	672	DO_2D( 0, 0, 0, 0 )
[13427]	673	ua_e(ji,jj) = ua_e(ji,jj) / ( 1._wp - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj) )
	674	va_e(ji,jj) = va_e(ji,jj) / ( 1._wp - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj) )
[12377]	675	END_2D
[10272]	676	ENDIF
[11536]	677
[13216]	678	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
[13427]	679	hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
	680	hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) )
	681	hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
	682	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]	683	ENDIF
	684	!
	685	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
[11536]	686	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp &
[12072]	687	& , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp &
	688	& , hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp )
[11536]	689	ELSE
	690	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )
[1438]	691	ENDIF
[13289]	692	! ! open boundaries
	693	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
	694	#if defined key_agrif
	695	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
	696	#endif
[4292]	697	! !* Swap
	698	! ! ----
[7753]	699	ubb_e (:,:) = ub_e (:,:)
	700	ub_e (:,:) = un_e (:,:)
	701	un_e (:,:) = ua_e (:,:)
	702	!
	703	vbb_e (:,:) = vb_e (:,:)
	704	vb_e (:,:) = vn_e (:,:)
	705	vn_e (:,:) = va_e (:,:)
	706	!
	707	sshbb_e(:,:) = sshb_e(:,:)
	708	sshb_e (:,:) = sshn_e(:,:)
	709	sshn_e (:,:) = ssha_e(:,:)
[4292]	710
	711	! !* Sum over whole bt loop
	712	! ! ----------------------
	713	za1 = wgtbtp1(jn)
[6140]	714	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
[12377]	715	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:)
	716	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:)
[9023]	717	ELSE ! Sum transports
	718	IF ( .NOT.ln_wd_dl ) THEN
[12377]	719	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:)
	720	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:)
[9023]	721	ELSE
[12377]	722	puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
	723	pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
[9023]	724	END IF
[4292]	725	ENDIF
[9023]	726	! ! Sum sea level
[12377]	727	pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
[9023]	728
[358]	729	! ! ==================== !
	730	END DO ! end loop !
	731	! ! ==================== !
[1438]	732	! -----------------------------------------------------------------------------
[1502]	733	! Phase 3. update the general trend with the barotropic trend
[1438]	734	! -----------------------------------------------------------------------------
[1502]	735	!
[4292]	736	! Set advection velocity correction:
[9023]	737	IF (ln_bt_fw) THEN
[12489]	738	IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
[13295]	739	DO_2D( 1, 1, 1, 1 )
[12377]	740	zun_save = un_adv(ji,jj)
	741	zvn_save = vn_adv(ji,jj)
	742	! ! apply the previously computed correction
[12489]	743	un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) )
	744	vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) )
[12377]	745	! ! Update corrective fluxes for next time step
[12489]	746	un_bf(ji,jj) = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
	747	vn_bf(ji,jj) = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
[12377]	748	! ! Save integrated transport for next computation
	749	ub2_b(ji,jj) = zun_save
	750	vb2_b(ji,jj) = zvn_save
	751	END_2D
[9023]	752	ELSE
[11536]	753	un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero
	754	vn_bf(:,:) = 0._wp
	755	ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation
	756	vb2_b(:,:) = vn_adv(:,:)
	757	END IF
[4292]	758	ENDIF
[9023]	759
	760
[4292]	761	!
	762	! Update barotropic trend:
[6140]	763	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
[4292]	764	DO jk=1,jpkm1
[12489]	765	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b
	766	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b
[4292]	767	END DO
	768	ELSE
[12377]	769	! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
[13427]	770	#if defined key_qcoTest_FluxForm
	771	! ! 'key_qcoTest_FluxForm' : simple ssh average
[13295]	772	DO_2D( 1, 0, 1, 0 )
[13427]	773	zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj ,Kaa) ) * ssumask(ji,jj)
	774	zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji ,jj+1,Kaa) ) * ssvmask(ji,jj)
[12377]	775	END_2D
[13427]	776	#else
	777	DO_2D( 1, 0, 1, 0 )
	778	zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) &
	779	& + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj)
	780	zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) &
	781	& + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj)
	782	END_2D
	783	#endif
[10425]	784	CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
[5930]	785	!
[4292]	786	DO jk=1,jpkm1
[13427]	787	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) &
	788	& * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b
	789	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) &
	790	& * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b
[4292]	791	END DO
	792	! Save barotropic velocities not transport:
[12377]	793	puu_b(:,:,Kaa) = puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
	794	pvv_b(:,:,Kaa) = pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
[4292]	795	ENDIF
[9023]	796
	797
	798	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
[4292]	799	DO jk = 1, jpkm1
[12377]	800	puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) umask(:,:,jk)
	801	pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) vmask(:,:,jk)
[358]	802	END DO
[9023]	803
	804	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
	805	DO jk = 1, jpkm1
[12377]	806	puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
	807	& + zuwdav2(:,:)(puu(:,:,jk,Kmm) - un_adv(:,:)r1_hu(:,:,Kmm)) ) * umask(:,:,jk)
	808	pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) &
	809	& + zvwdav2(:,:)(pvv(:,:,jk,Kmm) - vn_adv(:,:)r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)
[9023]	810	END DO
	811	END IF
	812
	813
[12377]	814	CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current
	815	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current
[1502]	816	!
[4486]	817	#if defined key_agrif
	818	! Save time integrated fluxes during child grid integration
[5656]	819	! (used to update coarse grid transports at next time step)
[4486]	820	!
[6140]	821	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
	822	IF( Agrif_NbStepint() == 0 ) THEN
[7753]	823	ub2_i_b(:,:) = 0._wp
	824	vb2_i_b(:,:) = 0._wp
[4486]	825	END IF
	826	!
	827	za1 = 1._wp / REAL(Agrif_rhot(), wp)
[7753]	828	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
	829	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
[4486]	830	ENDIF
	831	#endif
[1502]	832	! !* write time-spliting arrays in the restart
[6140]	833	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
[508]	834	!
[9023]	835	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
	836	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
[1662]	837	!
[12377]	838	CALL iom_put( "baro_u" , puu_b(:,:,Kmm) ) ! Barotropic U Velocity
	839	CALL iom_put( "baro_v" , pvv_b(:,:,Kmm) ) ! Barotropic V Velocity
[2715]	840	!
[508]	841	END SUBROUTINE dyn_spg_ts
	842
[6140]	843
[4292]	844	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
	845	!!---------------------------------------------------------------------
	846	!! * ROUTINE ts_wgt *
	847	!!
	848	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
	849	!!----------------------------------------------------------------------
	850	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
	851	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
	852	INTEGER, INTENT(inout) :: jpit ! cycle length
[12489]	853	REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, & ! Primary weights
[4292]	854	zwgt2 ! Secondary weights
	855
	856	INTEGER :: jic, jn, ji ! temporary integers
	857	REAL(wp) :: za1, za2
	858	!!----------------------------------------------------------------------
[508]	859
[4292]	860	zwgt1(:) = 0._wp
	861	zwgt2(:) = 0._wp
	862
	863	! Set time index when averaged value is requested
	864	IF (ll_fw) THEN
[12489]	865	jic = nn_e
[4292]	866	ELSE
[12489]	867	jic = 2 * nn_e
[4292]	868	ENDIF
	869
	870	! Set primary weights:
	871	IF (ll_av) THEN
	872	! Define simple boxcar window for primary weights
[12489]	873	! (width = nn_e, centered around jic)
[4292]	874	SELECT CASE ( nn_bt_flt )
	875	CASE( 0 ) ! No averaging
	876	zwgt1(jic) = 1._wp
	877	jpit = jic
	878
[12489]	879	CASE( 1 ) ! Boxcar, width = nn_e
	880	DO jn = 1, 3*nn_e
	881	za1 = ABS(float(jn-jic))/float(nn_e)
[4292]	882	IF (za1 < 0.5_wp) THEN
	883	zwgt1(jn) = 1._wp
	884	jpit = jn
	885	ENDIF
	886	ENDDO
	887
[12489]	888	CASE( 2 ) ! Boxcar, width = 2 * nn_e
	889	DO jn = 1, 3*nn_e
	890	za1 = ABS(float(jn-jic))/float(nn_e)
[4292]	891	IF (za1 < 1._wp) THEN
	892	zwgt1(jn) = 1._wp
	893	jpit = jn
	894	ENDIF
	895	ENDDO
	896	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
	897	END SELECT
	898
	899	ELSE ! No time averaging
	900	zwgt1(jic) = 1._wp
	901	jpit = jic
	902	ENDIF
	903
	904	! Set secondary weights
	905	DO jn = 1, jpit
	906	DO ji = jn, jpit
	907	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
	908	END DO
	909	END DO
	910
	911	! Normalize weigths:
	912	za1 = 1._wp / SUM(zwgt1(1:jpit))
	913	za2 = 1._wp / SUM(zwgt2(1:jpit))
	914	DO jn = 1, jpit
	915	zwgt1(jn) = zwgt1(jn) * za1
	916	zwgt2(jn) = zwgt2(jn) * za2
	917	END DO
	918	!
	919	END SUBROUTINE ts_wgt
	920
[6140]	921
[508]	922	SUBROUTINE ts_rst( kt, cdrw )
	923	!!---------------------------------------------------------------------
	924	!! * ROUTINE ts_rst *
	925	!!
	926	!! ** Purpose : Read or write time-splitting arrays in restart file
	927	!!----------------------------------------------------------------------
[9528]	928	INTEGER , INTENT(in) :: kt ! ocean time-step
	929	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
[508]	930	!!----------------------------------------------------------------------
	931	!
[9506]	932	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
	933	! ! ---------------
[13895]	934	IF( ln_rstart .AND. ln_bt_fw .AND. .NOT.l_1st_euler ) THEN !* Read the restart file
[14018]	935	CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp )
	936	CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp )
	937	CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp )
	938	CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp )
[9506]	939	IF( .NOT.ln_bt_av ) THEN
[14018]	940	CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp )
	941	CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp )
	942	CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp )
	943	CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp )
	944	CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp )
	945	CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp )
[9506]	946	ENDIF
[4486]	947	#if defined key_agrif
[9506]	948	! Read time integrated fluxes
	949	IF ( .NOT.Agrif_Root() ) THEN
[14018]	950	CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp )
	951	CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp )
[13998]	952	ELSE
	953	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
[9506]	954	ENDIF
	955	#endif
	956	ELSE !* Start from rest
	957	IF(lwp) WRITE(numout,*)
	958	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
[13998]	959	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
	960	un_adv (:,:) = 0._wp ; vn_adv (:,:) = 0._wp ! used in the 1st interpol of agrif
	961	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
[9506]	962	#if defined key_agrif
[13998]	963	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
[9506]	964	#endif
[4486]	965	ENDIF
[9506]	966	!
	967	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
	968	! ! -------------------
	969	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
[14018]	970	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
	971	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
	972	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:) )
	973	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:) )
[4292]	974	!
	975	IF (.NOT.ln_bt_av) THEN
[14018]	976	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
	977	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
	978	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
	979	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
	980	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
	981	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
[4292]	982	ENDIF
[4486]	983	#if defined key_agrif
	984	! Save time integrated fluxes
	985	IF ( .NOT.Agrif_Root() ) THEN
[14018]	986	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
	987	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
[4486]	988	ENDIF
	989	#endif
[4292]	990	ENDIF
	991	!
	992	END SUBROUTINE ts_rst
[2528]	993
[6140]	994
	995	SUBROUTINE dyn_spg_ts_init
[4292]	996	!!---------------------------------------------------------------------
	997	!! * ROUTINE dyn_spg_ts_init *
	998	!!
	999	!! ** Purpose : Set time splitting options
	1000	!!----------------------------------------------------------------------
[6140]	1001	INTEGER :: ji ,jj ! dummy loop indices
	1002	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
[9019]	1003	REAL(wp), DIMENSION(jpi,jpj) :: zcu
[4292]	1004	!!----------------------------------------------------------------------
[4370]	1005	!
[5930]	1006	! Max courant number for ext. grav. waves
[4370]	1007	!
[13295]	1008	DO_2D( 0, 0, 0, 0 )
[12377]	1009	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
	1010	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
	1011	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
	1012	END_2D
[5930]	1013	!
[13286]	1014	zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) )
[10425]	1015	CALL mpp_max( 'dynspg_ts', zcmax )
[2528]	1016
[4370]	1017	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
[12489]	1018	IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
[4292]	1019
[12489]	1020	rDt_e = rn_Dt / REAL( nn_e , wp )
	1021	zcmax = zcmax * rDt_e
[9023]	1022	! Print results
[4292]	1023	IF(lwp) WRITE(numout,*)
[9169]	1024	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
	1025	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
[5930]	1026	IF( ln_bt_auto ) THEN
[12489]	1027	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e '
[4370]	1028	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
[4292]	1029	ELSE
[12489]	1030	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e
[358]	1031	ENDIF
[4292]	1032
	1033	IF(ln_bt_av) THEN
[12489]	1034	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on '
[4292]	1035	ELSE
[9169]	1036	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
[4292]	1037	ENDIF
[508]	1038	!
[4292]	1039	!
	1040	IF(ln_bt_fw) THEN
[4370]	1041	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
[4292]	1042	ELSE
[4370]	1043	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
[4292]	1044	ENDIF
	1045	!
[4486]	1046	#if defined key_agrif
	1047	! Restrict the use of Agrif to the forward case only
[9023]	1048	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
[4486]	1049	#endif
	1050	!
[4370]	1051	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
[4292]	1052	SELECT CASE ( nn_bt_flt )
[6140]	1053	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
[12489]	1054	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e'
	1055	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_e'
[9169]	1056	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
[4292]	1057	END SELECT
	1058	!
[4370]	1059	IF(lwp) WRITE(numout,*) ' '
[12489]	1060	IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e
	1061	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rDt_e
[4370]	1062	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
	1063	!
[9023]	1064	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
	1065	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
	1066	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
	1067	ENDIF
	1068	!
[6140]	1069	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
[4292]	1070	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
	1071	ENDIF
[6140]	1072	IF( zcmax>0.9_wp ) THEN
[12489]	1073	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )
[4292]	1074	ENDIF
	1075	!
[9124]	1076	! ! Allocate time-splitting arrays
	1077	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
	1078	!
	1079	! ! read restart when needed
[9506]	1080	CALL ts_rst( nit000, 'READ' )
[9124]	1081	!
[4292]	1082	END SUBROUTINE dyn_spg_ts_init
[508]	1083
[11536]	1084
[12377]	1085	SUBROUTINE dyn_cor_2D_init( Kmm )
[11536]	1086	!!---------------------------------------------------------------------
[12377]	1087	!! * ROUTINE dyn_cor_2D_init *
[11536]	1088	!!
	1089	!! ** Purpose : Set time splitting options
	1090	!! Set arrays to remove/compute coriolis trend.
	1091	!! Do it once during initialization if volume is fixed, else at each long time step.
	1092	!! Note that these arrays are also used during barotropic loop. These are however frozen
	1093	!! although they should be updated in the variable volume case. Not a big approximation.
	1094	!! To remove this approximation, copy lines below inside barotropic loop
[13696]	1095	!! and update depths at T- points (ht) at each barotropic time step
[11536]	1096	!!
	1097	!! Compute zwz = f / ( height of the water colomn )
	1098	!!----------------------------------------------------------------------
[12377]	1099	INTEGER, INTENT(in) :: Kmm ! Time index
[11536]	1100	INTEGER :: ji ,jj, jk ! dummy loop indices
	1101	REAL(wp) :: z1_ht
	1102	!!----------------------------------------------------------------------
	1103	!
	1104	SELECT CASE( nvor_scheme )
[13696]	1105	CASE( np_EEN, np_ENE, np_ENS , np_MIX ) != schemes using the same e3f definition
	1106	SELECT CASE( nn_e3f_typ ) !* ff_f/e3 at F-point
[11536]	1107	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
[13696]	1108	DO_2D( 0, 0, 0, 0 )
[13427]	1109	zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) &
	1110	& + ht(ji,jj ) + ht(ji+1,jj ) ) * 0.25_wp
[12377]	1111	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
	1112	END_2D
[11536]	1113	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
[13696]	1114	DO_2D( 0, 0, 0, 0 )
[13427]	1115	zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) &
	1116	& + ht(ji,jj ) + ht(ji+1,jj ) ) &
	1117	& / ( MAX(ssmask(ji,jj+1) + ssmask(ji+1,jj+1) &
	1118	& + ssmask(ji,jj ) + ssmask(ji+1,jj ) , 1._wp ) )
[12377]	1119	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
	1120	END_2D
[11536]	1121	END SELECT
	1122	CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
[13696]	1123	END SELECT
	1124	!
	1125	SELECT CASE( nvor_scheme )
	1126	CASE( np_EEN )
[11536]	1127	!
[13427]	1128	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[13295]	1129	DO_2D( 0, 1, 0, 1 )
[12377]	1130	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	1131	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	1132	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	1133	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	1134	END_2D
[11536]	1135	!
[13696]	1136	CASE( np_EET ) != EEN scheme using e3t energy conserving scheme
[13427]	1137	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[13295]	1138	DO_2D( 0, 1, 0, 1 )
[12377]	1139	z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
	1140	ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht
	1141	ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht
	1142	ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
	1143	ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht
	1144	END_2D
[11536]	1145	!
	1146	END SELECT
	1147
	1148	END SUBROUTINE dyn_cor_2d_init
	1149
	1150
[13237]	1151	SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )
[11536]	1152	!!---------------------------------------------------------------------
	1153	!! * ROUTINE dyn_cor_2d *
	1154	!!
	1155	!! ** Purpose : Compute u and v coriolis trends
	1156	!!----------------------------------------------------------------------
	1157	INTEGER :: ji ,jj ! dummy loop indices
	1158	REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - -
[13237]	1159	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhU, zhV
[11536]	1160	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd
	1161	!!----------------------------------------------------------------------
	1162	SELECT CASE( nvor_scheme )
	1163	CASE( np_ENT ) ! enstrophy conserving scheme (f-point)
[13295]	1164	DO_2D( 0, 0, 0, 0 )
[12377]	1165	z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
	1166	z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
	1167	zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu &
[13237]	1168	& * ( e1e2t(ji+1,jj)pht(ji+1,jj)ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) &
	1169	& + e1e2t(ji ,jj)pht(ji ,jj)ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) )
[12377]	1170	!
	1171	zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv &
[13237]	1172	& * ( e1e2t(ji,jj+1)pht(ji,jj+1)ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) &
	1173	& + e1e2t(ji,jj )pht(ji,jj )ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) )
[12377]	1174	END_2D
[11536]	1175	!
	1176	CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX
[13295]	1177	DO_2D( 0, 0, 0, 0 )
[12377]	1178	zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
	1179	zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
	1180	zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
	1181	zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
	1182	! energy conserving formulation for planetary vorticity term
	1183	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	1184	zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
	1185	END_2D
[11536]	1186	!
	1187	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
[13295]	1188	DO_2D( 0, 0, 0, 0 )
[12377]	1189	zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) &
	1190	& + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
	1191	zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) &
	1192	& + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
	1193	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	1194	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
	1195	END_2D
[11536]	1196	!
	1197	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
[13295]	1198	DO_2D( 0, 0, 0, 0 )
[12377]	1199	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) &
	1200	& + ftnw(ji+1,jj) * zhV(ji+1,jj ) &
	1201	& + ftse(ji,jj ) * zhV(ji ,jj-1) &
	1202	& + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
	1203	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
	1204	& + ftse(ji,jj+1) * zhU(ji ,jj+1) &
	1205	& + ftnw(ji,jj ) * zhU(ji-1,jj ) &
	1206	& + ftne(ji,jj ) * zhU(ji ,jj ) )
	1207	END_2D
[11536]	1208	!
	1209	END SELECT
	1210	!
	1211	END SUBROUTINE dyn_cor_2D
	1212
	1213
	1214	SUBROUTINE wad_tmsk( pssh, ptmsk )
	1215	!!----------------------------------------------------------------------
	1216	!! * ROUTINE wad_lmt *
	1217	!!
	1218	!! ** Purpose : set wetting & drying mask at tracer points
	1219	!! for the current barotropic sub-step
	1220	!!
	1221	!! ** Method : ???
	1222	!!
	1223	!! ** Action : ptmsk : wetting & drying t-mask
	1224	!!----------------------------------------------------------------------
	1225	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh !
	1226	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk !
	1227	!
	1228	INTEGER :: ji, jj ! dummy loop indices
	1229	!!----------------------------------------------------------------------
	1230	!
	1231	IF( ln_wd_dl_rmp ) THEN
[13295]	1232	DO_2D( 1, 1, 1, 1 )
[12377]	1233	IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
	1234	! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
	1235	ptmsk(ji,jj) = 1._wp
	1236	ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
	1237	ptmsk(ji,jj) = TANH( 50._wp( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1) )
	1238	ELSE
	1239	ptmsk(ji,jj) = 0._wp
	1240	ENDIF
	1241	END_2D
[11536]	1242	ELSE
[13295]	1243	DO_2D( 1, 1, 1, 1 )
[12377]	1244	IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp
	1245	ELSE ; ptmsk(ji,jj) = 0._wp
	1246	ENDIF
	1247	END_2D
[11536]	1248	ENDIF
	1249	!
	1250	END SUBROUTINE wad_tmsk
	1251
	1252
	1253	SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
	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 ) :: pTmsk ! W & D t-mask
	1265	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
	1266	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
	1267	!
	1268	INTEGER :: ji, jj ! dummy loop indices
	1269	!!----------------------------------------------------------------------
	1270	!
[13998]	1271	DO_2D( 1, 1, 1, 0 ) ! not jpi-column
[12377]	1272	IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj)
	1273	ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj)
	1274	ENDIF
	1275	phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
	1276	pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
	1277	END_2D
[11536]	1278	!
[13998]	1279	DO_2D( 1, 0, 1, 1 ) ! not jpj-row
[12377]	1280	IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj )
	1281	ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1)
	1282	ENDIF
	1283	phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
	1284	pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
	1285	END_2D
[11536]	1286	!
	1287	END SUBROUTINE wad_Umsk
	1288
	1289
[12377]	1290	SUBROUTINE wad_spg( pshn, zcpx, zcpy )
[11536]	1291	!!---------------------------------------------------------------------
	1292	!! * ROUTINE wad_sp *
	1293	!!
	1294	!! ** Purpose :
	1295	!!----------------------------------------------------------------------
	1296	INTEGER :: ji ,jj ! dummy loop indices
	1297	LOGICAL :: ll_tmp1, ll_tmp2
[12377]	1298	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn
[11536]	1299	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
	1300	!!----------------------------------------------------------------------
[13295]	1301	DO_2D( 0, 0, 0, 0 )
[12377]	1302	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > &
	1303	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
	1304	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji+1,jj) + ht_0(ji+1,jj) ) &
	1305	& > rn_wdmin1 + rn_wdmin2
	1306	ll_tmp2 = ( ABS( pshn(ji+1,jj) - pshn(ji ,jj)) > 1.E-12 ).AND.( &
	1307	& MAX( pshn(ji,jj) , pshn(ji+1,jj) ) > &
	1308	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
	1309	IF(ll_tmp1) THEN
	1310	zcpx(ji,jj) = 1.0_wp
	1311	ELSEIF(ll_tmp2) THEN
	1312	! no worries about pshn(ji+1,jj) - pshn(ji ,jj) = 0, it won't happen ! here
	1313	zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
	1314	& / (pshn(ji+1,jj) - pshn(ji ,jj)) )
	1315	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
	1316	ELSE
	1317	zcpx(ji,jj) = 0._wp
	1318	ENDIF
	1319	!
	1320	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji,jj+1) ) > &
	1321	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
	1322	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji,jj+1) + ht_0(ji,jj+1) ) &
	1323	& > rn_wdmin1 + rn_wdmin2
	1324	ll_tmp2 = ( ABS( pshn(ji,jj) - pshn(ji,jj+1)) > 1.E-12 ).AND.( &
	1325	& MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > &
	1326	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
	1327
	1328	IF(ll_tmp1) THEN
	1329	zcpy(ji,jj) = 1.0_wp
	1330	ELSE IF(ll_tmp2) THEN
	1331	! no worries about pshn(ji,jj+1) - pshn(ji,jj ) = 0, it won't happen ! here
	1332	zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
	1333	& / (pshn(ji,jj+1) - pshn(ji,jj )) )
	1334	zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) )
	1335	ELSE
	1336	zcpy(ji,jj) = 0._wp
	1337	ENDIF
	1338	END_2D
[11536]	1339
	1340	END SUBROUTINE wad_spg
	1341
	1342
[12377]	1343	SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
[11536]	1344	!!----------------------------------------------------------------------
	1345	!! * ROUTINE dyn_drg_init *
	1346	!!
	1347	!! ** Purpose : - add the baroclinic top/bottom drag contribution to
	1348	!! the baroclinic part of the barotropic RHS
	1349	!! - compute the barotropic drag coefficients
	1350	!!
	1351	!! ** Method : computation done over the INNER domain only
	1352	!!----------------------------------------------------------------------
[12377]	1353	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
	1354	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation
	1355	REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels
	1356	REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS
	1357	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients
[11536]	1358	!
	1359	INTEGER :: ji, jj ! dummy loop indices
	1360	INTEGER :: ikbu, ikbv, iktu, iktv
	1361	REAL(wp) :: zztmp
	1362	REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i
	1363	!!----------------------------------------------------------------------
	1364	!
	1365	! !== Set the barotropic drag coef. ==!
	1366	!
[13998]	1367	IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities)
[11536]	1368
[13295]	1369	DO_2D( 0, 0, 0, 0 )
[12377]	1370	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
	1371	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
	1372	END_2D
[11536]	1373	ELSE ! bottom friction only
[13295]	1374	DO_2D( 0, 0, 0, 0 )
[12377]	1375	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
	1376	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
	1377	END_2D
[11536]	1378	ENDIF
	1379	!
	1380	! !== BOTTOM stress contribution from baroclinic velocities ==!
	1381	!
	1382	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
	1383
[13295]	1384	DO_2D( 0, 0, 0, 0 )
[12377]	1385	ikbu = mbku(ji,jj)
	1386	ikbv = mbkv(ji,jj)
	1387	zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
	1388	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
	1389	END_2D
[11536]	1390	ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
	1391
[13295]	1392	DO_2D( 0, 0, 0, 0 )
[12377]	1393	ikbu = mbku(ji,jj)
	1394	ikbv = mbkv(ji,jj)
	1395	zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
	1396	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
	1397	END_2D
[11536]	1398	ENDIF
	1399	!
	1400	IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please !
[12489]	1401	zztmp = -1._wp / rDt_e
[13295]	1402	DO_2D( 0, 0, 0, 0 )
[12377]	1403	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
	1404	& r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp )
	1405	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &
	1406	& r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
	1407	END_2D
[11536]	1408	ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
	1409
[13295]	1410	DO_2D( 0, 0, 0, 0 )
[12377]	1411	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)
	1412	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)
	1413	END_2D
[11536]	1414	END IF
	1415	!
	1416	! !== TOP stress contribution from baroclinic velocities ==! (no W/D case)
	1417	!
[13998]	1418	IF( ln_isfcav.OR.ln_drgice_imp ) THEN
[11536]	1419	!
	1420	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
	1421
[13295]	1422	DO_2D( 0, 0, 0, 0 )
[12377]	1423	iktu = miku(ji,jj)
	1424	iktv = mikv(ji,jj)
	1425	zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
	1426	zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
	1427	END_2D
[11536]	1428	ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
	1429
[13295]	1430	DO_2D( 0, 0, 0, 0 )
[12377]	1431	iktu = miku(ji,jj)
	1432	iktv = mikv(ji,jj)
	1433	zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
	1434	zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
	1435	END_2D
[11536]	1436	ENDIF
	1437	!
	1438	! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
	1439
[13295]	1440	DO_2D( 0, 0, 0, 0 )
[12377]	1441	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)
	1442	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)
	1443	END_2D
[11536]	1444	!
	1445	ENDIF
	1446	!
	1447	END SUBROUTINE dyn_drg_init
	1448
	1449	SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in
	1450	& za0, za1, za2, za3 ) ! ==> out
	1451	!!----------------------------------------------------------------------
	1452	INTEGER ,INTENT(in ) :: jn ! index of sub time step
	1453	LOGICAL ,INTENT(in ) :: ll_init !
	1454	REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient
	1455	!
	1456	REAL(wp) :: zepsilon, zgamma ! - -
	1457	!!----------------------------------------------------------------------
	1458	! ! set Half-step back interpolation coefficient
	1459	IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward
	1460	za0 = 1._wp
	1461	za1 = 0._wp
	1462	za2 = 0._wp
	1463	za3 = 0._wp
	1464	ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
	1465	za0 = 1.0833333333333_wp ! za0 = 1-gam-eps
	1466	za1 =-0.1666666666666_wp ! za1 = gam
	1467	za2 = 0.0833333333333_wp ! za2 = eps
	1468	za3 = 0._wp
	1469	ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
	1470	IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion
	1471	za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps
	1472	za1 = 0.285_wp ! za1 = 1/2 - 2gam - 3eps
	1473	za2 = 0.088_wp ! za2 = gam
	1474	za3 = 0.013_wp ! za3 = eps
	1475	ELSE ! no time diffusion
	1476	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
	1477	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
	1478	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
	1479	za1 = 1._wp - za0 - zgamma - zepsilon
	1480	za2 = zgamma
	1481	za3 = zepsilon
	1482	ENDIF
	1483	ENDIF
	1484	END SUBROUTINE ts_bck_interp
	1485
	1486
[358]	1487	!!======================================================================
	1488	END MODULE dynspg_ts

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

Download in other formats: