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

Context Navigation

source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5211

Last change on this file since 5211 was 5211, checked in by hliu, 9 years ago
add bottom friction limiter in dynspg_ts.F90
Property svn:keywords set to `Id`
File size: 61.0 KB

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

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

Download in other formats: