Context Navigation

dynspg_ts.F90 @ 15814

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/trunk/src/OCE/DYN/dynspg_ts.F90 @ 15814

Download in other formats: