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 @ 13427

Last change on this file since 13427 was 13427, checked in by techene, 4 years ago
debug in order to pass non linear SETTE test when agrif in not activated see #2385
Property svn:keywords set to `Id`
File size: 80.9 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
[13295]	266	DO_2D( 0, 0, 0, 0 )
[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	!
[13295]	281	DO_2D( 0, 0, 0, 0 )
[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
[13295]	479	DO_2D( 1, 1, 1, 0 )
[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
	487	DO_2D( 1, 1, 1, 0 )
[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
[13295]	492	DO_2D( 1, 0, 1, 1 )
[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	! ! ---------------
[12489]	934	IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN !* Read the restart file
[13286]	935	CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	936	CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
	937	CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	938	CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
[9506]	939	IF( .NOT.ln_bt_av ) THEN
[13286]	940	CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp, ldxios = lrxios )
	941	CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	942	CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
	943	CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp, ldxios = lrxios )
	944	CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	945	CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
[9506]	946	ENDIF
[4486]	947	#if defined key_agrif
[9506]	948	! Read time integrated fluxes
	949	IF ( .NOT.Agrif_Root() ) THEN
[13286]	950	CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios )
	951	CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios )
[9506]	952	ENDIF
	953	#endif
	954	ELSE !* Start from rest
	955	IF(lwp) WRITE(numout,*)
	956	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
	957	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
	958	un_adv(:,:) = 0._wp ; vn_adv(:,:) = 0._wp ! used in the 1st interpol of agrif
	959	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
	960	#if defined key_agrif
	961	IF ( .NOT.Agrif_Root() ) THEN
	962	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
	963	ENDIF
	964	#endif
[4486]	965	ENDIF
[9506]	966	!
	967	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
	968	! ! -------------------
	969	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
[9367]	970	IF( lwxios ) CALL iom_swap( cwxios_context )
	971	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios )
	972	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios )
	973	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios )
	974	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios )
[4292]	975	!
	976	IF (.NOT.ln_bt_av) THEN
[9367]	977	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:), ldxios = lwxios )
	978	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:), ldxios = lwxios )
	979	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:), ldxios = lwxios )
	980	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:), ldxios = lwxios )
	981	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:), ldxios = lwxios )
	982	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:), ldxios = lwxios )
[4292]	983	ENDIF
[4486]	984	#if defined key_agrif
	985	! Save time integrated fluxes
	986	IF ( .NOT.Agrif_Root() ) THEN
[9367]	987	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lwxios )
	988	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lwxios )
[4486]	989	ENDIF
	990	#endif
[9367]	991	IF( lwxios ) CALL iom_swap( cxios_context )
[4292]	992	ENDIF
	993	!
	994	END SUBROUTINE ts_rst
[2528]	995
[6140]	996
	997	SUBROUTINE dyn_spg_ts_init
[4292]	998	!!---------------------------------------------------------------------
	999	!! * ROUTINE dyn_spg_ts_init *
	1000	!!
	1001	!! ** Purpose : Set time splitting options
	1002	!!----------------------------------------------------------------------
[6140]	1003	INTEGER :: ji ,jj ! dummy loop indices
	1004	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
[9019]	1005	REAL(wp), DIMENSION(jpi,jpj) :: zcu
[4292]	1006	!!----------------------------------------------------------------------
[4370]	1007	!
[5930]	1008	! Max courant number for ext. grav. waves
[4370]	1009	!
[13295]	1010	DO_2D( 0, 0, 0, 0 )
[12377]	1011	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
	1012	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
	1013	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
	1014	END_2D
[5930]	1015	!
[13286]	1016	zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) )
[10425]	1017	CALL mpp_max( 'dynspg_ts', zcmax )
[2528]	1018
[4370]	1019	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
[12489]	1020	IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
[4292]	1021
[12489]	1022	rDt_e = rn_Dt / REAL( nn_e , wp )
	1023	zcmax = zcmax * rDt_e
[9023]	1024	! Print results
[4292]	1025	IF(lwp) WRITE(numout,*)
[9169]	1026	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
	1027	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
[5930]	1028	IF( ln_bt_auto ) THEN
[12489]	1029	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e '
[4370]	1030	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
[4292]	1031	ELSE
[12489]	1032	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e
[358]	1033	ENDIF
[4292]	1034
	1035	IF(ln_bt_av) THEN
[12489]	1036	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on '
[4292]	1037	ELSE
[9169]	1038	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
[4292]	1039	ENDIF
[508]	1040	!
[4292]	1041	!
	1042	IF(ln_bt_fw) THEN
[4370]	1043	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
[4292]	1044	ELSE
[4370]	1045	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
[4292]	1046	ENDIF
	1047	!
[4486]	1048	#if defined key_agrif
	1049	! Restrict the use of Agrif to the forward case only
[9023]	1050	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
[4486]	1051	#endif
	1052	!
[4370]	1053	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
[4292]	1054	SELECT CASE ( nn_bt_flt )
[6140]	1055	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
[12489]	1056	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e'
	1057	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_e'
[9169]	1058	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
[4292]	1059	END SELECT
	1060	!
[4370]	1061	IF(lwp) WRITE(numout,*) ' '
[12489]	1062	IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e
	1063	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rDt_e
[4370]	1064	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
	1065	!
[9023]	1066	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
	1067	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
	1068	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
	1069	ENDIF
	1070	!
[6140]	1071	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
[4292]	1072	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
	1073	ENDIF
[6140]	1074	IF( zcmax>0.9_wp ) THEN
[12489]	1075	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )
[4292]	1076	ENDIF
	1077	!
[9124]	1078	! ! Allocate time-splitting arrays
	1079	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
	1080	!
	1081	! ! read restart when needed
[9506]	1082	CALL ts_rst( nit000, 'READ' )
[9124]	1083	!
[9367]	1084	IF( lwxios ) THEN
	1085	! define variables in restart file when writing with XIOS
	1086	CALL iom_set_rstw_var_active('ub2_b')
	1087	CALL iom_set_rstw_var_active('vb2_b')
	1088	CALL iom_set_rstw_var_active('un_bf')
	1089	CALL iom_set_rstw_var_active('vn_bf')
	1090	!
	1091	IF (.NOT.ln_bt_av) THEN
	1092	CALL iom_set_rstw_var_active('sshbb_e')
	1093	CALL iom_set_rstw_var_active('ubb_e')
	1094	CALL iom_set_rstw_var_active('vbb_e')
	1095	CALL iom_set_rstw_var_active('sshb_e')
	1096	CALL iom_set_rstw_var_active('ub_e')
	1097	CALL iom_set_rstw_var_active('vb_e')
	1098	ENDIF
	1099	#if defined key_agrif
	1100	! Save time integrated fluxes
	1101	IF ( .NOT.Agrif_Root() ) THEN
	1102	CALL iom_set_rstw_var_active('ub2_i_b')
	1103	CALL iom_set_rstw_var_active('vb2_i_b')
	1104	ENDIF
	1105	#endif
	1106	ENDIF
	1107	!
[4292]	1108	END SUBROUTINE dyn_spg_ts_init
[508]	1109
[11536]	1110
[12377]	1111	SUBROUTINE dyn_cor_2D_init( Kmm )
[11536]	1112	!!---------------------------------------------------------------------
[12377]	1113	!! * ROUTINE dyn_cor_2D_init *
[11536]	1114	!!
	1115	!! ** Purpose : Set time splitting options
	1116	!! Set arrays to remove/compute coriolis trend.
	1117	!! Do it once during initialization if volume is fixed, else at each long time step.
	1118	!! Note that these arrays are also used during barotropic loop. These are however frozen
	1119	!! although they should be updated in the variable volume case. Not a big approximation.
	1120	!! To remove this approximation, copy lines below inside barotropic loop
	1121	!! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
	1122	!!
	1123	!! Compute zwz = f / ( height of the water colomn )
	1124	!!----------------------------------------------------------------------
[12377]	1125	INTEGER, INTENT(in) :: Kmm ! Time index
[11536]	1126	INTEGER :: ji ,jj, jk ! dummy loop indices
	1127	REAL(wp) :: z1_ht
	1128	REAL(wp), DIMENSION(jpi,jpj) :: zhf
	1129	!!----------------------------------------------------------------------
	1130	!
	1131	SELECT CASE( nvor_scheme )
[13237]	1132	CASE( np_EEN ) != EEN scheme using e3f energy & enstrophy scheme
[11536]	1133	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
	1134	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
[13295]	1135	DO_2D( 1, 0, 1, 0 )
[13427]	1136	zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) &
	1137	& + ht(ji,jj ) + ht(ji+1,jj ) ) * 0.25_wp
[12377]	1138	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
	1139	END_2D
[11536]	1140	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
[13295]	1141	DO_2D( 1, 0, 1, 0 )
[13427]	1142	zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) &
	1143	& + ht(ji,jj ) + ht(ji+1,jj ) ) &
	1144	& / ( MAX(ssmask(ji,jj+1) + ssmask(ji+1,jj+1) &
	1145	& + ssmask(ji,jj ) + ssmask(ji+1,jj ) , 1._wp ) )
[12377]	1146	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
	1147	END_2D
[11536]	1148	END SELECT
	1149	CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
	1150	!
[13427]	1151	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[13295]	1152	DO_2D( 0, 1, 0, 1 )
[12377]	1153	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	1154	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	1155	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	1156	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	1157	END_2D
[11536]	1158	!
[13237]	1159	CASE( np_EET ) != EEN scheme using e3t energy conserving scheme
[13427]	1160	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[13295]	1161	DO_2D( 0, 1, 0, 1 )
[12377]	1162	z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
	1163	ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht
	1164	ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht
	1165	ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
	1166	ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht
	1167	END_2D
[11536]	1168	!
	1169	CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !
	1170	!
	1171	zwz(:,:) = 0._wp
	1172	zhf(:,:) = 0._wp
	1173
	1174	!!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
	1175	!!gm A priori a better value should be something like :
	1176	!!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
	1177	!!gm divided by the sum of the corresponding mask
	1178	!!gm
	1179	!!
	1180	IF( .NOT.ln_sco ) THEN
	1181
	1182	!!gm agree the JC comment : this should be done in a much clear way
	1183
	1184	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
	1185	! Set it to zero for the time being
	1186	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
	1187	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
	1188	! ENDIF
	1189	! zhf(:,:) = gdepw_0(:,:,jk+1)
	1190	!
	1191	ELSE
	1192	!
	1193	!zhf(:,:) = hbatf(:,:)
[13295]	1194	DO_2D( 1, 0, 1, 0 )
[12377]	1195	zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &
	1196	& + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &
	1197	& / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &
	1198	& + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )
	1199	END_2D
[11536]	1200	ENDIF
	1201	!
[13427]	1202	DO jj = 1, jpjm1 ! keep only the value at the coast if sco
[11536]	1203	zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
	1204	END DO
	1205	!
[13427]	1206	DO jk = 1, jpkm1 ! ocean point : sum of masked e3f
[11536]	1207	DO jj = 1, jpjm1
[12377]	1208	zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
[11536]	1209	END DO
	1210	END DO
	1211	CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
	1212	! JC: TBC. hf should be greater than 0
[13295]	1213	DO_2D( 1, 1, 1, 1 )
[12377]	1214	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj)
	1215	END_2D
[11536]	1216	zwz(:,:) = ff_f(:,:) * zwz(:,:)
	1217	END SELECT
	1218
	1219	END SUBROUTINE dyn_cor_2d_init
	1220
	1221
	1222
[13237]	1223	SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )
[11536]	1224	!!---------------------------------------------------------------------
	1225	!! * ROUTINE dyn_cor_2d *
	1226	!!
	1227	!! ** Purpose : Compute u and v coriolis trends
	1228	!!----------------------------------------------------------------------
	1229	INTEGER :: ji ,jj ! dummy loop indices
	1230	REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - -
[13237]	1231	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhU, zhV
[11536]	1232	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd
	1233	!!----------------------------------------------------------------------
	1234	SELECT CASE( nvor_scheme )
	1235	CASE( np_ENT ) ! enstrophy conserving scheme (f-point)
[13295]	1236	DO_2D( 0, 0, 0, 0 )
[12377]	1237	z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
	1238	z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
	1239	zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu &
[13237]	1240	& * ( e1e2t(ji+1,jj)pht(ji+1,jj)ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) &
	1241	& + e1e2t(ji ,jj)pht(ji ,jj)ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) )
[12377]	1242	!
	1243	zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv &
[13237]	1244	& * ( e1e2t(ji,jj+1)pht(ji,jj+1)ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) &
	1245	& + e1e2t(ji,jj )pht(ji,jj )ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) )
[12377]	1246	END_2D
[11536]	1247	!
	1248	CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX
[13295]	1249	DO_2D( 0, 0, 0, 0 )
[12377]	1250	zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
	1251	zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
	1252	zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
	1253	zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
	1254	! energy conserving formulation for planetary vorticity term
	1255	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	1256	zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
	1257	END_2D
[11536]	1258	!
	1259	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
[13295]	1260	DO_2D( 0, 0, 0, 0 )
[12377]	1261	zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) &
	1262	& + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)
	1263	zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) &
	1264	& + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)
	1265	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	1266	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
	1267	END_2D
[11536]	1268	!
	1269	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
[13295]	1270	DO_2D( 0, 0, 0, 0 )
[12377]	1271	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) &
	1272	& + ftnw(ji+1,jj) * zhV(ji+1,jj ) &
	1273	& + ftse(ji,jj ) * zhV(ji ,jj-1) &
	1274	& + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
	1275	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
	1276	& + ftse(ji,jj+1) * zhU(ji ,jj+1) &
	1277	& + ftnw(ji,jj ) * zhU(ji-1,jj ) &
	1278	& + ftne(ji,jj ) * zhU(ji ,jj ) )
	1279	END_2D
[11536]	1280	!
	1281	END SELECT
	1282	!
	1283	END SUBROUTINE dyn_cor_2D
	1284
	1285
	1286	SUBROUTINE wad_tmsk( pssh, ptmsk )
	1287	!!----------------------------------------------------------------------
	1288	!! * ROUTINE wad_lmt *
	1289	!!
	1290	!! ** Purpose : set wetting & drying mask at tracer points
	1291	!! for the current barotropic sub-step
	1292	!!
	1293	!! ** Method : ???
	1294	!!
	1295	!! ** Action : ptmsk : wetting & drying t-mask
	1296	!!----------------------------------------------------------------------
	1297	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh !
	1298	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk !
	1299	!
	1300	INTEGER :: ji, jj ! dummy loop indices
	1301	!!----------------------------------------------------------------------
	1302	!
	1303	IF( ln_wd_dl_rmp ) THEN
[13295]	1304	DO_2D( 1, 1, 1, 1 )
[12377]	1305	IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
	1306	! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
	1307	ptmsk(ji,jj) = 1._wp
	1308	ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
	1309	ptmsk(ji,jj) = TANH( 50._wp( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1) )
	1310	ELSE
	1311	ptmsk(ji,jj) = 0._wp
	1312	ENDIF
	1313	END_2D
[11536]	1314	ELSE
[13295]	1315	DO_2D( 1, 1, 1, 1 )
[12377]	1316	IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp
	1317	ELSE ; ptmsk(ji,jj) = 0._wp
	1318	ENDIF
	1319	END_2D
[11536]	1320	ENDIF
	1321	!
	1322	END SUBROUTINE wad_tmsk
	1323
	1324
	1325	SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
	1326	!!----------------------------------------------------------------------
	1327	!! * ROUTINE wad_lmt *
	1328	!!
	1329	!! ** Purpose : set wetting & drying mask at tracer points
	1330	!! for the current barotropic sub-step
	1331	!!
	1332	!! ** Method : ???
	1333	!!
	1334	!! ** Action : ptmsk : wetting & drying t-mask
	1335	!!----------------------------------------------------------------------
	1336	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
	1337	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
	1338	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
	1339	!
	1340	INTEGER :: ji, jj ! dummy loop indices
	1341	!!----------------------------------------------------------------------
	1342	!
[13295]	1343	DO_2D( 1, 1, 1, 0 )
[12377]	1344	IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj)
	1345	ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj)
	1346	ENDIF
	1347	phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
	1348	pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
	1349	END_2D
[11536]	1350	!
[13295]	1351	DO_2D( 1, 0, 1, 1 )
[12377]	1352	IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj )
	1353	ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1)
	1354	ENDIF
	1355	phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
	1356	pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
	1357	END_2D
[11536]	1358	!
	1359	END SUBROUTINE wad_Umsk
	1360
	1361
[12377]	1362	SUBROUTINE wad_spg( pshn, zcpx, zcpy )
[11536]	1363	!!---------------------------------------------------------------------
	1364	!! * ROUTINE wad_sp *
	1365	!!
	1366	!! ** Purpose :
	1367	!!----------------------------------------------------------------------
	1368	INTEGER :: ji ,jj ! dummy loop indices
	1369	LOGICAL :: ll_tmp1, ll_tmp2
[12377]	1370	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn
[11536]	1371	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
	1372	!!----------------------------------------------------------------------
[13295]	1373	DO_2D( 0, 0, 0, 0 )
[12377]	1374	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > &
	1375	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
	1376	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji+1,jj) + ht_0(ji+1,jj) ) &
	1377	& > rn_wdmin1 + rn_wdmin2
	1378	ll_tmp2 = ( ABS( pshn(ji+1,jj) - pshn(ji ,jj)) > 1.E-12 ).AND.( &
	1379	& MAX( pshn(ji,jj) , pshn(ji+1,jj) ) > &
	1380	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
	1381	IF(ll_tmp1) THEN
	1382	zcpx(ji,jj) = 1.0_wp
	1383	ELSEIF(ll_tmp2) THEN
	1384	! no worries about pshn(ji+1,jj) - pshn(ji ,jj) = 0, it won't happen ! here
	1385	zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
	1386	& / (pshn(ji+1,jj) - pshn(ji ,jj)) )
	1387	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
	1388	ELSE
	1389	zcpx(ji,jj) = 0._wp
	1390	ENDIF
	1391	!
	1392	ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji,jj+1) ) > &
	1393	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
	1394	& MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji,jj+1) + ht_0(ji,jj+1) ) &
	1395	& > rn_wdmin1 + rn_wdmin2
	1396	ll_tmp2 = ( ABS( pshn(ji,jj) - pshn(ji,jj+1)) > 1.E-12 ).AND.( &
	1397	& MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > &
	1398	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
	1399
	1400	IF(ll_tmp1) THEN
	1401	zcpy(ji,jj) = 1.0_wp
	1402	ELSE IF(ll_tmp2) THEN
	1403	! no worries about pshn(ji,jj+1) - pshn(ji,jj ) = 0, it won't happen ! here
	1404	zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
	1405	& / (pshn(ji,jj+1) - pshn(ji,jj )) )
	1406	zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) )
	1407	ELSE
	1408	zcpy(ji,jj) = 0._wp
	1409	ENDIF
	1410	END_2D
[11536]	1411
	1412	END SUBROUTINE wad_spg
	1413
	1414
	1415
[12377]	1416	SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
[11536]	1417	!!----------------------------------------------------------------------
	1418	!! * ROUTINE dyn_drg_init *
	1419	!!
	1420	!! ** Purpose : - add the baroclinic top/bottom drag contribution to
	1421	!! the baroclinic part of the barotropic RHS
	1422	!! - compute the barotropic drag coefficients
	1423	!!
	1424	!! ** Method : computation done over the INNER domain only
	1425	!!----------------------------------------------------------------------
[12377]	1426	INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices
	1427	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation
	1428	REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels
	1429	REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS
	1430	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients
[11536]	1431	!
	1432	INTEGER :: ji, jj ! dummy loop indices
	1433	INTEGER :: ikbu, ikbv, iktu, iktv
	1434	REAL(wp) :: zztmp
	1435	REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i
	1436	!!----------------------------------------------------------------------
	1437	!
	1438	! !== Set the barotropic drag coef. ==!
	1439	!
	1440	IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)
	1441
[13295]	1442	DO_2D( 0, 0, 0, 0 )
[12377]	1443	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
	1444	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
	1445	END_2D
[11536]	1446	ELSE ! bottom friction only
[13295]	1447	DO_2D( 0, 0, 0, 0 )
[12377]	1448	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
	1449	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
	1450	END_2D
[11536]	1451	ENDIF
	1452	!
	1453	! !== BOTTOM stress contribution from baroclinic velocities ==!
	1454	!
	1455	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
	1456
[13295]	1457	DO_2D( 0, 0, 0, 0 )
[12377]	1458	ikbu = mbku(ji,jj)
	1459	ikbv = mbkv(ji,jj)
	1460	zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
	1461	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
	1462	END_2D
[11536]	1463	ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
	1464
[13295]	1465	DO_2D( 0, 0, 0, 0 )
[12377]	1466	ikbu = mbku(ji,jj)
	1467	ikbv = mbkv(ji,jj)
	1468	zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
	1469	zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
	1470	END_2D
[11536]	1471	ENDIF
	1472	!
	1473	IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please !
[12489]	1474	zztmp = -1._wp / rDt_e
[13295]	1475	DO_2D( 0, 0, 0, 0 )
[12377]	1476	pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
	1477	& r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp )
	1478	pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &
	1479	& r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
	1480	END_2D
[11536]	1481	ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
	1482
[13295]	1483	DO_2D( 0, 0, 0, 0 )
[12377]	1484	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)
	1485	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)
	1486	END_2D
[11536]	1487	END IF
	1488	!
	1489	! !== TOP stress contribution from baroclinic velocities ==! (no W/D case)
	1490	!
	1491	IF( ln_isfcav ) THEN
	1492	!
	1493	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
	1494
[13295]	1495	DO_2D( 0, 0, 0, 0 )
[12377]	1496	iktu = miku(ji,jj)
	1497	iktv = mikv(ji,jj)
	1498	zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
	1499	zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
	1500	END_2D
[11536]	1501	ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
	1502
[13295]	1503	DO_2D( 0, 0, 0, 0 )
[12377]	1504	iktu = miku(ji,jj)
	1505	iktv = mikv(ji,jj)
	1506	zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
	1507	zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
	1508	END_2D
[11536]	1509	ENDIF
	1510	!
	1511	! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
	1512
[13295]	1513	DO_2D( 0, 0, 0, 0 )
[12377]	1514	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)
	1515	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)
	1516	END_2D
[11536]	1517	!
	1518	ENDIF
	1519	!
	1520	END SUBROUTINE dyn_drg_init
	1521
	1522	SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in
	1523	& za0, za1, za2, za3 ) ! ==> out
	1524	!!----------------------------------------------------------------------
	1525	INTEGER ,INTENT(in ) :: jn ! index of sub time step
	1526	LOGICAL ,INTENT(in ) :: ll_init !
	1527	REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient
	1528	!
	1529	REAL(wp) :: zepsilon, zgamma ! - -
	1530	!!----------------------------------------------------------------------
	1531	! ! set Half-step back interpolation coefficient
	1532	IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward
	1533	za0 = 1._wp
	1534	za1 = 0._wp
	1535	za2 = 0._wp
	1536	za3 = 0._wp
	1537	ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
	1538	za0 = 1.0833333333333_wp ! za0 = 1-gam-eps
	1539	za1 =-0.1666666666666_wp ! za1 = gam
	1540	za2 = 0.0833333333333_wp ! za2 = eps
	1541	za3 = 0._wp
	1542	ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
	1543	IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion
	1544	za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps
	1545	za1 = 0.285_wp ! za1 = 1/2 - 2gam - 3eps
	1546	za2 = 0.088_wp ! za2 = gam
	1547	za3 = 0.013_wp ! za3 = eps
	1548	ELSE ! no time diffusion
	1549	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
	1550	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
	1551	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
	1552	za1 = 1._wp - za0 - zgamma - zepsilon
	1553	za2 = zgamma
	1554	za3 = zepsilon
	1555	ENDIF
	1556	ENDIF
	1557	END SUBROUTINE ts_bck_interp
	1558
	1559
[358]	1560	!!======================================================================
	1561	END MODULE dynspg_ts

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

Download in other formats: