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/trunk/src/OCE/DYN – NEMO

Context Navigation

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

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

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

Download in other formats: