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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynspg_ts.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg_ts.F90 @ 10919

Last change on this file since 10919 was 10919, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps :

Introduce ssh, uu_b, vv_b variables with time dimension (and pointers).
Convert dynspg modules (but no change to subcycling in dynspg_ts.F90).

Property svn:keywords set to Id

File size: 78.4 KB

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

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

Download in other formats: