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_momentum_trends/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynspg_ts.F90 @ 11613

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

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

Download in other formats: