Context Navigation

dynspg_ts.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Property svn:keywords set to Id

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

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

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

Context Navigation

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

Download in other formats: