Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90 @ 7019

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

Download in other formats: