Context Navigation

dynspg_ts.F90 @ 10272

Last change on this file since 10272 was 10272, checked in by deazer, 6 years ago

Bug fix for wd case with implicit bottom friction.
(Currently only done for ln_dynadv_vec)
Trivial, just move lines so that the are done in all advection cases

see Ticket:

https://forge.ipsl.jussieu.fr/nemo/ticket/2149

Property svn:keywords set to Id

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

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

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

Context Navigation

source: NEMO/trunk/src/OCE/DYN/dynspg_ts.F90 @ 10272

Download in other formats: