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

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

dynspg_ts.F90 in NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynspg_ts.F90 @ 10876

Last change on this file since 10876 was 10009, checked in by gm, 6 years ago
#1911 (ENHANCE-04): RK3 branch - step II.1 time-level dimension on ssh
Property svn:keywords set to `Id`
File size: 76.8 KB

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

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

Download in other formats: