Context Navigation

dynspg_ts.F90 @ 11917

Last change on this file since 11917 was 11917, checked in by davestorkey, 4 years ago

UKMO/NEMO_4.0_momentum_trends :

Refactor so that ZDF trend always includes all bottom friction contributions. The budget can now be closed by including the ZDF trend without any BFR trends.

Some tidying of trddyn.F90.

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

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

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

Context Navigation

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

Download in other formats: