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_MLF/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynspg_ts.F90 @ 9863

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

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

Download in other formats: