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

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

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

Download in other formats: