Context Navigation

dynspg_ts.F90 @ 6389

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 6389

Download in other formats: