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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9528

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

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

Download in other formats: