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/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/OCE/DYN/dynspg_ts.F90 @ 11219

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

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

Download in other formats: