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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynspg_ts.F90 in NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynspg_ts.F90 @ 14618

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

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

Download in other formats: