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

Context Navigation

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5842

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

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

Download in other formats: