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/UKMO/NEMO_4.0.4_generic_obs/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/DYN/dynspg_ts.F90 @ 15248

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

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

Download in other formats: