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/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DYN/dynspg_ts.F90 @ 11380

Last change on this file since 11380 was 11380, checked in by girrmann, 5 years ago
dev_r10984_HPC-13 : adding extra halos in dyn_spg_ts is now possible, only works with a single halo when used with tide or bdy, see #2308
Property svn:keywords set to `Id`
File size: 96.4 KB

Line
1	MODULE dynspg_ts
2
3	!! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out !
4
5	!!======================================================================
6	!! * MODULE dynspg_ts *
7	!! Ocean dynamics: surface pressure gradient trend, split-explicit scheme
8	!!======================================================================
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
15	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
16	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
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
19	!! 3.7 ! 2015-11 (J. Chanut) free surface simplification
20	!! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence
21	!! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90)
22	!!---------------------------------------------------------------------
23
24	!!----------------------------------------------------------------------
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
29	!!----------------------------------------------------------------------
30	USE oce ! ocean dynamics and tracers
31	USE dom_oce ! ocean space and time domain
32	USE sbc_oce ! surface boundary condition: ocean
33	USE zdf_oce ! vertical physics: variables
34	USE zdfdrg ! vertical physics: top/bottom drag coef.
35	USE sbcisf ! ice shelf variable (fwfisf)
36	USE sbcapr ! surface boundary condition: atmospheric pressure
37	USE dynadv , ONLY: ln_dynadv_vec
38	USE dynvor ! vortivity scheme indicators
39	USE phycst ! physical constants
40	USE dynvor ! vorticity term
41	USE wet_dry ! wetting/drying flux limter
42	USE bdy_oce ! open boundary
43	USE bdyvol ! open boundary volume conservation
44	USE bdytides ! open boundary condition data
45	USE bdydyn2d ! open boundary conditions on barotropic variables
46	USE sbctide ! tides
47	USE updtide ! tide potential
48	USE sbcwave ! surface wave
49	USE diatmb ! Top,middle,bottom output
50	#if defined key_agrif
51	USE agrif_oce_interp ! agrif
52	USE agrif_oce
53	#endif
54	#if defined key_asminc
55	USE asminc ! Assimilation increment
56	#endif
57	!
58	USE in_out_manager ! I/O manager
59	USE lib_mpp ! distributed memory computing library
60	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
61	USE prtctl ! Print control
62	USE iom ! IOM library
63	USE restart ! only for lrst_oce
64	USE diatmb ! Top,middle,bottom output
65
66
67	IMPLICIT NONE
68	PRIVATE
69
70	PUBLIC dyn_spg_ts ! called by dyn_spg
71	PUBLIC dyn_spg_ts_init ! - - dyn_spg_init
72
73	!! Time filtered arrays at baroclinic time step:
74	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step
75	!
76	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
77	REAL(wp),SAVE :: rdtbt ! Barotropic time step
78	!
79	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields
80	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points
81	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
82	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
83
84	!! Arrays at barotropic time step: ! befbefore! before ! now ! after !
85	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubb_e , ub_e , un_e , ua_e !: u-external velocity
86	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vbb_e , vb_e , vn_e , va_e !: v-external velocity
87	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, sshn_e, ssha_e !: external ssh
88	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_e !: external u-depth
89	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_e !: external v-depth
90	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hur_e !: inverse of u-depth
91	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hvr_e !: inverse of v-depth
92	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_b , vb2_b !: Half step fluxes (ln_bt_fw=T)
93	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_bf , vn_bf !: Asselin filtered half step fluxes (ln_bt_fw=T)
94
95	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0_xtd , hu_0_xtd , hv_0_xtd
96	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_n_xtd , hv_n_xtd
97	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_e1e2t_xtd, r1_e1e2u_xtd, r1_e1e2v_xtd
98	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2t_xtd
99	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: e2u_xtd , e1v_xtd
100	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_e1u_xtd, r1_e2v_xtd
101	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask_xtd, ssumask_xtd, ssvmask_xtd
102	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff_t_xtd ! used in ENT scheme
103
104	#if defined key_agrif
105	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_i_b, vb2_i_b !: Half step time integrated fluxes
106	#endif
107
108	REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios
109	REAL(wp) :: r1_8 = 0.125_wp !
110	REAL(wp) :: r1_4 = 0.25_wp !
111	REAL(wp) :: r1_2 = 0.5_wp !
112
113	!! * Substitutions
114	# include "vectopt_loop_substitute.h90"
115	!!----------------------------------------------------------------------
116	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
117	!! $Id$
118	!! Software governed by the CeCILL license (see ./LICENSE)
119	!!----------------------------------------------------------------------
120	CONTAINS
121
122	INTEGER FUNCTION dyn_spg_ts_alloc()
123	!!----------------------------------------------------------------------
124	!! * routine dyn_spg_ts_alloc *
125	!!----------------------------------------------------------------------
126	INTEGER :: ierr(6)
127	INTEGER :: idbi, idei, idbj, idej ! lower/upper bounds of extended arrays
128	!!----------------------------------------------------------------------
129	ierr(:) = 0
130	idbi = 1 - nn_hlts ; idbj = 1 - nn_hlts
131	idei = jpi + nn_hlts ; idej = jpj + nn_hlts
132	!
133	ALLOCATE( wgtbtp1(3nn_baro), wgtbtp2(3nn_baro), zwz(idbi:idei,idbj:idej), STAT=ierr(1) )
134	IF( ln_dynvor_een .OR. ln_dynvor_eeT ) THEN
135	ALLOCATE( ftnw(idbi:idei,idbj:idej) , ftne(idbi:idei,idbj:idej) , ftsw(idbi:idei,idbj:idej) , ftse(idbi:idei,idbj:idej) &
136	& , STAT=ierr(2) )
137	ELSEIF( ln_dynvor_enT ) THEN
138	ALLOCATE( ff_t_xtd(idbi:idei,idbj:idej), STAT=ierr(2) )
139	END IF
140	!
141	ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) )
142	ALLOCATE( ssha_e(idbi:idei,idbj:idej), sshn_e(idbi:idei,idbj:idej), sshb_e(idbi:idei,idbj:idej), sshbb_e(idbi:idei,idbj:idej) &
143	& , ua_e(idbi:idei,idbj:idej), un_e(idbi:idei,idbj:idej), ub_e(idbi:idei,idbj:idej), ubb_e(idbi:idei,idbj:idej) &
144	& , va_e(idbi:idei,idbj:idej), vn_e(idbi:idei,idbj:idej), vb_e(idbi:idei,idbj:idej), vbb_e(idbi:idei,idbj:idej) &
145	& , hu_e(idbi:idei,idbj:idej), hur_e(idbi:idei,idbj:idej), hv_e(idbi:idei,idbj:idej), hvr_e(idbi:idei,idbj:idej) &
146	& , STAT=ierr(4) )
147	!
148	ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_bf(jpi,jpj), vn_bf(jpi,jpj) , STAT=ierr(5) )
149
150	#if defined key_agrif
151	ALLOCATE( ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj) , STAT=ierr(5) )
152	#endif
153	ALLOCATE( ht_0_xtd (idbi:idei,idbj:idej), hu_0_xtd (idbi:idei,idbj:idej), hv_0_xtd (idbi:idei,idbj:idej) &
154	& , r1_e1e2t_xtd(idbi:idei,idbj:idej), r1_e1e2u_xtd(idbi:idei,idbj:idej), r1_e1e2v_xtd(idbi:idei,idbj:idej) &
155	& , ssmask_xtd (idbi:idei,idbj:idej), ssumask_xtd (idbi:idei,idbj:idej), ssvmask_xtd (idbi:idei,idbj:idej) &
156	& , e1e2t_xtd (idbi:idei,idbj:idej), e2u_xtd (idbi:idei,idbj:idej), e1v_xtd (idbi:idei,idbj:idej) &
157	& , r1_e1u_xtd (idbi:idei,idbj:idej), r1_e2v_xtd (idbi:idei,idbj:idej) &
158	& , STAT=ierr(6) )
159	!
160	dyn_spg_ts_alloc = MAXVAL( ierr(:) )
161	!
162	CALL mpp_sum( 'dynspg_ts', dyn_spg_ts_alloc )
163	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_stop( 'STOP', 'dyn_spg_ts_alloc: failed to allocate arrays' )
164	!
165	END FUNCTION dyn_spg_ts_alloc
166
167
168	SUBROUTINE dyn_spg_ts( kt )
169	!!----------------------------------------------------------------------
170	!!
171	!! ** Purpose : - Compute the now trend due to the explicit time stepping
172	!! of the quasi-linear barotropic system, and add it to the
173	!! general momentum trend.
174	!!
175	!! ** Method : - split-explicit schem (time splitting) :
176	!! Barotropic variables are advanced from internal time steps
177	!! "n" to "n+1" if ln_bt_fw=T
178	!! or from
179	!! "n-1" to "n+1" if ln_bt_fw=F
180	!! thanks to a generalized forward-backward time stepping (see ref. below).
181	!!
182	!! ** Action :
183	!! -Update the filtered free surface at step "n+1" : ssha
184	!! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
185	!! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv
186	!! These are used to advect tracers and are compliant with discrete
187	!! continuity equation taken at the baroclinic time steps. This
188	!! ensures tracers conservation.
189	!! - (ua, va) momentum trend updated with barotropic component.
190	!!
191	!! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
192	!!---------------------------------------------------------------------
193	INTEGER, INTENT(in) :: kt ! ocean time-step index
194	!
195	INTEGER :: ji, jj, jk, jm ! dummy loop indices
196	LOGICAL :: ll_fw_start ! =T : forward integration
197	LOGICAL :: ll_init ! =T : special startup of 2d equations
198	INTEGER :: noffset ! local integers : time offset for bdy update
199	REAL(wp) :: r1_2dt_b, z1_hu, z1_hv ! local scalars
200	REAL(wp) :: za0, za1, za2, za3 ! - -
201	REAL(wp) :: zmdi, zztmp, zldg ! - -
202	REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - -
203	REAL(wp) :: zun_save, zvn_save ! - -
204	!
205	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zu_trd, zssh_frc
206	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zv_trd
207	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zhU , zhV
208	!
209	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: hu_n_xtd, hv_n_xtd
210	!
211	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zu_spg , zv_spg
212	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zu_frc , zv_frc
213	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsshu_a, zhup2_e, zhtp2_e
214	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsshv_a, zhvp2_e, zsshp2_e
215	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zCdU_u , zCdU_v ! top/bottom stress at u- & v-points
216
217	!
218	REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True.
219
220	INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point
221
222	REAL(wp) :: zepsilon, zgamma ! - -
223	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef.
224	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
225	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask
226
227
228	INTEGER :: idbi, idei, idbj, idej ! lower/upper bounds of extended arrays
229	INTEGER :: ixtd ! number of halos over which the solution is currently correct
230	INTEGER :: ibi, iei, ibj, iej ! lower and upper bounds over which the solution is currently correct
231	!!----------------------------------------------------------------------
232	!
233
234
235	IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) )
236	! !* Allocate temporary arrays
237	IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj) )
238
239	idbi = 1 - nn_hlts ; idbj = 1 - nn_hlts
240	idei = jpi + nn_hlts ; idej = jpj + nn_hlts
241	! ! allocate local arrays
242	ALLOCATE( zu_spg (idbi:idei,idbj:idej), zv_spg (idbi:idei,idbj:idej) &
243	& , zsshu_a (idbi:idei,idbj:idej), zsshv_a (idbi:idei,idbj:idej) &
244	& , zhup2_e (idbi:idei,idbj:idej), zhvp2_e (idbi:idei,idbj:idej) &
245	& , zCdU_u (idbi:idei,idbj:idej), zCdU_v (idbi:idei,idbj:idej) &
246	& , zhtp2_e (idbi:idei,idbj:idej), zsshp2_e(idbi:idei,idbj:idej) &
247	& , zu_trd (idbi:idei,idbj:idej), zu_frc (idbi:idei,idbj:idej) &
248	& , zv_trd (idbi:idei,idbj:idej), zv_frc (idbi:idei,idbj:idej) &
249	& , zhU (idbi:idei,idbj:idej), zhV (idbi:idei,idbj:idej) &
250	& , zssh_frc(idbi:idei,idbj:idej) )
251	! ! allocate redundant arrays
252	ALLOCATE( hu_n_xtd(idbi:idei,idbj:idej), hv_n_xtd(idbi:idei,idbj:idej) )
253	!
254	zmdi=1.e+20 ! missing data indicator for masking
255	!
256	zwdramp = r_rn_wdmin1 ! simplest ramp
257	! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
258	! ! inverse of baroclinic time step
259	IF( kt == nit000 .AND. neuler == 0 ) THEN ; r1_2dt_b = 1._wp / ( rdt )
260	ELSE ; r1_2dt_b = 1._wp / ( 2._wp * rdt )
261	ENDIF
262	!
263	ll_init = ln_bt_av ! if no time averaging, then no specific restart
264	ll_fw_start = .FALSE.
265	! ! time offset in steps for bdy data update
266	IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro
267	ELSE ; noffset = 0
268	ENDIF
269	!
270	IF( kt == nit000 ) THEN !* initialisation
271	!
272	IF(lwp) WRITE(numout,*)
273	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
274	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
275	IF(lwp) WRITE(numout,*)
276	!
277	IF( neuler == 0 ) ll_init=.TRUE.
278	!
279	IF( ln_bt_fw .OR. neuler == 0 ) THEN
280	ll_fw_start =.TRUE.
281	noffset = 0
282	ELSE
283	ll_fw_start =.FALSE.
284	ENDIF
285	! ! Set averaging weights and cycle length:
286	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
287	!
288	ENDIF
289	!
290	! If forward start at previous time step, and centered integration,
291	! then update averaging weights:
292	IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
293	ll_fw_start=.FALSE.
294	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
295	ENDIF
296	!
297
298	! -----------------------------------------------------------------------------
299	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
300	! -----------------------------------------------------------------------------
301	!
302	!
303	! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends)
304	! ! --------------------------- !
305	zu_frc(1:jpi,1:jpj) = SUM( e3u_n(:,:,:) * ua(:,:,:) * umask(:,:,:) , DIM=3 ) * r1_hu_n(:,:)
306	zv_frc(1:jpi,1:jpj) = SUM( e3v_n(:,:,:) * va(:,:,:) * vmask(:,:,:) , DIM=3 ) * r1_hv_n(:,:)
307	!
308	!
309	! != Ua => baroclinic trend =! (remove its vertical mean)
310	DO jk = 1, jpkm1 ! ------------------------ !
311	ua(:,:,jk) = ( ua(:,:,jk) - zu_frc(1:jpi,1:jpj) ) * umask(:,:,jk)
312	va(:,:,jk) = ( va(:,:,jk) - zv_frc(1:jpi,1:jpj) ) * vmask(:,:,jk)
313	END DO
314
315	!!gm Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
316	!!gm Is it correct to do so ? I think so...
317
318	! != remove 2D Coriolis and pressure gradient trends =!
319	! ! ------------------------------------------------- !
320	!
321	! Set zwz, the barotropic Coriolis force coefficient
322	! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes
323	IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2d_init( idbi, idei, idbj, idej )
324	!
325	! !* 2D Coriolis trends
326	zhU(1:jpi,1:jpj) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes
327	zhV(1:jpi,1:jpj) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) ! NB: FULL domain : put a value in last row and column
328	!
329	! ! ht_n, hu_n, hv_n, un_b, vn_b are of size 1:jpi 1:jpj
330	! ! zhU, zhV, zu_trd, zv_trd are of size idbi:idei idbj:idej
331	CALL dyn_cor_2d( ht_n, hu_n, hv_n, un_b, vn_b, zhU, zhV, 1 , jpi , 1 , jpj & ! <<== in
332	& , zu_trd, zv_trd, idbi, idei, idbj, idej ) ! ==>> out
333	!
334	IF( .NOT.ln_linssh ) THEN !* surface pressure gradient (variable volume only)
335	!
336	IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg
337	CALL wad_spg( sshn, zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy
338	DO jj = 2, jpjm1
339	DO ji = 2, jpim1 ! SPG with the application of W/D gravity filters
340	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) &
341	& * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth
342	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) &
343	& * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth
344	END DO
345	END DO
346	ELSE ! now suface pressure gradient
347	DO jj = 2, jpjm1
348	DO ji = fs_2, fs_jpim1 ! vector opt.
349	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj)
350	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj)
351	END DO
352	END DO
353	ENDIF
354	!
355	ENDIF
356	!
357	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
358	DO ji = fs_2, fs_jpim1
359	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
360	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
361	END DO
362	END DO
363	!
364	! != Add bottom stress contribution from baroclinic velocities =!
365	! ! ----------------------------------------------------------- !
366	CALL drg_init( idbi, idei, idbj, idej, zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients
367	! ! arrays are computed on inner domain
368	!
369	! != Add atmospheric pressure forcing =!
370	! ! ---------------------------------- !
371	IF( ln_apr_dyn ) THEN
372	IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
373	DO jj = 2, jpjm1
374	DO ji = fs_2, fs_jpim1 ! vector opt.
375	zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
376	zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
377	END DO
378	END DO
379	ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
380	zztmp = grav * r1_2
381	DO jj = 2, jpjm1
382	DO ji = fs_2, fs_jpim1 ! vector opt.
383	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
384	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
385	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
386	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
387	END DO
388	END DO
389	ENDIF
390	ENDIF
391	!
392	! != Add atmospheric pressure forcing =!
393	! ! ---------------------------------- !
394	IF( ln_bt_fw ) THEN ! Add wind forcing
395	DO jj = 2, jpjm1
396	DO ji = fs_2, fs_jpim1 ! vector opt.
397	zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)
398	zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)
399	END DO
400	END DO
401	ELSE
402	zztmp = r1_rau0 * r1_2
403	DO jj = 2, jpjm1
404	DO ji = fs_2, fs_jpim1 ! vector opt.
405	zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)
406	zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
407	END DO
408	END DO
409	ENDIF
410	!
411	! !----------------!
412	! !== sssh_frc ==! Right-Hand-Side of the barotropic ssh equation (over the FULL domain)
413	! !----------------!
414	! != Net water flux forcing applied to a water column =!
415	! ! --------------------------------------------------- !
416	IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
417	zssh_frc(1:jpi,1:jpj) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
418	ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
419	zztmp = r1_rau0 * r1_2
420	zssh_frc(1:jpi,1:jpj) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) + fwfisf(:,:) + fwfisf_b(:,:) )
421	ENDIF
422	! != Add Stokes drift divergence =! (if exist)
423	IF( ln_sdw ) THEN ! ----------------------------- !
424	zssh_frc(1:jpi,1:jpj) = zssh_frc(1:jpi,1:jpj) + div_sd(:,:)
425	ENDIF
426	!
427	#if defined key_asminc
428	! != Add the IAU weighted SSH increment =!
429	! ! ------------------------------------ !
430	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
431	zssh_frc(1:jpi,1:jpj) = zssh_frc(1:jpi,1:jpj) - ssh_iau(:,:)
432	ENDIF
433	#endif
434	! != Fill boundary data arrays for AGRIF
435	! ! ------------------------------------
436	#if defined key_agrif
437	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
438	#endif
439
440
441	!
442	! -----------------------------------------------------------------------
443	! Phase 2 : Integration of the barotropic equations
444	! -----------------------------------------------------------------------
445	!
446	! ! ==================== !
447	! ! Initialisations !
448	! ! ==================== !
449	! Initialize barotropic variables:
450	IF( ll_init ) THEN
451	sshbb_e(:,:) = 0._wp
452	ubb_e (:,:) = 0._wp
453	vbb_e (:,:) = 0._wp
454	sshb_e (:,:) = 0._wp
455	ub_e (:,:) = 0._wp
456	vb_e (:,:) = 0._wp
457	ENDIF
458	!
459	IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0)
460	zhtp2_e(1:jpi,1:jpj) = ht_n(1:jpi,1:jpj)
461	zhup2_e(1:jpi,1:jpj) = hu_n(1:jpi,1:jpj)
462	zhvp2_e(1:jpi,1:jpj) = hv_n(1:jpi,1:jpj)
463	ENDIF
464	!
465	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
466	sshn_e(1:jpi,1:jpj) = sshn(1:jpi,1:jpj)
467	un_e (1:jpi,1:jpj) = un_b(1:jpi,1:jpj)
468	vn_e (1:jpi,1:jpj) = vn_b(1:jpi,1:jpj)
469	!
470	hu_e (1:jpi,1:jpj) = hu_n(1:jpi,1:jpj)
471	hv_e (1:jpi,1:jpj) = hv_n(1:jpi,1:jpj)
472	hur_e (1:jpi,1:jpj) = r1_hu_n(1:jpi,1:jpj)
473	hvr_e (1:jpi,1:jpj) = r1_hv_n(1:jpi,1:jpj)
474	ELSE ! CENTRED integration: start from BEFORE fields
475	sshn_e(1:jpi,1:jpj) = sshb(1:jpi,1:jpj)
476	un_e (1:jpi,1:jpj) = ub_b(1:jpi,1:jpj)
477	vn_e (1:jpi,1:jpj) = vb_b(1:jpi,1:jpj)
478	!
479	hu_e (1:jpi,1:jpj) = hu_b(1:jpi,1:jpj)
480	hv_e (1:jpi,1:jpj) = hv_b(1:jpi,1:jpj)
481	hur_e (1:jpi,1:jpj) = r1_hu_b(1:jpi,1:jpj)
482	hvr_e (1:jpi,1:jpj) = r1_hv_b(1:jpi,1:jpj)
483	ENDIF
484	!
485	hu_n_xtd(1:jpi,1:jpj) = hu_n(1:jpi,1:jpj)
486	hv_n_xtd(1:jpi,1:jpj) = hv_n(1:jpi,1:jpj)
487	!
488	!
489	! ! Extend arrays
490	! ! --------------
491	!
492	IF( ln_linssh ) THEN
493	CALL lbc_lnk_multi( 'dynspg_ts', hu_n_xtd, 'U', -1._wp, hv_n_xtd, 'V', -1._wp &
494	& , zCdU_u , 'U', -1._wp, zCdU_v , 'V', -1._wp &
495	& , zu_frc , 'U', -1._wp, zv_frc , 'V', -1._wp &
496	& , un_e , 'U', -1._wp, vn_e , 'V', -1._wp &
497	& , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp &
498	& , hur_e , 'U', -1._wp, hvr_e , 'V', -1._wp &
499	& , zhtp2_e , 'T', 1._wp, zhup2_e, 'U', -1._wp, zhvp2_e , 'V', -1._wp &
500	& , zssh_frc, 'T', 1._wp, sshn_e , 'T', 1._wp, khlcom = nn_hls+nn_hlts )
501	ELSE
502	CALL lbc_lnk_multi( 'dynspg_ts', hu_n_xtd, 'U', -1._wp, hv_n_xtd, 'V', -1._wp &
503	& , zCdU_u , 'U', -1._wp, zCdU_v , 'V', -1._wp &
504	& , zu_frc , 'U', -1._wp, zv_frc , 'V', -1._wp &
505	& , un_e , 'U', -1._wp, vn_e , 'V', -1._wp &
506	& , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp &
507	& , hur_e , 'U', -1._wp, hvr_e , 'V', -1._wp &
508	& , zssh_frc, 'T', 1._wp, sshn_e , 'T', 1._wp, khlcom = nn_hls+nn_hlts )
509	END IF
510	!
511	! Initialize sums:
512	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
513	va_b (:,:) = 0._wp
514	ssha (:,:) = 0._wp ! Sum for after averaged sea level
515	un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop
516	vn_adv(:,:) = 0._wp
517	!
518	IF( ln_wd_dl ) THEN
519	zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary)
520	zvwdmask(:,:) = 0._wp !
521	zuwdav2 (:,:) = 0._wp
522	zvwdav2 (:,:) = 0._wp
523	END IF
524
525	ixtd = nn_hls + nn_hlts ! solution is now correct over the whole domain (interior + regular halos + time splitting halos)
526	ibi = 1 - nn_hlts ; ibj = 1 - nn_hlts
527	iei = jpi + nn_hlts ; iej = jpj + nn_hlts
528	! ! ==================== !
529	DO jm = 1, icycle ! sub-time-step loop !
530	! ! ==================== !
531	!
532	l_full_nf_update = jm == icycle ! false: disable full North fold update (performances) for jm = 1 to icycle-1
533	!
534	! !== Update the forcing ==! (BDY and tides)
535	!
536	IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jm, kt_offset= noffset+1 )
537	IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jm, kt_offset= noffset )
538	!
539	! !== extrapolation at mid-step ==! (jm+1/2)
540	!
541	! !* Set extrapolation coefficients for predictor step:
542	IF( (jm<3) .AND. ll_init ) THEN ! Forward
543	za1 = 1._wp
544	za2 = 0._wp
545	za3 = 0._wp
546	ELSE ! AB3-AM4 Coefficients: bet=0.281105
547	za1 = 1.781105_wp ! za1 = 3/2 + bet
548	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
549	za3 = 0.281105_wp ! za3 = bet
550	ENDIF
551	!
552	! !* Extrapolate barotropic velocities at mid-step (jm+1/2)
553	!-- m+1/2 m m-1 m-2 --!
554	!-- u = (3/2+beta) u -(1/2+2beta) u + beta u --!
555	!-------------------------------------------------------------------------!
556	ua_e(ibi:iei,ibj:iej) = za1 * un_e(ibi:iei,ibj:iej) + za2 * ub_e(ibi:iei,ibj:iej) + za3 * ubb_e(ibi:iei,ibj:iej)
557	va_e(ibi:iei,ibj:iej) = za1 * vn_e(ibi:iei,ibj:iej) + za2 * vb_e(ibi:iei,ibj:iej) + za3 * vbb_e(ibi:iei,ibj:iej)
558
559	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
560	! ! ------------------
561	! !* Extrapolate Sea Level at mid-step (jm+1/2)
562	!-- m+1/2 m m-1 m-2 --!
563	!-- ssh = (3/2+beta) ssh -(1/2+2beta) ssh + beta ssh --!
564	!--------------------------------------------------------------------------------!
565	zsshp2_e(ibi:iei,ibj:iej) = za1 * sshn_e (ibi:iei,ibj:iej) + za2 * sshb_e(ibi:iei,ibj:iej) &
566	& + za3 * sshbb_e(ibi:iei,ibj:iej)
567	! set wetting & drying mask at tracer points for this barotropic mid-step
568	IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask )
569	!
570	! ! ocean t-depth at mid-step
571	zhtp2_e(ibi:iei,ibj:iej) = ht_0_xtd(ibi:iei,ibj:iej) + zsshp2_e(ibi:iei,ibj:iej)
572	!
573	! ! ocean u- and v-depth at mid-step
574	DO jj = ibj, iej-1 ! not last column, not last row
575	DO ji = ibi, iei-1
576	zhup2_e(ji,jj) = hu_0_xtd(ji,jj) + r1_2 * r1_e1e2u_xtd(ji,jj) &
577	& * ( e1e2t_xtd(ji ,jj) * zsshp2_e(ji ,jj) &
578	& + e1e2t_xtd(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask_xtd(ji,jj)
579	zhvp2_e(ji,jj) = hv_0_xtd(ji,jj) + r1_2 * r1_e1e2v_xtd(ji,jj) &
580	& * ( e1e2t_xtd(ji,jj ) * zsshp2_e(ji,jj ) &
581	& + e1e2t_xtd(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask_xtd(ji,jj)
582	END DO
583	END DO
584	!
585	ENDIF
586	!
587	! !== after SSH ==! (jm+1)
588	!
589	! ! update (ua_e,va_e) to enforce volume conservation at open boundaries
590	! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
591	IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jm, ua_e, va_e, zhup2_e, zhvp2_e )
592	!
593	! ! resulting flux at mid-step (not over the full domain)
594	zhU(ibi:iei-1,ibj:iej-1) = e2u_xtd(ibi:iei-1,ibj:iej-1) * ua_e(ibi:iei-1,ibj:iej-1) * zhup2_e(ibi:iei-1,ibj:iej-1)
595	zhV(ibi:iei-1,ibj:iej-1) = e1v_xtd(ibi:iei-1,ibj:iej-1) * va_e(ibi:iei-1,ibj:iej-1) * zhvp2_e(ibi:iei-1,ibj:iej-1)
596	!
597	#if defined key_agrif
598	! Set fluxes during predictor step to ensure volume conservation
599	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
600	IF((nbondi == -1).OR.(nbondi == 2)) THEN
601	DO jj = 1, jpj
602	zhU(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)
603	zhV(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)
604	END DO
605	ENDIF
606	IF((nbondi == 1).OR.(nbondi == 2)) THEN
607	DO jj=1,jpj
608	zhU(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
609	zhV(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj)
610	END DO
611	ENDIF
612	IF((nbondj == -1).OR.(nbondj == 2)) THEN
613	DO ji=1,jpi
614	zhV(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)
615	zhU(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)
616	END DO
617	ENDIF
618	IF((nbondj == 1).OR.(nbondj == 2)) THEN
619	DO ji=1,jpi
620	zhV(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
621	zhU(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1)
622	END DO
623	ENDIF
624	ENDIF
625	#endif
626	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
627
628	IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where
629	! ! the direction of the flow is from dry cells
630	CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V
631	!
632	ENDIF
633	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True)
634	IF ( ln_wd_dl_bc ) THEN
635	zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column
636	zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row
637	END IF
638	!
639	! Sum over sub-time-steps to compute advective velocities (only correct on interior domain)
640	za2 = wgtbtp2(jm) ! zhU, zhV hold fluxes extrapolated at jm+1/2
641	un_adv(1:jpi,1:jpj) = un_adv(1:jpi,1:jpj) + za2 * zhU(1:jpi,1:jpj) * r1_e2u(1:jpi,1:jpj)
642	vn_adv(1:jpi,1:jpj) = vn_adv(1:jpi,1:jpj) + za2 * zhV(1:jpi,1:jpj) * r1_e1v(1:jpi,1:jpj)
643	!
644	!
645	! Compute Sea Level at step jit+1
646	!-- m+1 m m+1/2 --!
647	!-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --!
648	!-------------------------------------------------------------------------!
649	! correct domain reduction
650	ixtd = ixtd - 1
651	ibi = ibi + 1 ; ibj = ibj + 1
652	iei = iei - 1 ; iej = iej - 1
653	DO jj = ibj, iej
654	DO ji = ibi, iei
655	zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t_xtd(ji,jj)
656	ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask_xtd(ji,jj)
657	END DO
658	END DO
659	!
660	IF( nn_hlts == 0 ) THEN
661	CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )
662	ixtd = nn_hls + nn_hlts ! solution is now correct over the whole domain
663	ibi = 1 - nn_hlts ; ibj = 1 - nn_hlts
664	iei = jpi + nn_hlts ; iej = jpj + nn_hlts
665	END IF
666
667	!
668	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
669	IF( ln_bdy ) THEN
670	CALL swap_bdyptr ! bdy treatment is now done on extended domain
671	CALL bdy_ssh( ssha_e, idbi, idei, idbj, idej, ldcomall=.true., pmask=ssmask_xtd, khlcom=nn_hls+nn_hlts )
672	CALL swap_bdyptr ! bdy treatment is now done on regular domain
673	END IF
674
675	#if defined key_agrif
676	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jm )
677	#endif
678	!
679	! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
680	!-- m+1/2 m+1 m m-1 m-2 --!
681	!-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --!
682	!-----------------------------------------------------------------------------------------!
683	CALL ts_bck_interp( jm, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation
684	zsshp2_e(ibi:iei,ibj:iej) = za0 * ssha_e(ibi:iei,ibj:iej) + za1 * sshn_e (ibi:iei,ibj:iej) &
685	& + za2 * sshb_e(ibi:iei,ibj:iej) + za3 * sshbb_e(ibi:iei,ibj:iej)
686	!
687	!
688	! Sea Surface Height at u-,v-points (vvl case only)
689	IF( .NOT.ln_linssh ) THEN
690	DO jj = ibj, iej-1
691	DO ji = ibi, iei-1
692	zsshu_a(ji,jj) = r1_2 * ssumask_xtd(ji,jj) * r1_e1e2u_xtd(ji,jj) &
693	& * ( e1e2t_xtd(ji ,jj ) * ssha_e(ji ,jj ) &
694	& + e1e2t_xtd(ji+1,jj ) * ssha_e(ji+1,jj ) )
695	zsshv_a(ji,jj) = r1_2 * ssvmask_xtd(ji,jj) * r1_e1e2v_xtd(ji,jj) &
696	& * ( e1e2t_xtd(ji ,jj ) * ssha_e(ji ,jj ) &
697	& + e1e2t_xtd(ji ,jj+1) * ssha_e(ji ,jj+1) )
698	END DO
699	END DO
700	ENDIF
701	!
702	! ! Surface pressure gradient
703	zldg = ( 1._wp - rn_scal_load ) * grav ! local factor
704	DO jj = ibj, iej-1
705	DO ji = ibi, iei-1
706	zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u_xtd(ji,jj)
707	zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v_xtd(ji,jj)
708	END DO
709	END DO
710	IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient
711	CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters
712	zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
713	zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
714	ENDIF
715	!
716	! Add Coriolis trend:
717	! - zwz array used in dyn_cor_2d or triads normally depend on sea level with ln_linssh=F and should be updated
718	! at each time step. We however keep them constant here for optimization.
719	! - Recall that zhU and zhV hold fluxes at jm+1/2 (extrapolated not backward interpolated)
720	! - zu_trd_xtd and zv_trd_xtd are only correct on (ibi+1:iei-1,ibj+1:iej-1)
721	! NOTE : input flux arguments have to be correct (ibi:iei,ibj:iej) -> a lbc call between input arguments computation
722	! and this call without fluxes (typically after ssh at step m+1 computation) would not yield correct results
723	CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV , idbi, idei, idbj, idej &
724	& , zu_trd, zv_trd , idbi, idei, idbj, idej )
725	!
726	! Add tidal astronomical forcing if defined
727	! pot_astro is correct on 1:jpi,1:jpj
728	IF ( ln_tide .AND. ln_tide_pot ) THEN
729	DO jj = 2, jpjm1
730	DO ji = fs_2, fs_jpim1 ! vector opt.
731	zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
732	zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
733	END DO
734	END DO
735	ENDIF
736	!
737	! Add bottom stresses:
738	!jth do implicitly instead
739	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
740	DO jj = ibj, iej
741	DO ji = ibi, iei
742	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
743	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
744	END DO
745	END DO
746	ENDIF
747	!
748	! Set next velocities:
749	! Compute barotropic speeds at step jit+1 (h : total height of the water colomn)
750	!-- VECTOR FORM
751	!-- m+1 m / m+1/2 \ --!
752	!-- u = u + delta_t' * \ (1-r)g grad_x( ssh') - f * k vect u + frc / --!
753	!-- --!
754	!-- FLUX FORM --!
755	!-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --!
756	!-- u = m+1 \| h * u + delta_t' * \ h * (1-r)g grad_x( ssh') - h * f * k vect u + h * frc / \| --!
757	!-- h \ / --!
758	!------------------------------------------------------------------------------------------------------------------------!
759	! correct domain reduction
760	ixtd = ixtd - 1
761	ibi = ibi + 1 ; ibj = ibj + 1
762	iei = iei - 1 ; iej = iej - 1
763	!
764	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
765	DO jj = ibj, iej
766	DO ji = ibi, iei
767	ua_e(ji,jj) = ( un_e(ji,jj) &
768	& + rdtbt * ( zu_spg(ji,jj) &
769	& + zu_trd(ji,jj) &
770	& + zu_frc(ji,jj) ) &
771	& ) * ssumask_xtd(ji,jj)
772
773	va_e(ji,jj) = ( vn_e(ji,jj) &
774	& + rdtbt * ( zv_spg(ji,jj) &
775	& + zv_trd(ji,jj) &
776	& + zv_frc(ji,jj) ) &
777	& ) * ssvmask_xtd(ji,jj)
778	END DO
779	END DO
780	!
781	ELSE !* Flux form
782	DO jj = ibj, iej
783	DO ji = ibi, iei
784	! ! hu_e, hv_e hold depth at jm, zhup2_e, zhvp2_e hold extrapolated depth at jm+1/2
785	! ! backward interpolated depth used in spg terms at jm+1/2
786	zhu_bck = hu_0_xtd(ji,jj) + r1_2r1_e1e2u_xtd(ji,jj) &
787	& ( e1e2t_xtd(ji ,jj) * zsshp2_e(ji ,jj) &
788	& + e1e2t_xtd(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask_xtd(ji,jj)
789	zhv_bck = hv_0_xtd(ji,jj) + r1_2r1_e1e2v_xtd(ji,jj) &
790	& ( e1e2t_xtd(ji,jj ) * zsshp2_e(ji,jj ) &
791	& + e1e2t_xtd(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask_xtd(ji,jj)
792	! ! inverse depth at jm+1
793	z1_hu = ssumask_xtd(ji,jj) / ( hu_0_xtd(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask_xtd(ji,jj) )
794	z1_hv = ssvmask_xtd(ji,jj) / ( hv_0_xtd(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask_xtd(ji,jj) )
795	!
796	ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) &
797	& + rdtbt * ( zhu_bck * zu_spg(ji,jj) & !
798	& + zhup2_e(ji,jj) * zu_trd(ji,jj) & !
799	& + hu_n_xtd (ji,jj) * zu_frc(ji,jj) ) ) * z1_hu
800	!
801	va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) &
802	& + rdtbt * ( zhv_bck * zv_spg(ji,jj) & !
803	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) & !
804	& + hv_n_xtd (ji,jj) * zv_frc(ji,jj) ) ) * z1_hv
805	END DO
806	END DO
807	ENDIF
808	!jth implicit bottom friction:
809	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
810	DO jj = 2, jpjm1
811	DO ji = fs_2, fs_jpim1 ! vector opt.
812	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
813	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
814	END DO
815	END DO
816	ENDIF
817
818	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) and inverse depth
819	hu_e (ibi:iei,ibj:iej) = hu_0_xtd(ibi:iei,ibj:iej) + zsshu_a(ibi:iei,ibj:iej)
820	hv_e (ibi:iei,ibj:iej) = hv_0_xtd(ibi:iei,ibj:iej) + zsshv_a(ibi:iei,ibj:iej)
821	hur_e(ibi:iei,ibj:iej) = ssumask_xtd(ibi:iei,ibj:iej) / ( hu_e(ibi:iei,ibj:iej) + 1._wp - ssumask_xtd(ibi:iei,ibj:iej) )
822	hvr_e(ibi:iei,ibj:iej) = ssvmask_xtd(ibi:iei,ibj:iej) / ( hv_e(ibi:iei,ibj:iej) + 1._wp - ssvmask_xtd(ibi:iei,ibj:iej) )
823	ENDIF
824
825	IF( ixtd == 0 ) THEN
826	IF( .NOT. ln_linssh ) THEN
827	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & ! after
828	& , un_e , 'U', -1._wp, vn_e , 'V', -1._wp & ! now
829	& , ub_e , 'U', -1._wp, vb_e , 'V', -1._wp & ! before
830	& , ubb_e , 'U', -1._wp, vbb_e , 'V', -1._wp & ! before before
831	& , ssha_e, 'T', 1._wp, sshn_e , 'T', 1._wp & ! after, now
832	& , sshb_e, 'T', 1._wp, sshbb_e, 'T', 1._wp & ! before, before before
833	& , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp &
834	& , hur_e , 'U', -1._wp, hvr_e , 'V', -1._wp &
835	& , khlcom = nn_hls+nn_hlts )
836	ELSE
837	CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & ! after
838	& , un_e , 'U', -1._wp, vn_e , 'V', -1._wp & ! now
839	& , ub_e , 'U', -1._wp, vb_e , 'V', -1._wp & ! before
840	& , ubb_e , 'U', -1._wp, vbb_e , 'V', -1._wp & ! before before
841	& , ssha_e, 'T', 1._wp, sshn_e , 'T', 1._wp & ! after, now
842	& , sshb_e, 'T', 1._wp, sshbb_e, 'T', 1._wp & ! before, before before
843	& , khlcom = nn_hls+nn_hlts )
844	END IF
845	ixtd = nn_hls + nn_hlts ! solution is now correct over the whole domain
846	ibi = 1 - nn_hlts ; ibj = 1 - nn_hlts
847	iei = jpi + nn_hlts ; iej = jpj + nn_hlts
848	END IF
849	!
850	!
851	! ! open boundaries
852	! ! bdy treatment is here done on regular domain (nn_hlts forced to 1 if ln_bdy or ln_tides)
853	IF( ln_bdy ) CALL bdy_dyn2d( jm, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e, idbi, idei, idbj, idej &
854	& , ldcomall=.true., pumask=ssumask_xtd, pvmask=ssvmask_xtd, khlcom=nn_hls+nn_hlts )
855	#if defined key_agrif
856	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jm ) ! Agrif
857	#endif
858	! !* Swap
859	! ! ----
860	ubb_e (:,:) = ub_e (:,:)
861	ub_e (:,:) = un_e (:,:)
862	un_e (:,:) = ua_e (:,:)
863	!
864	vbb_e (:,:) = vb_e (:,:)
865	vb_e (:,:) = vn_e (:,:)
866	vn_e (:,:) = va_e (:,:)
867	!
868	sshbb_e(:,:) = sshb_e(:,:)
869	sshb_e (:,:) = sshn_e(:,:)
870	sshn_e (:,:) = ssha_e(:,:)
871
872	! !* Sum over whole bt loop
873	! ! ----------------------
874	za1 = wgtbtp1(jm)
875	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
876	ua_b(1:jpi,1:jpj) = ua_b(1:jpi,1:jpj) + za1 * ua_e(1:jpi,1:jpj)
877	va_b(1:jpi,1:jpj) = va_b(1:jpi,1:jpj) + za1 * va_e(1:jpi,1:jpj)
878	ELSE ! Sum transports
879	IF ( .NOT.ln_wd_dl ) THEN
880	ua_b(1:jpi,1:jpj) = ua_b(1:jpi,1:jpj) + za1 * ua_e(1:jpi,1:jpj) * hu_e(1:jpi,1:jpj)
881	va_b(1:jpi,1:jpj) = va_b(1:jpi,1:jpj) + za1 * va_e(1:jpi,1:jpj) * hv_e(1:jpi,1:jpj)
882	ELSE
883	ua_b(1:jpi,1:jpj) = ua_b(1:jpi,1:jpj) + za1 * ua_e(1:jpi,1:jpj) * hu_e(1:jpi,1:jpj) * zuwdmask(1:jpi,1:jpj)
884	va_b(1:jpi,1:jpj) = va_b(1:jpi,1:jpj) + za1 * va_e(1:jpi,1:jpj) * hv_e(1:jpi,1:jpj) * zvwdmask(1:jpi,1:jpj)
885	END IF
886	ENDIF
887	! ! Sum sea level
888	ssha(1:jpi,1:jpj) = ssha(1:jpi,1:jpj) + za1 * ssha_e(1:jpi,1:jpj)
889
890	! ! ==================== !
891	END DO ! end loop !
892	! ! ==================== !
893	! -----------------------------------------------------------------------------
894	! Phase 3. update the general trend with the barotropic trend
895	! -----------------------------------------------------------------------------
896	! Correction on regular halos
897	CALL lbc_lnk_multi( 'dynspg_ts', un_adv, 'U', -1._wp, vn_adv, 'V', -1._wp &
898	& , ua_b , 'U', -1._wp, va_b , 'V', -1._wp &
899	& , ssha , 'T', -1._wp )
900	!
901	! Set advection velocity correction:
902	IF (ln_bt_fw) THEN
903	IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
904	DO jj = 1, jpj
905	DO ji = 1, jpi
906	zun_save = un_adv(ji,jj)
907	zvn_save = vn_adv(ji,jj)
908	! ! apply the previously computed correction
909	un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) )
910	vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) )
911	! ! Update corrective fluxes for next time step
912	un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
913	vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
914	! ! Save integrated transport for next computation
915	ub2_b(ji,jj) = zun_save
916	vb2_b(ji,jj) = zvn_save
917	END DO
918	END DO
919	ELSE
920	un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero
921	vn_bf(:,:) = 0._wp
922	ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation
923	vb2_b(:,:) = vn_adv(:,:)
924	END IF
925	ENDIF
926
927
928	!
929	! Update barotropic trend:
930	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
931	DO jk=1,jpkm1
932	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
933	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
934	END DO
935	ELSE
936	! At this stage, ssha has been corrected: compute new depths at velocity points
937	DO jj = 1, jpjm1
938	DO ji = 1, jpim1 ! NO Vector Opt.
939	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
940	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
941	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
942	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
943	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
944	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
945	END DO
946	END DO ! Boundary conditions
947	CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp, khlcom=nn_hls+nn_hlts ) ! change array used?
948	!
949	DO jk=1,jpkm1
950	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
951	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
952	END DO
953	! Save barotropic velocities not transport:
954	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(1:jpi,1:jpj) + 1._wp - ssumask(:,:) )
955	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(1:jpi,1:jpj) + 1._wp - ssvmask(:,:) )
956	ENDIF
957
958
959	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
960	DO jk = 1, jpkm1
961	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)r1_hu_n(:,:) - un_b(:,:) ) umask(:,:,jk)
962	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)r1_hv_n(:,:) - vn_b(:,:) ) vmask(:,:,jk)
963	END DO
964
965	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
966	DO jk = 1, jpkm1
967	un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
968	& + zuwdav2(:,:)(un(:,:,jk) - un_adv(:,:)r1_hu_n(:,:)) ) * umask(:,:,jk)
969	vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) &
970	& + zvwdav2(:,:)(vn(:,:,jk) - vn_adv(:,:)r1_hv_n(:,:)) ) * vmask(:,:,jk)
971	END DO
972	END IF
973
974
975	CALL iom_put( "ubar", un_adv(:,:)*r1_hu_n(:,:) ) ! barotropic i-current
976	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv_n(:,:) ) ! barotropic i-current
977	!
978	#if defined key_agrif
979	! Save time integrated fluxes during child grid integration
980	! (used to update coarse grid transports at next time step)
981	!
982	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
983	IF( Agrif_NbStepint() == 0 ) THEN
984	ub2_i_b(:,:) = 0._wp
985	vb2_i_b(:,:) = 0._wp
986	END IF
987	!
988	za1 = 1._wp / REAL(Agrif_rhot(), wp)
989	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
990	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
991	ENDIF
992	#endif
993	! !* write time-spliting arrays in the restart
994	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
995	!
996	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
997	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
998	!
999	IF( ln_diatmb ) THEN
1000	CALL iom_put( "baro_u" , un_bssumask(:,:)+zmdi(1.-ssumask(:,:) ) ) ! Barotropic U Velocity
1001	CALL iom_put( "baro_v" , vn_bssvmask(:,:)+zmdi(1.-ssvmask(:,:) ) ) ! Barotropic V Velocity
1002	ENDIF
1003	!
1004
1005	! deallocate temporary arrays
1006	DEALLOCATE( zu_trd , zv_trd &
1007	& , zu_frc , zv_frc &
1008	& , zu_spg , zv_spg &
1009	& , zsshu_a , zsshv_a &
1010	& , zhup2_e , zhvp2_e &
1011	& , zCdU_u , zCdU_v &
1012	& , zhU , zhV &
1013	& , zssh_frc, zsshp2_e &
1014	& , zhtp2_e &
1015	& , hu_n_xtd, hv_n_xtd )
1016	!
1017	END SUBROUTINE dyn_spg_ts
1018
1019
1020	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1021	!!---------------------------------------------------------------------
1022	!! * ROUTINE ts_wgt *
1023	!!
1024	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1025	!!----------------------------------------------------------------------
1026	LOGICAL, INTENT(in ) :: ll_av ! temporal averaging=.true.
1027	LOGICAL, INTENT(in ) :: ll_fw ! forward time splitting =.true.
1028	INTEGER, INTENT(inout) :: jpit ! cycle length
1029	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
1030	zwgt2 ! Secondary weights
1031
1032	INTEGER :: jic, jm, ji ! temporary integers
1033	REAL(wp) :: za1, za2
1034	!!----------------------------------------------------------------------
1035
1036	zwgt1(:) = 0._wp
1037	zwgt2(:) = 0._wp
1038
1039	! Set time index when averaged value is requested
1040	IF (ll_fw) THEN
1041	jic = nn_baro
1042	ELSE
1043	jic = 2 * nn_baro
1044	ENDIF
1045
1046	! Set primary weights:
1047	IF (ll_av) THEN
1048	! Define simple boxcar window for primary weights
1049	! (width = nn_baro, centered around jic)
1050	SELECT CASE ( nn_bt_flt )
1051	CASE( 0 ) ! No averaging
1052	zwgt1(jic) = 1._wp
1053	jpit = jic
1054
1055	CASE( 1 ) ! Boxcar, width = nn_baro
1056	DO jm = 1, 3*nn_baro
1057	za1 = ABS(float(jm-jic))/float(nn_baro)
1058	IF (za1 < 0.5_wp) THEN
1059	zwgt1(jm) = 1._wp
1060	jpit = jm
1061	ENDIF
1062	ENDDO
1063
1064	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
1065	DO jm = 1, 3*nn_baro
1066	za1 = ABS(float(jm-jic))/float(nn_baro)
1067	IF (za1 < 1._wp) THEN
1068	zwgt1(jm) = 1._wp
1069	jpit = jm
1070	ENDIF
1071	ENDDO
1072	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1073	END SELECT
1074
1075	ELSE ! No time averaging
1076	zwgt1(jic) = 1._wp
1077	jpit = jic
1078	ENDIF
1079
1080	! Set secondary weights
1081	DO jm = 1, jpit
1082	DO ji = jm, jpit
1083	zwgt2(jm) = zwgt2(jm) + zwgt1(ji)
1084	END DO
1085	END DO
1086
1087	! Normalize weigths:
1088	za1 = 1._wp / SUM(zwgt1(1:jpit))
1089	za2 = 1._wp / SUM(zwgt2(1:jpit))
1090	DO jm = 1, jpit
1091	zwgt1(jm) = zwgt1(jm) * za1
1092	zwgt2(jm) = zwgt2(jm) * za2
1093	END DO
1094	!
1095	END SUBROUTINE ts_wgt
1096
1097
1098	SUBROUTINE ts_rst( kt, cdrw )
1099	!!---------------------------------------------------------------------
1100	!! * ROUTINE ts_rst *
1101	!!
1102	!! ** Purpose : Read or write time-splitting arrays in restart file
1103	!!----------------------------------------------------------------------
1104	INTEGER , INTENT(in) :: kt ! ocean time-step
1105	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1106	!!----------------------------------------------------------------------
1107	!
1108	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
1109	! ! ---------------
1110	IF( ln_rstart .AND. ln_bt_fw .AND. (neuler/=0) ) THEN !* Read the restart file
1111	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios )
1112	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios )
1113	CALL iom_get( numror, jpdom_autoglo, 'un_bf' , un_bf (:,:), ldxios = lrxios )
1114	CALL iom_get( numror, jpdom_autoglo, 'vn_bf' , vn_bf (:,:), ldxios = lrxios )
1115	IF( .NOT.ln_bt_av ) THEN
1116	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:), ldxios = lrxios )
1117	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:), ldxios = lrxios )
1118	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:), ldxios = lrxios )
1119	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:), ldxios = lrxios )
1120	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:), ldxios = lrxios )
1121	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:), ldxios = lrxios )
1122	ENDIF
1123	#if defined key_agrif
1124	! Read time integrated fluxes
1125	IF ( .NOT.Agrif_Root() ) THEN
1126	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lrxios )
1127	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lrxios )
1128	ENDIF
1129	#endif
1130	ELSE !* Start from rest
1131	IF(lwp) WRITE(numout,*)
1132	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
1133	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
1134	un_adv(:,:) = 0._wp ; vn_adv(:,:) = 0._wp ! used in the 1st interpol of agrif
1135	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
1136	#if defined key_agrif
1137	IF ( .NOT.Agrif_Root() ) THEN
1138	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
1139	ENDIF
1140	#endif
1141	ENDIF
1142	!
1143	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
1144	! ! -------------------
1145	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
1146	IF( lwxios ) CALL iom_swap( cwxios_context )
1147	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios )
1148	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios )
1149	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios )
1150	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios )
1151	!
1152	IF (.NOT.ln_bt_av) THEN
1153	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:), ldxios = lwxios )
1154	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:), ldxios = lwxios )
1155	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:), ldxios = lwxios )
1156	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:), ldxios = lwxios )
1157	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:), ldxios = lwxios )
1158	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:), ldxios = lwxios )
1159	ENDIF
1160	#if defined key_agrif
1161	! Save time integrated fluxes
1162	IF ( .NOT.Agrif_Root() ) THEN
1163	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lwxios )
1164	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lwxios )
1165	ENDIF
1166	#endif
1167	IF( lwxios ) CALL iom_swap( cxios_context )
1168	ENDIF
1169	!
1170	END SUBROUTINE ts_rst
1171
1172
1173	SUBROUTINE dyn_spg_ts_init
1174	!!---------------------------------------------------------------------
1175	!! * ROUTINE dyn_spg_ts_init *
1176	!!
1177	!! ** Purpose : Set time splitting options
1178	!!----------------------------------------------------------------------
1179	INTEGER :: ji ,jj ! dummy loop indices
1180	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
1181	REAL(wp), DIMENSION(jpi,jpj) :: zcu
1182	INTEGER :: inum
1183	!!----------------------------------------------------------------------
1184	!
1185	! Max courant number for ext. grav. waves
1186	!
1187	DO jj = 1, jpj
1188	DO ji =1, jpi
1189	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1190	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1191	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1192	END DO
1193	END DO
1194	!
1195	zcmax = MAXVAL( zcu(:,:) )
1196	CALL mpp_max( 'dynspg_ts', zcmax )
1197
1198	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1199	IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1200
1201	rdtbt = rdt / REAL( nn_baro , wp )
1202	zcmax = zcmax * rdtbt
1203	! Print results
1204	IF(lwp) WRITE(numout,*)
1205	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1206	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1207	IF( ln_bt_auto ) THEN
1208	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_baro '
1209	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1210	ELSE
1211	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist nn_baro = ', nn_baro
1212	ENDIF
1213
1214	IF(ln_bt_av) THEN
1215	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_baro time steps is on '
1216	ELSE
1217	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
1218	ENDIF
1219	!
1220	!
1221	IF(ln_bt_fw) THEN
1222	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1223	ELSE
1224	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1225	ENDIF
1226	!
1227	#if defined key_agrif
1228	! Restrict the use of Agrif to the forward case only
1229	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1230	#endif
1231	!
1232	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1233	SELECT CASE ( nn_bt_flt )
1234	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1235	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1236	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1237	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1238	END SELECT
1239	!
1240	IF(lwp) WRITE(numout,*) ' '
1241	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1242	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1243	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1244	!
1245	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1246	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1247	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1248	ENDIF
1249	!
1250	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1251	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1252	ENDIF
1253	IF( zcmax>0.9_wp ) THEN
1254	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1255	ENDIF
1256	!
1257	! ! Allocate time-splitting arrays
1258	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
1259	!
1260	! ! read restart when needed
1261	CALL ts_rst( nit000, 'READ' )
1262	!
1263	IF( lwxios ) THEN
1264	! define variables in restart file when writing with XIOS
1265	CALL iom_set_rstw_var_active('ub2_b')
1266	CALL iom_set_rstw_var_active('vb2_b')
1267	CALL iom_set_rstw_var_active('un_bf')
1268	CALL iom_set_rstw_var_active('vn_bf')
1269	!
1270	IF (.NOT.ln_bt_av) THEN
1271	CALL iom_set_rstw_var_active('sshbb_e')
1272	CALL iom_set_rstw_var_active('ubb_e')
1273	CALL iom_set_rstw_var_active('vbb_e')
1274	CALL iom_set_rstw_var_active('sshb_e')
1275	CALL iom_set_rstw_var_active('ub_e')
1276	CALL iom_set_rstw_var_active('vb_e')
1277	ENDIF
1278	#if defined key_agrif
1279	! Save time integrated fluxes
1280	IF ( .NOT.Agrif_Root() ) THEN
1281	CALL iom_set_rstw_var_active('ub2_i_b')
1282	CALL iom_set_rstw_var_active('vb2_i_b')
1283	ENDIF
1284	#endif
1285	ENDIF
1286	!
1287	! initialize extended scale factors
1288	ht_0_xtd (1:jpi,1:jpj) = ht_0 (1:jpi,1:jpj)
1289	hu_0_xtd (1:jpi,1:jpj) = hu_0 (1:jpi,1:jpj)
1290	hv_0_xtd (1:jpi,1:jpj) = hv_0 (1:jpi,1:jpj)
1291	r1_e1e2t_xtd(1:jpi,1:jpj) = r1_e1e2t(1:jpi,1:jpj)
1292	r1_e1e2u_xtd(1:jpi,1:jpj) = r1_e1e2u(1:jpi,1:jpj)
1293	r1_e1e2v_xtd(1:jpi,1:jpj) = r1_e1e2v(1:jpi,1:jpj)
1294	e1e2t_xtd (1:jpi,1:jpj) = e1e2t (1:jpi,1:jpj)
1295	ssmask_xtd (1:jpi,1:jpj) = ssmask (1:jpi,1:jpj)
1296	ssumask_xtd (1:jpi,1:jpj) = ssumask (1:jpi,1:jpj)
1297	ssvmask_xtd (1:jpi,1:jpj) = ssvmask (1:jpi,1:jpj)
1298	e2u_xtd (1:jpi,1:jpj) = e2u (1:jpi,1:jpj)
1299	e1v_xtd (1:jpi,1:jpj) = e1v (1:jpi,1:jpj)
1300	r1_e1u_xtd (1:jpi,1:jpj) = r1_e1u (1:jpi,1:jpj)
1301	r1_e2v_xtd (1:jpi,1:jpj) = r1_e2v (1:jpi,1:jpj)
1302	!
1303	CALL lbc_lnk_multi( 'dynspg_ts', ht_0_xtd , 'T', 1._wp, hu_0_xtd , 'U', -1._wp, hv_0_xtd , 'V', -1._wp &
1304	& , r1_e1e2t_xtd, 'T', 1._wp, r1_e1e2u_xtd, 'U', -1._wp, r1_e1e2v_xtd, 'V', -1._wp &
1305	& , ssmask_xtd , 'T', 1._wp, ssumask_xtd , 'U', -1._wp, ssvmask_xtd , 'V', -1._wp &
1306	& , e1e2t_xtd , 'T', 1._wp, e2u_xtd , 'U', -1._wp, e1v_xtd , 'V', -1._wp &
1307	& , r1_e1u_xtd , 'U', -1._wp, r1_e2v_xtd , 'V', -1._wp &
1308	& , khlcom = nn_hls+nn_hlts )
1309	IF( ln_dynvor_enT ) THEN
1310	ff_t_xtd (1:jpi,1:jpj) = ff_t (1:jpi,1:jpj)
1311	CALL lbc_lnk_multi( 'dynspg_ts', ff_t_xtd , 'F', -1._wp, khlcom = nn_hls+nn_hlts )
1312	END IF
1313
1314	!
1315	END SUBROUTINE dyn_spg_ts_init
1316
1317
1318	SUBROUTINE dyn_cor_2d_init( kdbi, kdei, kdbj, kdej )
1319	!!---------------------------------------------------------------------
1320	!! * ROUTINE dyn_cor_2d_init *
1321	!!
1322	!! ** Purpose : Set time splitting options
1323	!! Set arrays to remove/compute coriolis trend.
1324	!! Do it once during initialization if volume is fixed, else at each long time step.
1325	!! Note that these arrays are also used during barotropic loop. These are however frozen
1326	!! although they should be updated in the variable volume case. Not a big approximation.
1327	!! To remove this approximation, copy lines below inside barotropic loop
1328	!! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
1329	!!
1330	!! Compute zwz = f / ( height of the water colomn )
1331	!!----------------------------------------------------------------------
1332	INTEGER , INTENT(in ) :: kdbi, kdei, kdbj, kdej
1333	INTEGER :: ji ,jj, jk ! dummy loop indices
1334	REAL(wp) :: z1_ht
1335	REAL(wp), DIMENSION(jpi,jpj) :: zhf
1336	!!----------------------------------------------------------------------
1337	!
1338	SELECT CASE( nvor_scheme )
1339	CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)
1340	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
1341	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
1342	DO jj = 1, jpj-1
1343	DO ji = 1, jpi-1
1344	zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &
1345	& ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp
1346	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1347	END DO
1348	END DO
1349	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
1350	DO jj = 1, jpj-1
1351	DO ji = 1, jpi-1
1352	zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) &
1353	& + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) &
1354	& / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) &
1355	& + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) )
1356	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
1357	END DO
1358	END DO
1359	END SELECT
1360	!
1361	DO jj = 2, jpj-1
1362	DO ji = 2, jpi-1
1363	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
1364	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
1365	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
1366	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
1367	END DO
1368	END DO
1369	CALL lbc_lnk_multi( 'dynspg_ts', ftne, 'F', 1._wp, ftnw, 'F', 1._wp &
1370	& , ftse, 'F', 1._wp, ftsw, 'F', 1._wp, khlcom = nn_hls+nn_hlts )
1371	!
1372	CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme)
1373	DO jj = 2, jpj
1374	DO ji = 2, jpi
1375	z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )
1376	ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht
1377	ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht
1378	ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
1379	ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht
1380	END DO
1381	END DO
1382	CALL lbc_lnk_multi( 'dynspg_ts', ftne, 'F', 1._wp, ftnw, 'F', 1._wp &
1383	& , ftse, 'F', 1._wp, ftsw, 'F', 1._wp, khlcom = nn_hls+nn_hlts )
1384	!
1385	CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !
1386	!
1387	zwz(:,:) = 0._wp
1388	zhf(:,:) = 0._wp
1389
1390	!!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed
1391	!!gm A priori a better value should be something like :
1392	!!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)
1393	!!gm divided by the sum of the corresponding mask
1394	!!gm
1395	!!
1396	IF( .NOT.ln_sco ) THEN
1397
1398	!!gm agree the JC comment : this should be done in a much clear way
1399
1400	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
1401	! Set it to zero for the time being
1402	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
1403	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
1404	! ENDIF
1405	! zhf(:,:) = gdepw_0(:,:,jk+1)
1406	!
1407	ELSE
1408	!
1409	!zhf(:,:) = hbatf(:,:)
1410	DO jj = 1, jpjm1
1411	DO ji = 1, jpim1
1412	zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &
1413	& + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &
1414	& / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &
1415	& + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )
1416	END DO
1417	END DO
1418	ENDIF
1419	!
1420	DO jj = 1, jpjm1
1421	zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
1422	END DO
1423	!
1424	DO jk = 1, jpkm1
1425	DO jj = 1, jpjm1
1426	zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
1427	END DO
1428	END DO
1429	! JC: TBC. hf should be greater than 0
1430	DO jj = 2, jpjm1
1431	DO ji = 2, jpim1
1432	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj)
1433	END DO
1434	END DO
1435	zwz(2:jpim1,2:jpjm1) = ff_f(2:jpim1,2:jpjm1) * zwz(2:jpim1,2:jpjm1)
1436	CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp, khlcom = nn_hls+nn_hlts )
1437	END SELECT
1438
1439	END SUBROUTINE dyn_cor_2d_init
1440
1441
1442
1443	SUBROUTINE dyn_cor_2d( phgtt, phgtu, phgtv, pun, pvn, phU, phV, kdbi , kdei , kdbj , kdej , pu_trd, pv_trd &
1444	& , kdbi2, kdei2, kdbj2, kdej2 )
1445	!!---------------------------------------------------------------------
1446	!! * ROUTINE dyn_cor_2d *
1447	!!
1448	!! ** Purpose : Compute u and v coriolis trends
1449	!!
1450	!! kdXX2 are useful in the initialisation where some arrays are not over the whole domain
1451	!! and some are
1452	!!----------------------------------------------------------------------
1453	REAL(wp), DIMENSION(kdbi :kdei ,kdbj :kdej ), INTENT(in ) :: phgtt, phgtu, phgtv, pun, pvn ! height, speed
1454	INTEGER , INTENT(in ) :: kdbi , kdei , kdbj , kdej ! arrays size
1455	REAL(wp), DIMENSION(kdbi2:kdei2,kdbj2:kdej2), INTENT(in ) :: phU, phV ! flux
1456	REAL(wp), DIMENSION(kdbi2:kdei2,kdbj2:kdej2), INTENT( out) :: pu_trd, pv_trd
1457	INTEGER , INTENT(in ) :: kdbi2, kdei2, kdbj2, kdej2 ! arrays size
1458	INTEGER :: ji, jj ! dummy loop indices
1459	REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! local integer
1460	!!----------------------------------------------------------------------
1461	SELECT CASE( nvor_scheme )
1462	CASE( np_ENT ) ! enstrophy conserving scheme (f-point)
1463	DO jj = kdbj+1, kdej-1
1464	DO ji = kdbi+1, kdei-1
1465	z1_hu = ssumask_xtd(ji,jj) / ( phgtu(ji,jj) + 1._wp - ssumask_xtd(ji,jj) )
1466	z1_hv = ssvmask_xtd(ji,jj) / ( phgtv(ji,jj) + 1._wp - ssvmask_xtd(ji,jj) )
1467	pu_trd(ji,jj) = + r1_4 * r1_e1e2u_xtd(ji,jj) * z1_hu &
1468	& * ( e1e2t_xtd(ji+1,jj)phgtt(ji+1,jj)ff_t_xtd(ji+1,jj) * ( pvn(ji+1,jj) + pvn(ji+1,jj-1) ) &
1469	& + e1e2t_xtd(ji ,jj)phgtt(ji ,jj)ff_t_xtd(ji ,jj) * ( pvn(ji ,jj) + pvn(ji ,jj-1) ) )
1470	!
1471	pv_trd(ji,jj) = - r1_4 * r1_e1e2v_xtd(ji,jj) * z1_hv &
1472	& * ( e1e2t_xtd(ji,jj+1)phgtt(ji,jj+1)ff_t_xtd(ji,jj+1) * ( pun(ji,jj+1) + pun(ji-1,jj+1) ) &
1473	& + e1e2t_xtd(ji,jj )phgtt(ji,jj )ff_t_xtd(ji,jj ) * ( pun(ji,jj ) + pun(ji-1,jj ) ) )
1474	END DO
1475	END DO
1476	!
1477	CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX
1478	DO jj = kdbj+1, kdej-1
1479	DO ji = kdbi+1, kdei-1
1480	zy1 = ( phV(ji,jj-1) + phV(ji+1,jj-1) ) * r1_e1u_xtd(ji,jj)
1481	zy2 = ( phV(ji,jj ) + phV(ji+1,jj ) ) * r1_e1u_xtd(ji,jj)
1482	zx1 = ( phU(ji-1,jj) + phU(ji-1,jj+1) ) * r1_e2v_xtd(ji,jj)
1483	zx2 = ( phU(ji ,jj) + phU(ji ,jj+1) ) * r1_e2v_xtd(ji,jj)
1484	! energy conserving formulation for planetary vorticity term
1485	pu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1486	pv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
1487	END DO
1488	END DO
1489	!
1490	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
1491	DO jj = kdbj+1, kdej-1
1492	DO ji = kdbi+1, kdei-1
1493	zy1 = r1_8 * ( phV(ji ,jj-1) + phV(ji+1,jj-1) &
1494	& + phV(ji ,jj ) + phV(ji+1,jj ) ) * r1_e1u_xtd(ji,jj)
1495	zx1 = - r1_8 * ( phU(ji-1,jj ) + phU(ji-1,jj+1) &
1496	& + phU(ji ,jj ) + phU(ji ,jj+1) ) * r1_e2v_xtd(ji,jj)
1497	pu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
1498	pv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
1499	END DO
1500	END DO
1501	!
1502	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
1503	DO jj = kdbj+1, kdej-1
1504	DO ji = kdbi+1, kdei-1
1505	pu_trd(ji,jj) = + r1_12 * r1_e1u_xtd(ji,jj) * ( ftne(ji,jj ) * phV(ji ,jj ) &
1506	& + ftnw(ji+1,jj) * phV(ji+1,jj ) &
1507	& + ftse(ji,jj ) * phV(ji ,jj-1) &
1508	& + ftsw(ji+1,jj) * phV(ji+1,jj-1) )
1509	pv_trd(ji,jj) = - r1_12 * r1_e2v_xtd(ji,jj) * ( ftsw(ji,jj+1) * phU(ji-1,jj+1) &
1510	& + ftse(ji,jj+1) * phU(ji ,jj+1) &
1511	& + ftnw(ji,jj ) * phU(ji-1,jj ) &
1512	& + ftne(ji,jj ) * phU(ji ,jj ) )
1513	END DO
1514	END DO
1515	!
1516	END SELECT
1517	!
1518	END SUBROUTINE dyn_cor_2d
1519
1520
1521	SUBROUTINE wad_tmsk( pssh, ptmsk )
1522	!!----------------------------------------------------------------------
1523	!! * ROUTINE wad_lmt *
1524	!!
1525	!! ** Purpose : set wetting & drying mask at tracer points
1526	!! for the current barotropic sub-step
1527	!!
1528	!! ** Method : ???
1529	!!
1530	!! ** Action : ptmsk : wetting & drying t-mask
1531	!!----------------------------------------------------------------------
1532	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh !
1533	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk !
1534	!
1535	INTEGER :: ji, jj ! dummy loop indices
1536	!!----------------------------------------------------------------------
1537	!
1538	IF( ln_wd_dl_rmp ) THEN
1539	DO jj = 1, jpj
1540	DO ji = 1, jpi
1541	IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN
1542	! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN
1543	ptmsk(ji,jj) = 1._wp
1544	ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN
1545	ptmsk(ji,jj) = TANH( 50._wp( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )r_rn_wdmin1) )
1546	ELSE
1547	ptmsk(ji,jj) = 0._wp
1548	ENDIF
1549	END DO
1550	END DO
1551	ELSE
1552	DO jj = 1, jpj
1553	DO ji = 1, jpi
1554	IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp
1555	ELSE ; ptmsk(ji,jj) = 0._wp
1556	ENDIF
1557	END DO
1558	END DO
1559	ENDIF
1560	!
1561	END SUBROUTINE wad_tmsk
1562
1563
1564	SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
1565	!!----------------------------------------------------------------------
1566	!! * ROUTINE wad_lmt *
1567	!!
1568	!! ** Purpose : set wetting & drying mask at tracer points
1569	!! for the current barotropic sub-step
1570	!!
1571	!! ** Method : ???
1572	!!
1573	!! ** Action : ptmsk : wetting & drying t-mask
1574	!!----------------------------------------------------------------------
1575	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask
1576	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports
1577	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask
1578	!
1579	INTEGER :: ji, jj ! dummy loop indices
1580	!!----------------------------------------------------------------------
1581	!
1582	DO jj = 1, jpj
1583	DO ji = 1, jpim1 ! not jpi-column
1584	IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj)
1585	ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj)
1586	ENDIF
1587	phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
1588	pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
1589	END DO
1590	END DO
1591	!
1592	DO jj = 1, jpjm1 ! not jpj-row
1593	DO ji = 1, jpi
1594	IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj )
1595	ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1)
1596	ENDIF
1597	phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj)
1598	pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
1599	END DO
1600	END DO
1601	!
1602	END SUBROUTINE wad_Umsk
1603
1604
1605	SUBROUTINE wad_spg( sshn, zcpx, zcpy )
1606	!!---------------------------------------------------------------------
1607	!! * ROUTINE wad_sp *
1608	!!
1609	!! ** Purpose :
1610	!!----------------------------------------------------------------------
1611	INTEGER :: ji ,jj ! dummy loop indices
1612	LOGICAL :: ll_tmp1, ll_tmp2
1613	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: sshn
1614	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
1615	!!----------------------------------------------------------------------
1616	DO jj = 2, jpjm1
1617	DO ji = 2, jpim1
1618	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > &
1619	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
1620	& MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) &
1621	& > rn_wdmin1 + rn_wdmin2
1622	ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( &
1623	& MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > &
1624	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1625	IF(ll_tmp1) THEN
1626	zcpx(ji,jj) = 1.0_wp
1627	ELSEIF(ll_tmp2) THEN
1628	! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here
1629	zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
1630	& / (sshn(ji+1,jj) - sshn(ji ,jj)) )
1631	zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1632	ELSE
1633	zcpx(ji,jj) = 0._wp
1634	ENDIF
1635	!
1636	ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > &
1637	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
1638	& MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) &
1639	& > rn_wdmin1 + rn_wdmin2
1640	ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( &
1641	& MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > &
1642	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1643
1644	IF(ll_tmp1) THEN
1645	zcpy(ji,jj) = 1.0_wp
1646	ELSE IF(ll_tmp2) THEN
1647	! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here
1648	zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
1649	& / (sshn(ji,jj+1) - sshn(ji,jj )) )
1650	zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) )
1651	ELSE
1652	zcpy(ji,jj) = 0._wp
1653	ENDIF
1654	END DO
1655	END DO
1656
1657	END SUBROUTINE wad_spg
1658
1659
1660
1661	SUBROUTINE drg_init( kdbi, kdei, kdbj, kdej, pu_frc, pv_frc, pCdU_u, pCdU_v )
1662	!!----------------------------------------------------------------------
1663	!! * ROUTINE drg_init *
1664	!!
1665	!! ** Purpose : - add the baroclinic top/bottom drag contribution to
1666	!! the baroclinic part of the barotropic RHS
1667	!! - compute the barotropic drag coefficients
1668	!!
1669	!! ** Method : computation done over the INNER domain only
1670	!!----------------------------------------------------------------------
1671	INTEGER , INTENT(in ) :: kdbi, kdei, kdbj, kdej ! arrays size
1672	REAL(wp), DIMENSION(kdbi:kdei,kdbj:kdej), INTENT(inout) :: pu_frc, pv_frc ! baroclinic part of the barotropic RHS
1673	REAL(wp), DIMENSION(kdbi:kdei,kdbj:kdej), INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients
1674	!
1675	INTEGER :: ji, jj ! dummy loop indices
1676	INTEGER :: ikbu, ikbv, iktu, iktv
1677	REAL(wp) :: zztmp
1678	REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i
1679	!!----------------------------------------------------------------------
1680	!
1681	! !== Set the barotropic drag coef. ==!
1682	!
1683	IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)
1684
1685	DO jj = 2, jpjm1
1686	DO ji = 2, jpim1 ! INNER domain
1687	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
1688	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
1689	END DO
1690	END DO
1691	ELSE ! bottom friction only
1692	DO jj = 2, jpjm1
1693	DO ji = 2, jpim1 ! INNER domain
1694	pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
1695	pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
1696	END DO
1697	END DO
1698	ENDIF
1699	!
1700	! !== BOTTOM stress contribution from baroclinic velocities ==!
1701	!
1702	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
1703
1704	DO jj = 2, jpjm1
1705	DO ji = 2, jpim1 ! INNER domain
1706	ikbu = mbku(ji,jj)
1707	ikbv = mbkv(ji,jj)
1708	zu_i(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj)
1709	zv_i(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
1710	END DO
1711	END DO
1712	ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
1713
1714	DO jj = 2, jpjm1
1715	DO ji = 2, jpim1 ! INNER domain
1716	ikbu = mbku(ji,jj)
1717	ikbv = mbkv(ji,jj)
1718	zu_i(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj)
1719	zv_i(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
1720	END DO
1721	END DO
1722	ENDIF
1723	!
1724	IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please !
1725	zztmp = -1._wp / rdtbt
1726	DO jj = 2, jpjm1
1727	DO ji = 2, jpim1 ! INNER domain
1728	pu_frc(ji,jj) = pu_frc(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &
1729	& r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp )
1730	pv_frc(ji,jj) = pv_frc(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &
1731	& r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
1732	END DO
1733	END DO
1734	ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
1735
1736	DO jj = 2, jpjm1
1737	DO ji = 2, jpim1 ! INNER domain
1738	pu_frc(ji,jj) = pu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) zu_i(ji,jj)
1739	pv_frc(ji,jj) = pv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) zv_i(ji,jj)
1740	END DO
1741	END DO
1742	END IF
1743	!
1744	! !== TOP stress contribution from baroclinic velocities ==! (no W/D case)
1745	!
1746	IF( ln_isfcav ) THEN
1747	!
1748	IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
1749
1750	DO jj = 2, jpjm1
1751	DO ji = 2, jpim1 ! INNER domain
1752	iktu = miku(ji,jj)
1753	iktv = mikv(ji,jj)
1754	zu_i(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj)
1755	zv_i(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
1756	END DO
1757	END DO
1758	ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
1759
1760	DO jj = 2, jpjm1
1761	DO ji = 2, jpim1 ! INNER domain
1762	iktu = miku(ji,jj)
1763	iktv = mikv(ji,jj)
1764	zu_i(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj)
1765	zv_i(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
1766	END DO
1767	END DO
1768	ENDIF
1769	!
1770	! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
1771
1772	DO jj = 2, jpjm1
1773	DO ji = 2, jpim1 ! INNER domain
1774	pu_frc(ji,jj) = pu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) zu_i(ji,jj)
1775	pv_frc(ji,jj) = pv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) zv_i(ji,jj)
1776	END DO
1777	END DO
1778	!
1779	ENDIF
1780	!
1781	END SUBROUTINE drg_init
1782
1783
1784	SUBROUTINE ts_bck_interp( km, ld_init, & ! <== in
1785	& za0, za1, za2, za3 ) ! ==> out
1786	!!----------------------------------------------------------------------
1787	INTEGER ,INTENT(in ) :: km ! index of sub time step
1788	LOGICAL ,INTENT(in ) :: ld_init !
1789	REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient
1790	!
1791	REAL(wp) :: zepsilon, zgamma ! - -
1792	!!----------------------------------------------------------------------
1793	! ! set Half-step back interpolation coefficient
1794	IF ( km==1 .AND. ld_init ) THEN !* Forward-backward
1795	za0 = 1._wp
1796	za1 = 0._wp
1797	za2 = 0._wp
1798	za3 = 0._wp
1799	ELSEIF( km==2 .AND. ld_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
1800	za0 = 1.0833333333333_wp ! za0 = 1-gam-eps
1801	za1 =-0.1666666666666_wp ! za1 = gam
1802	za2 = 0.0833333333333_wp ! za2 = eps
1803	za3 = 0._wp
1804	ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
1805	IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion
1806	za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps
1807	za1 = 0.285_wp ! za1 = 1/2 - 2gam - 3eps
1808	za2 = 0.088_wp ! za2 = gam
1809	za3 = 0.013_wp ! za3 = eps
1810	ELSE ! no time diffusion
1811	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
1812	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
1813	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
1814	za1 = 1._wp - za0 - zgamma - zepsilon
1815	za2 = zgamma
1816	za3 = zepsilon
1817	ENDIF
1818	ENDIF
1819	END SUBROUTINE ts_bck_interp
1820
1821
1822	!!======================================================================
1823	END MODULE dynspg_ts

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

Download in other formats: