Context Navigation

dynspg_ts.F90 @ 5914

Last change on this file since 5914 was 5914, checked in by rfurner, 8 years ago
pulled over bug fix from trunk, see nemo ticket #1633
Property svn:keywords set to `Id`
File size: 63.4 KB

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

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/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5914

Download in other formats: