Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 62.5 KB

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

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

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

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7910

Download in other formats: