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_r8624_ENHANCE4_FREESURFACE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2017/dev_r8624_ENHANCE4_FREESURFACE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 8842

Last change on this file since 8842 was 8842, checked in by jchanut, 7 years ago
Enhance4-freesurface. step 3: diffusive barotropic time stepping - #1959
Property svn:keywords set to `Id`
File size: 65.0 KB

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

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

Download in other formats: