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 branches/2017/dev_METO_2017/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2017/dev_METO_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8985

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

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

Download in other formats: