Context Navigation

dynspg_ts.F90 @ 9568

Last change on this file since 9568 was 9568, checked in by davestorkey, 6 years ago
Update branch to be relative to rev 9565 of dev_merge_2017.
File size: 76.3 KB

Line
1	MODULE dynspg_ts
2
3	!! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out !
4
5	!!======================================================================
6	!! * MODULE dynspg_ts *
7	!! Ocean dynamics: surface pressure gradient trend, split-explicit scheme
8	!!======================================================================
9	!! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code
10	!! - ! 2005-11 (V. Garnier, G. Madec) optimization
11	!! - ! 2006-08 (S. Masson) distributed restart using iom
12	!! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines
13	!! - ! 2008-01 (R. Benshila) change averaging method
14	!! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
15	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
16	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
17	!! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping
18	!! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility
19	!! 3.7 ! 2015-11 (J. Chanut) free surface simplification
20	!! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence
21	!! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90)
22	!!---------------------------------------------------------------------
23
24	!!----------------------------------------------------------------------
25	!! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme
26	!! dyn_spg_ts_init: initialisation of the time-splitting scheme
27	!! ts_wgt : set time-splitting weights for temporal averaging (or not)
28	!! ts_rst : read/write time-splitting fields in restart file
29	!!----------------------------------------------------------------------
30	USE oce ! ocean dynamics and tracers
31	USE dom_oce ! ocean space and time domain
32	USE sbc_oce ! surface boundary condition: ocean
33	USE zdf_oce ! vertical physics: variables
34	USE zdfdrg ! vertical physics: top/bottom drag coef.
35	USE sbcisf ! ice shelf variable (fwfisf)
36	USE sbcapr ! surface boundary condition: atmospheric pressure
37	USE dynadv , ONLY: ln_dynadv_vec
38	USE dynvor ! vortivity scheme indicators
39	USE phycst ! physical constants
40	USE dynvor ! vorticity term
41	USE wet_dry ! wetting/drying flux limter
42	USE bdy_oce ! open boundary
43	USE bdytides ! open boundary condition data
44	USE bdydyn2d ! open boundary conditions on barotropic variables
45	USE sbctide ! tides
46	USE updtide ! tide potential
47	USE sbcwave ! surface wave
48	USE diatmb ! Top,middle,bottom output
49	#if defined key_agrif
50	USE agrif_opa_interp ! agrif
51	USE agrif_oce
52	#endif
53	#if defined key_asminc
54	USE asminc ! Assimilation increment
55	#endif
56	!
57	USE in_out_manager ! I/O manager
58	USE lib_mpp ! distributed memory computing library
59	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
60	USE prtctl ! Print control
61	USE iom ! IOM library
62	USE restart ! only for lrst_oce
63	USE diatmb ! Top,middle,bottom output
64
65	IMPLICIT NONE
66	PRIVATE
67
68	PUBLIC dyn_spg_ts ! called by dyn_spg
69	PUBLIC dyn_spg_ts_init ! - - dyn_spg_init
70
71	!! Time filtered arrays at baroclinic time step:
72	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step
73	!
74	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
75	REAL(wp),SAVE :: rdtbt ! Barotropic time step
76	!
77	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields
78	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points
79	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
80	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
81
82	REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios
83	REAL(wp) :: r1_8 = 0.125_wp !
84	REAL(wp) :: r1_4 = 0.25_wp !
85	REAL(wp) :: r1_2 = 0.5_wp !
86
87	!! * Substitutions
88	# include "vectopt_loop_substitute.h90"
89	!!----------------------------------------------------------------------
90	!! NEMO/OPA 4.0 , NEMO Consortium (2017)
91	!! $Id$
92	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
93	!!----------------------------------------------------------------------
94	CONTAINS
95
96	INTEGER FUNCTION dyn_spg_ts_alloc()
97	!!----------------------------------------------------------------------
98	!! * routine dyn_spg_ts_alloc *
99	!!----------------------------------------------------------------------
100	INTEGER :: ierr(3)
101	!!----------------------------------------------------------------------
102	ierr(:) = 0
103	!
104	ALLOCATE( wgtbtp1(3nn_baro), wgtbtp2(3nn_baro), zwz(jpi,jpj), STAT=ierr(1) )
105	!
106	IF( ln_dynvor_een .OR. ln_dynvor_eeT ) &
107	& ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &
108	& ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) )
109	!
110	ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) )
111	!
112	dyn_spg_ts_alloc = MAXVAL( ierr(:) )
113	!
114	IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc )
115	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays')
116	!
117	END FUNCTION dyn_spg_ts_alloc
118
119
120	SUBROUTINE dyn_spg_ts( kt )
121	!!----------------------------------------------------------------------
122	!!
123	!! ** Purpose : - Compute the now trend due to the explicit time stepping
124	!! of the quasi-linear barotropic system, and add it to the
125	!! general momentum trend.
126	!!
127	!! ** Method : - split-explicit schem (time splitting) :
128	!! Barotropic variables are advanced from internal time steps
129	!! "n" to "n+1" if ln_bt_fw=T
130	!! or from
131	!! "n-1" to "n+1" if ln_bt_fw=F
132	!! thanks to a generalized forward-backward time stepping (see ref. below).
133	!!
134	!! ** Action :
135	!! -Update the filtered free surface at step "n+1" : ssha
136	!! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
137	!! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv
138	!! These are used to advect tracers and are compliant with discrete
139	!! continuity equation taken at the baroclinic time steps. This
140	!! ensures tracers conservation.
141	!! - (ua, va) momentum trend updated with barotropic component.
142	!!
143	!! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.
144	!!---------------------------------------------------------------------
145	INTEGER, INTENT(in) :: kt ! ocean time-step index
146	!
147	INTEGER :: ji, jj, jk, jn ! dummy loop indices
148	LOGICAL :: ll_fw_start ! =T : forward integration
149	LOGICAL :: ll_init ! =T : special startup of 2d equations
150	LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D
151	INTEGER :: ikbu, iktu, noffset ! local integers
152	INTEGER :: ikbv, iktv ! - -
153	REAL(wp) :: r1_2dt_b, z2dt_bf ! local scalars
154	REAL(wp) :: zx1, zx2, zu_spg, zhura, z1_hu ! - -
155	REAL(wp) :: zy1, zy2, zv_spg, zhvra, z1_hv ! - -
156	REAL(wp) :: za0, za1, za2, za3 ! - -
157	REAL(wp) :: zmdi, zztmp , z1_ht ! - -
158	REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf
159	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc
160	REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv
161	REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e
162	REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e
163	REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points
164	!
165	REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True.
166
167	INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point
168
169	REAL(wp) :: zepsilon, zgamma ! - -
170	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef.
171	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
172	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask
173	!!----------------------------------------------------------------------
174	!
175	IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) )
176	! !* Allocate temporary arrays
177	IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj))
178	!
179	zmdi=1.e+20 ! missing data indicator for masking
180	!
181	zwdramp = r_rn_wdmin1 ! simplest ramp
182	! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
183	! ! reciprocal of baroclinic time step
184	IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt
185	ELSE ; z2dt_bf = 2.0_wp * rdt
186	ENDIF
187	r1_2dt_b = 1.0_wp / z2dt_bf
188	!
189	ll_init = ln_bt_av ! if no time averaging, then no specific restart
190	ll_fw_start = .FALSE.
191	! ! time offset in steps for bdy data update
192	IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro
193	ELSE ; noffset = 0
194	ENDIF
195	!
196	IF( kt == nit000 ) THEN !* initialisation
197	!
198	IF(lwp) WRITE(numout,*)
199	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
200	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
201	IF(lwp) WRITE(numout,*)
202	!
203	IF( neuler == 0 ) ll_init=.TRUE.
204	!
205	IF( ln_bt_fw .OR. neuler == 0 ) THEN
206	ll_fw_start =.TRUE.
207	noffset = 0
208	ELSE
209	ll_fw_start =.FALSE.
210	ENDIF
211	!
212	! Set averaging weights and cycle length:
213	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
214	!
215	ENDIF
216	!
217	IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)
218	DO jj = 2, jpjm1
219	DO ji = fs_2, fs_jpim1 ! vector opt.
220	zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
221	zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
222	END DO
223	END DO
224	ELSE ! bottom friction only
225	DO jj = 2, jpjm1
226	DO ji = fs_2, fs_jpim1 ! vector opt.
227	zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
228	zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
229	END DO
230	END DO
231	ENDIF
232	!
233	! Set arrays to remove/compute coriolis trend.
234	! Do it once at kt=nit000 if volume is fixed, else at each long time step.
235	! Note that these arrays are also used during barotropic loop. These are however frozen
236	! although they should be updated in the variable volume case. Not a big approximation.
237	! To remove this approximation, copy lines below inside barotropic loop
238	! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
239	!
240	IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN
241	!
242	SELECT CASE( nvor_scheme )
243	CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)
244	SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point
245	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
246	DO jj = 1, jpjm1
247	DO ji = 1, jpim1
248	zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &
249	& ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp
250	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
251	END DO
252	END DO
253	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
254	DO jj = 1, jpjm1
255	DO ji = 1, jpim1
256	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	END DO
804	ENDIF
805	IF((nbondi == 1).OR.(nbondi == 2)) THEN
806	DO jj=1,jpj
807	zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
808	END DO
809	ENDIF
810	IF((nbondj == -1).OR.(nbondj == 2)) THEN
811	DO ji=1,jpi
812	zwy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)
813	END DO
814	ENDIF
815	IF((nbondj == 1).OR.(nbondj == 2)) THEN
816	DO ji=1,jpi
817	zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
818	END DO
819	ENDIF
820	ENDIF
821	#endif
822	IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
823
824	IF ( ln_wd_dl ) THEN
825	!
826	! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells
827	!
828	DO jj = 1, jpjm1
829	DO ji = 1, jpim1
830	IF ( zwx(ji,jj) > 0.0 ) THEN
831	zuwdmask(ji, jj) = ztwdmask(ji ,jj)
832	ELSE
833	zuwdmask(ji, jj) = ztwdmask(ji+1,jj)
834	END IF
835	zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj)
836	un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj)
837
838	IF ( zwy(ji,jj) > 0.0 ) THEN
839	zvwdmask(ji, jj) = ztwdmask(ji, jj )
840	ELSE
841	zvwdmask(ji, jj) = ztwdmask(ji, jj+1)
842	END IF
843	zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj)
844	vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj)
845	END DO
846	END DO
847	!
848	ENDIF
849
850	! Sum over sub-time-steps to compute advective velocities
851	za2 = wgtbtp2(jn)
852	un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
853	vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
854
855	! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True)
856	IF ( ln_wd_dl_bc ) THEN
857	zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:)
858	zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:)
859	END IF
860
861	! Set next sea level:
862	DO jj = 2, jpjm1
863	DO ji = fs_2, fs_jpim1 ! vector opt.
864	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
865	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj)
866	END DO
867	END DO
868	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
869
870	CALL lbc_lnk( ssha_e, 'T', 1._wp )
871
872	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
873	IF( ln_bdy ) CALL bdy_ssh( ssha_e )
874	#if defined key_agrif
875	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
876	#endif
877	!
878	! Sea Surface Height at u-,v-points (vvl case only)
879	IF( .NOT.ln_linssh ) THEN
880	DO jj = 2, jpjm1
881	DO ji = 2, jpim1 ! NO Vector Opt.
882	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
883	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
884	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
885	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
886	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
887	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
888	END DO
889	END DO
890	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
891	ENDIF
892	!
893	! Half-step back interpolation of SSH for surface pressure computation:
894	!----------------------------------------------------------------------
895	IF ((jn==1).AND.ll_init) THEN
896	za0=1._wp ! Forward-backward
897	za1=0._wp
898	za2=0._wp
899	za3=0._wp
900	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
901	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
902	za1=-0.1666666666666_wp ! za1 = gam
903	za2= 0.0833333333333_wp ! za2 = eps
904	za3= 0._wp
905	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
906	IF (rn_bt_alpha==0._wp) THEN
907	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
908	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
909	za2=0.088_wp ! za2 = gam
910	za3=0.013_wp ! za3 = eps
911	ELSE
912	zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
913	zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
914	za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
915	za1 = 1._wp - za0 - zgamma - zepsilon
916	za2 = zgamma
917	za3 = zepsilon
918	ENDIF
919	ENDIF
920	!
921	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
922	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
923
924	IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters
925	DO jj = 2, jpjm1
926	DO ji = 2, jpim1
927	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
928	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
929	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) &
930	& > rn_wdmin1 + rn_wdmin2
931	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( &
932	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > &
933	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
934
935	IF(ll_tmp1) THEN
936	zcpx(ji,jj) = 1.0_wp
937	ELSE IF(ll_tmp2) THEN
938	! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here
939	zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
940	& / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) )
941	ELSE
942	zcpx(ji,jj) = 0._wp
943	ENDIF
944	!
945	ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
946	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
947	& MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) &
948	& > rn_wdmin1 + rn_wdmin2
949	ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( &
950	& MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > &
951	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
952
953	IF(ll_tmp1) THEN
954	zcpy(ji,jj) = 1.0_wp
955	ELSEIF(ll_tmp2) THEN
956	! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here
957	zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
958	& / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) )
959	ELSE
960	zcpy(ji,jj) = 0._wp
961	ENDIF
962	END DO
963	END DO
964	ENDIF
965	!
966	! Compute associated depths at U and V points:
967	IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form
968	!
969	DO jj = 2, jpjm1
970	DO ji = 2, jpim1
971	zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) &
972	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) &
973	& + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
974	zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) &
975	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) &
976	& + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
977	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
978	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
979	END DO
980	END DO
981	!
982	ENDIF
983	!
984	! Add Coriolis trend:
985	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
986	! at each time step. We however keep them constant here for optimization.
987	! Recall that zwx and zwy arrays hold fluxes at this stage:
988	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
989	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
990	!
991	SELECT CASE( nvor_scheme )
992	CASE( np_ENT ) ! energy conserving scheme (t-point)
993	DO jj = 2, jpjm1
994	DO ji = 2, jpim1 ! vector opt.
995
996	z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) )
997	z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) )
998
999	zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu &
1000	& * ( 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) ) &
1001	& + e1e2t(ji ,jj)zhtp2_e(ji ,jj)ff_t(ji ,jj) * ( va_e(ji ,jj) + va_e(ji ,jj-1) ) )
1002	!
1003	zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv &
1004	& * ( 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) ) &
1005	& + e1e2t(ji,jj )zhtp2_e(ji,jj )ff_t(ji,jj ) * ( ua_e(ji,jj ) + ua_e(ji-1,jj ) ) )
1006	END DO
1007	END DO
1008	!
1009	CASE( np_ENE, np_MIX ) ! energy conserving scheme (f-point)
1010	DO jj = 2, jpjm1
1011	DO ji = fs_2, fs_jpim1 ! vector opt.
1012	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
1013	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
1014	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
1015	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
1016	zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
1017	zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
1018	END DO
1019	END DO
1020	!
1021	CASE( np_ENS ) ! enstrophy conserving scheme (f-point)
1022	DO jj = 2, jpjm1
1023	DO ji = fs_2, fs_jpim1 ! vector opt.
1024	zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
1025	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
1026	zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
1027	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
1028	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
1029	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
1030	END DO
1031	END DO
1032	!
1033	CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)
1034	DO jj = 2, jpjm1
1035	DO ji = fs_2, fs_jpim1 ! vector opt.
1036	zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
1037	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
1038	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
1039	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
1040	zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
1041	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
1042	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
1043	& + ftne(ji,jj ) * zwx(ji ,jj ) )
1044	END DO
1045	END DO
1046	!
1047	END SELECT
1048	!
1049	! Add tidal astronomical forcing if defined
1050	IF ( ln_tide .AND. ln_tide_pot ) THEN
1051	DO jj = 2, jpjm1
1052	DO ji = fs_2, fs_jpim1 ! vector opt.
1053	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
1054	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
1055	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
1056	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
1057	END DO
1058	END DO
1059	ENDIF
1060	!
1061	! Add bottom stresses:
1062	!jth do implicitly instead
1063	IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs
1064	DO jj = 2, jpjm1
1065	DO ji = fs_2, fs_jpim1 ! vector opt.
1066	zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)
1067	zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)
1068	END DO
1069	END DO
1070	ENDIF
1071	!
1072	! Surface pressure trend:
1073	IF( ln_wd_il ) THEN
1074	DO jj = 2, jpjm1
1075	DO ji = 2, jpim1
1076	! Add surface pressure gradient
1077	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1078	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1079	zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj)
1080	zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj)
1081	END DO
1082	END DO
1083	ELSE
1084	DO jj = 2, jpjm1
1085	DO ji = fs_2, fs_jpim1 ! vector opt.
1086	! Add surface pressure gradient
1087	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
1088	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
1089	zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg
1090	zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg
1091	END DO
1092	END DO
1093	END IF
1094
1095	!
1096	! Set next velocities:
1097	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
1098	DO jj = 2, jpjm1
1099	DO ji = fs_2, fs_jpim1 ! vector opt.
1100	ua_e(ji,jj) = ( un_e(ji,jj) &
1101	& + rdtbt * ( zwx(ji,jj) &
1102	& + zu_trd(ji,jj) &
1103	& + zu_frc(ji,jj) ) &
1104	& ) * ssumask(ji,jj)
1105
1106	va_e(ji,jj) = ( vn_e(ji,jj) &
1107	& + rdtbt * ( zwy(ji,jj) &
1108	& + zv_trd(ji,jj) &
1109	& + zv_frc(ji,jj) ) &
1110	& ) * ssvmask(ji,jj)
1111
1112	!jth implicit bottom friction:
1113	IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs
1114	ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj))
1115	va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj))
1116	ENDIF
1117
1118	END DO
1119	END DO
1120	!
1121	ELSE !* Flux form
1122	DO jj = 2, jpjm1
1123	DO ji = fs_2, fs_jpim1 ! vector opt.
1124
1125	zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1126	zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1127
1128	zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
1129	zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
1130
1131	ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) &
1132	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
1133	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
1134	& + hu_n(ji,jj) * zu_frc(ji,jj) ) &
1135	& ) * zhura
1136
1137	va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) &
1138	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
1139	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
1140	& + hv_n(ji,jj) * zv_frc(ji,jj) ) &
1141	& ) * zhvra
1142	END DO
1143	END DO
1144	ENDIF
1145
1146
1147	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
1148	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1149	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
1150	hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
1151	hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
1152	!
1153	ENDIF
1154	! !* domain lateral boundary
1155	CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
1156	!
1157	! ! open boundaries
1158	IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
1159	#if defined key_agrif
1160	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
1161	#endif
1162	! !* Swap
1163	! ! ----
1164	ubb_e (:,:) = ub_e (:,:)
1165	ub_e (:,:) = un_e (:,:)
1166	un_e (:,:) = ua_e (:,:)
1167	!
1168	vbb_e (:,:) = vb_e (:,:)
1169	vb_e (:,:) = vn_e (:,:)
1170	vn_e (:,:) = va_e (:,:)
1171	!
1172	sshbb_e(:,:) = sshb_e(:,:)
1173	sshb_e (:,:) = sshn_e(:,:)
1174	sshn_e (:,:) = ssha_e(:,:)
1175
1176	! !* Sum over whole bt loop
1177	! ! ----------------------
1178	za1 = wgtbtp1(jn)
1179	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
1180	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
1181	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
1182	ELSE ! Sum transports
1183	IF ( .NOT.ln_wd_dl ) THEN
1184	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
1185	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
1186	ELSE
1187	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
1188	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
1189	END IF
1190	ENDIF
1191	! ! Sum sea level
1192	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1193
1194	! ! ==================== !
1195	END DO ! end loop !
1196	! ! ==================== !
1197	! -----------------------------------------------------------------------------
1198	! Phase 3. update the general trend with the barotropic trend
1199	! -----------------------------------------------------------------------------
1200	!
1201	! Set advection velocity correction:
1202	IF (ln_bt_fw) THEN
1203	zwx(:,:) = un_adv(:,:)
1204	zwy(:,:) = vn_adv(:,:)
1205	IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
1206	un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
1207	vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
1208	!
1209	! Update corrective fluxes for next time step:
1210	un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
1211	vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
1212	ELSE
1213	un_bf(:,:) = 0._wp
1214	vn_bf(:,:) = 0._wp
1215	END IF
1216	! Save integrated transport for next computation
1217	ub2_b(:,:) = zwx(:,:)
1218	vb2_b(:,:) = zwy(:,:)
1219	ENDIF
1220
1221
1222	!
1223	! Update barotropic trend:
1224	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1225	DO jk=1,jpkm1
1226	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b
1227	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b
1228	END DO
1229	ELSE
1230	! At this stage, ssha has been corrected: compute new depths at velocity points
1231	DO jj = 1, jpjm1
1232	DO ji = 1, jpim1 ! NO Vector Opt.
1233	zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
1234	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
1235	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1236	zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
1237	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
1238	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1239	END DO
1240	END DO
1241	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1242	!
1243	DO jk=1,jpkm1
1244	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b
1245	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b
1246	END DO
1247	! Save barotropic velocities not transport:
1248	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1249	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1250	ENDIF
1251
1252
1253	! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)
1254	DO jk = 1, jpkm1
1255	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)r1_hu_n(:,:) - un_b(:,:) ) umask(:,:,jk)
1256	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)r1_hv_n(:,:) - vn_b(:,:) ) vmask(:,:,jk)
1257	END DO
1258
1259	IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
1260	DO jk = 1, jpkm1
1261	un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &
1262	& + zuwdav2(:,:)(un(:,:,jk) - un_adv(:,:)r1_hu_n(:,:)) ) * umask(:,:,jk)
1263	vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) &
1264	& + zvwdav2(:,:)(vn(:,:,jk) - vn_adv(:,:)r1_hv_n(:,:)) ) * vmask(:,:,jk)
1265	END DO
1266	END IF
1267
1268
1269	CALL iom_put( "ubar", un_adv(:,:)*r1_hu_n(:,:) ) ! barotropic i-current
1270	CALL iom_put( "vbar", vn_adv(:,:)*r1_hv_n(:,:) ) ! barotropic i-current
1271	!
1272	#if defined key_agrif
1273	! Save time integrated fluxes during child grid integration
1274	! (used to update coarse grid transports at next time step)
1275	!
1276	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1277	IF( Agrif_NbStepint() == 0 ) THEN
1278	ub2_i_b(:,:) = 0._wp
1279	vb2_i_b(:,:) = 0._wp
1280	END IF
1281	!
1282	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1283	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1284	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1285	ENDIF
1286	#endif
1287	! !* write time-spliting arrays in the restart
1288	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
1289	!
1290	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
1291	IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
1292	!
1293	IF( ln_diatmb ) THEN
1294	CALL iom_put( "baro_u" , un_bssumask(:,:)+zmdi(1.-ssumask(:,:) ) ) ! Barotropic U Velocity
1295	CALL iom_put( "baro_v" , vn_bssvmask(:,:)+zmdi(1.-ssvmask(:,:) ) ) ! Barotropic V Velocity
1296	ENDIF
1297	!
1298	END SUBROUTINE dyn_spg_ts
1299
1300
1301	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1302	!!---------------------------------------------------------------------
1303	!! * ROUTINE ts_wgt *
1304	!!
1305	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1306	!!----------------------------------------------------------------------
1307	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
1308	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
1309	INTEGER, INTENT(inout) :: jpit ! cycle length
1310	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
1311	zwgt2 ! Secondary weights
1312
1313	INTEGER :: jic, jn, ji ! temporary integers
1314	REAL(wp) :: za1, za2
1315	!!----------------------------------------------------------------------
1316
1317	zwgt1(:) = 0._wp
1318	zwgt2(:) = 0._wp
1319
1320	! Set time index when averaged value is requested
1321	IF (ll_fw) THEN
1322	jic = nn_baro
1323	ELSE
1324	jic = 2 * nn_baro
1325	ENDIF
1326
1327	! Set primary weights:
1328	IF (ll_av) THEN
1329	! Define simple boxcar window for primary weights
1330	! (width = nn_baro, centered around jic)
1331	SELECT CASE ( nn_bt_flt )
1332	CASE( 0 ) ! No averaging
1333	zwgt1(jic) = 1._wp
1334	jpit = jic
1335
1336	CASE( 1 ) ! Boxcar, width = nn_baro
1337	DO jn = 1, 3*nn_baro
1338	za1 = ABS(float(jn-jic))/float(nn_baro)
1339	IF (za1 < 0.5_wp) THEN
1340	zwgt1(jn) = 1._wp
1341	jpit = jn
1342	ENDIF
1343	ENDDO
1344
1345	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
1346	DO jn = 1, 3*nn_baro
1347	za1 = ABS(float(jn-jic))/float(nn_baro)
1348	IF (za1 < 1._wp) THEN
1349	zwgt1(jn) = 1._wp
1350	jpit = jn
1351	ENDIF
1352	ENDDO
1353	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1354	END SELECT
1355
1356	ELSE ! No time averaging
1357	zwgt1(jic) = 1._wp
1358	jpit = jic
1359	ENDIF
1360
1361	! Set secondary weights
1362	DO jn = 1, jpit
1363	DO ji = jn, jpit
1364	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1365	END DO
1366	END DO
1367
1368	! Normalize weigths:
1369	za1 = 1._wp / SUM(zwgt1(1:jpit))
1370	za2 = 1._wp / SUM(zwgt2(1:jpit))
1371	DO jn = 1, jpit
1372	zwgt1(jn) = zwgt1(jn) * za1
1373	zwgt2(jn) = zwgt2(jn) * za2
1374	END DO
1375	!
1376	END SUBROUTINE ts_wgt
1377
1378
1379	SUBROUTINE ts_rst( kt, cdrw )
1380	!!---------------------------------------------------------------------
1381	!! * ROUTINE ts_rst *
1382	!!
1383	!! ** Purpose : Read or write time-splitting arrays in restart file
1384	!!----------------------------------------------------------------------
1385	INTEGER , INTENT(in) :: kt ! ocean time-step
1386	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1387	!!----------------------------------------------------------------------
1388	!
1389	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
1390	! ! ---------------
1391	IF( ln_rstart .AND. ln_bt_fw ) THEN !* Read the restart file
1392	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios )
1393	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios )
1394	CALL iom_get( numror, jpdom_autoglo, 'un_bf' , un_bf (:,:), ldxios = lrxios )
1395	CALL iom_get( numror, jpdom_autoglo, 'vn_bf' , vn_bf (:,:), ldxios = lrxios )
1396	IF( .NOT.ln_bt_av ) THEN
1397	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:), ldxios = lrxios )
1398	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:), ldxios = lrxios )
1399	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:), ldxios = lrxios )
1400	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:), ldxios = lrxios )
1401	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:), ldxios = lrxios )
1402	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:), ldxios = lrxios )
1403	ENDIF
1404	#if defined key_agrif
1405	! Read time integrated fluxes
1406	IF ( .NOT.Agrif_Root() ) THEN
1407	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lrxios )
1408	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lrxios )
1409	ENDIF
1410	#endif
1411	ELSE !* Start from rest
1412	IF(lwp) WRITE(numout,*)
1413	IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
1414	ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
1415	un_adv(:,:) = 0._wp ; vn_adv(:,:) = 0._wp ! used in the 1st interpol of agrif
1416	un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
1417	#if defined key_agrif
1418	IF ( .NOT.Agrif_Root() ) THEN
1419	ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
1420	ENDIF
1421	#endif
1422	ENDIF
1423	!
1424	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
1425	! ! -------------------
1426	IF(lwp) WRITE(numout,*) '---- ts_rst ----'
1427	IF( lwxios ) CALL iom_swap( cwxios_context )
1428	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios )
1429	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios )
1430	CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios )
1431	CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios )
1432	!
1433	IF (.NOT.ln_bt_av) THEN
1434	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:), ldxios = lwxios )
1435	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:), ldxios = lwxios )
1436	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:), ldxios = lwxios )
1437	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:), ldxios = lwxios )
1438	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:), ldxios = lwxios )
1439	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:), ldxios = lwxios )
1440	ENDIF
1441	#if defined key_agrif
1442	! Save time integrated fluxes
1443	IF ( .NOT.Agrif_Root() ) THEN
1444	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lwxios )
1445	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lwxios )
1446	ENDIF
1447	#endif
1448	IF( lwxios ) CALL iom_swap( cxios_context )
1449	ENDIF
1450	!
1451	END SUBROUTINE ts_rst
1452
1453
1454	SUBROUTINE dyn_spg_ts_init
1455	!!---------------------------------------------------------------------
1456	!! * ROUTINE dyn_spg_ts_init *
1457	!!
1458	!! ** Purpose : Set time splitting options
1459	!!----------------------------------------------------------------------
1460	INTEGER :: ji ,jj ! dummy loop indices
1461	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
1462	REAL(wp), DIMENSION(jpi,jpj) :: zcu
1463	!!----------------------------------------------------------------------
1464	!
1465	! Max courant number for ext. grav. waves
1466	!
1467	DO jj = 1, jpj
1468	DO ji =1, jpi
1469	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1470	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1471	zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
1472	END DO
1473	END DO
1474	!
1475	zcmax = MAXVAL( zcu(:,:) )
1476	IF( lk_mpp ) CALL mpp_max( zcmax )
1477
1478	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1479	IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1480
1481	rdtbt = rdt / REAL( nn_baro , wp )
1482	zcmax = zcmax * rdtbt
1483	! Print results
1484	IF(lwp) WRITE(numout,*)
1485	IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
1486	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
1487	IF( ln_bt_auto ) THEN
1488	IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_baro '
1489	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1490	ELSE
1491	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist nn_baro = ', nn_baro
1492	ENDIF
1493
1494	IF(ln_bt_av) THEN
1495	IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_baro time steps is on '
1496	ELSE
1497	IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
1498	ENDIF
1499	!
1500	!
1501	IF(ln_bt_fw) THEN
1502	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1503	ELSE
1504	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1505	ENDIF
1506	!
1507	#if defined key_agrif
1508	! Restrict the use of Agrif to the forward case only
1509	!!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1510	#endif
1511	!
1512	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1513	SELECT CASE ( nn_bt_flt )
1514	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1515	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1516	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1517	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
1518	END SELECT
1519	!
1520	IF(lwp) WRITE(numout,*) ' '
1521	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1522	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1523	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1524	!
1525	IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
1526	IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
1527	CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
1528	ENDIF
1529	!
1530	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1531	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1532	ENDIF
1533	IF( zcmax>0.9_wp ) THEN
1534	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1535	ENDIF
1536	!
1537	! ! Allocate time-splitting arrays
1538	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
1539	!
1540	! ! read restart when needed
1541	CALL ts_rst( nit000, 'READ' )
1542	!
1543	IF( lwxios ) THEN
1544	! define variables in restart file when writing with XIOS
1545	CALL iom_set_rstw_var_active('ub2_b')
1546	CALL iom_set_rstw_var_active('vb2_b')
1547	CALL iom_set_rstw_var_active('un_bf')
1548	CALL iom_set_rstw_var_active('vn_bf')
1549	!
1550	IF (.NOT.ln_bt_av) THEN
1551	CALL iom_set_rstw_var_active('sshbb_e')
1552	CALL iom_set_rstw_var_active('ubb_e')
1553	CALL iom_set_rstw_var_active('vbb_e')
1554	CALL iom_set_rstw_var_active('sshb_e')
1555	CALL iom_set_rstw_var_active('ub_e')
1556	CALL iom_set_rstw_var_active('vb_e')
1557	ENDIF
1558	#if defined key_agrif
1559	! Save time integrated fluxes
1560	IF ( .NOT.Agrif_Root() ) THEN
1561	CALL iom_set_rstw_var_active('ub2_i_b')
1562	CALL iom_set_rstw_var_active('vb2_i_b')
1563	ENDIF
1564	#endif
1565	ENDIF
1566	!
1567	END SUBROUTINE dyn_spg_ts_init
1568
1569	!!======================================================================
1570	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: branches/UKMO/dev_merge_2017_restart_datestamp_GO6_mixing/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9568

Download in other formats: