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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynspg_ts.F90 in branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5066

Last change on this file since 5066 was 5066, checked in by rfurner, 9 years ago
added current state of wetting and drying code to test...note it does not work
Property svn:keywords set to `Id`
File size: 61.0 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=-2*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	IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient
374	IF(ln_wd) THEN ! Calculating and applying W/D gravity filters
375	DO jj = 2, jpjm1
376	DO ji = 2, jpim1
377	ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))
378	ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + &
379	& rn_wdmin1 + rn_wdmin2
380	IF(ll_tmp1) THEN
381	zcpx(ji,jj) = 1.0_wp
382	ELSE IF(ll_tmp2) THEN
383	! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
384	zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&
385	& (sshn(ji+1,jj) - sshn(ji,jj)))
386	ELSE
387	zcpx(ji,jj) = 0._wp
388	END IF
389
390	ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))
391	ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + &
392	& rn_wdmin1 + rn_wdmin2
393	IF(ll_tmp1) THEN
394	zcpy(ji,jj) = 1.0_wp
395	ELSE IF(ll_tmp2) THEN
396	! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
397	zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) + bathy(ji,jj)) /&
398	& (sshn(ji,jj+1) - sshn(ji,jj)))
399	ELSE
400	zcpy(ji,jj) = 0._wp
401	END IF
402	END DO
403	END DO
404	CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp )
405	ENDIF
406
407	IF(ln_wd) THEN ! Calculating and applying W/D gravity filters
408	DO jj = 2, jpjm1
409	DO ji = 2, jpim1
410	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / &
411	& e1u(ji,jj) * zcpx(ji,jj)
412	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / &
413	& e2v(ji,jj) * zcpy(ji,jj)
414	END DO
415	END DO
416	ELSE
417	DO jj = 2, jpjm1
418	DO ji = fs_2, fs_jpim1 ! vector opt.
419	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj)
420	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj)
421	END DO
422	END DO
423	END IF
424	ENDIF
425
426	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
427	DO ji = fs_2, fs_jpim1
428	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)
429	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
430	END DO
431	END DO
432	!
433	! ! Add bottom stress contribution from baroclinic velocities:
434	IF (ln_bt_fw) THEN
435	DO jj = 2, jpjm1
436	DO ji = fs_2, fs_jpim1 ! vector opt.
437	ikbu = mbku(ji,jj)
438	ikbv = mbkv(ji,jj)
439	zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
440	zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
441	END DO
442	END DO
443	ELSE
444	DO jj = 2, jpjm1
445	DO ji = fs_2, fs_jpim1 ! vector opt.
446	ikbu = mbku(ji,jj)
447	ikbv = mbkv(ji,jj)
448	zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
449	zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
450	END DO
451	END DO
452	ENDIF
453	!
454	! Note that the "unclipped" bottom friction parameter is used even with explicit drag
455	zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
456	zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
457	!
458	IF (ln_bt_fw) THEN ! Add wind forcing
459	zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
460	zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)
461	ELSE
462	zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)
463	zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
464	ENDIF
465	!
466	IF ( ln_apr_dyn ) THEN ! Add atm pressure forcing
467	IF (ln_bt_fw) THEN
468	DO jj = 2, jpjm1
469	DO ji = fs_2, fs_jpim1 ! vector opt.
470	zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
471	zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)
472	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
473	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
474	END DO
475	END DO
476	ELSE
477	DO jj = 2, jpjm1
478	DO ji = fs_2, fs_jpim1 ! vector opt.
479	zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
480	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj)
481	zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
482	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)
483	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
484	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
485	END DO
486	END DO
487	ENDIF
488	ENDIF
489	! !* Right-Hand-Side of the barotropic ssh equation
490	! ! -----------------------------------------------
491	! ! Surface net water flux and rivers
492	IF (ln_bt_fw) THEN
493	zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) )
494	ELSE
495	zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:))
496	ENDIF
497	#if defined key_asminc
498	! ! Include the IAU weighted SSH increment
499	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
500	zssh_frc(:,:) = zssh_frc(:,:) + ssh_iau(:,:)
501	ENDIF
502	#endif
503	! !* Fill boundary data arrays with AGRIF
504	! ! -------------------------------------
505	#if defined key_agrif
506	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
507	#endif
508	!
509	! -----------------------------------------------------------------------
510	! Phase 2 : Integration of the barotropic equations
511	! -----------------------------------------------------------------------
512	!
513	! ! ==================== !
514	! ! Initialisations !
515	! ! ==================== !
516	! Initialize barotropic variables:
517
518	IF(ln_wd) THEN !preserve the positivity of water depth
519	!ssh[b,n,a] should have already been processed for this
520	sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:))
521	sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:))
522	ENDIF
523
524	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
525	sshn_e(:,:) = sshn (:,:)
526	zun_e (:,:) = un_b (:,:)
527	zvn_e (:,:) = vn_b (:,:)
528	!
529	hu_e (:,:) = hu (:,:)
530	hv_e (:,:) = hv (:,:)
531	hur_e (:,:) = hur (:,:)
532	hvr_e (:,:) = hvr (:,:)
533	ELSE ! CENTRED integration: start from BEFORE fields
534	sshn_e(:,:) = sshb (:,:)
535	zun_e (:,:) = ub_b (:,:)
536	zvn_e (:,:) = vb_b (:,:)
537	!
538	hu_e (:,:) = hu_b (:,:)
539	hv_e (:,:) = hv_b (:,:)
540	hur_e (:,:) = hur_b(:,:)
541	hvr_e (:,:) = hvr_b(:,:)
542	ENDIF
543	!
544	!
545	!
546	! Initialize sums:
547	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
548	va_b (:,:) = 0._wp
549	ssha (:,:) = 0._wp ! Sum for after averaged sea level
550	zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop
551	zv_sum(:,:) = 0._wp
552	! ! ==================== !
553	DO jn = 1, icycle ! sub-time-step loop !
554	! ! ==================== !
555	! !* Update the forcing (BDY and tides)
556	! ! ------------------
557	! Update only tidal forcing at open boundaries
558	#if defined key_tide
559	IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
560	IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset )
561	#endif
562	!
563	! Set extrapolation coefficients for predictor step:
564	IF ((jn<3).AND.ll_init) THEN ! Forward
565	za1 = 1._wp
566	za2 = 0._wp
567	za3 = 0._wp
568	ELSE ! AB3-AM4 Coefficients: bet=0.281105
569	za1 = 1.781105_wp ! za1 = 3/2 + bet
570	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
571	za3 = 0.281105_wp ! za3 = bet
572	ENDIF
573
574	! Extrapolate barotropic velocities at step jit+0.5:
575	ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
576	va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
577
578	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
579	! ! ------------------
580	! Extrapolate Sea Level at step jit+0.5:
581	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
582	!
583	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
584	DO ji = 2, fs_jpim1 ! Vector opt.
585	zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
586	& * ( e12t(ji ,jj) * zsshp2_e(ji ,jj) &
587	& + e12t(ji+1,jj) * zsshp2_e(ji+1,jj) )
588	zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
589	& * ( e12t(ji,jj ) * zsshp2_e(ji,jj ) &
590	& + e12t(ji,jj+1) * zsshp2_e(ji,jj+1) )
591	END DO
592	END DO
593	CALL lbc_lnk( zwx, 'U', 1._wp ) ; CALL lbc_lnk( zwy, 'V', 1._wp )
594	!
595	zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
596	zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
597	IF(ln_wd) THEN
598	zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1)
599	zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1)
600	END IF
601	ELSE
602	zhup2_e (:,:) = hu(:,:)
603	zhvp2_e (:,:) = hv(:,:)
604	ENDIF
605	! !* after ssh
606	! ! -----------
607	! One should enforce volume conservation at open boundaries here
608	! considering fluxes below:
609	!
610	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
611	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
612	!
613	#if defined key_agrif
614	! Set fluxes during predictor step to ensure
615	! volume conservation
616	IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN
617	IF((nbondi == -1).OR.(nbondi == 2)) THEN
618	DO jj=1,jpj
619	zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
620	END DO
621	ENDIF
622	IF((nbondi == 1).OR.(nbondi == 2)) THEN
623	DO jj=1,jpj
624	zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
625	END DO
626	ENDIF
627	IF((nbondj == -1).OR.(nbondj == 2)) THEN
628	DO ji=1,jpi
629	zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
630	END DO
631	ENDIF
632	IF((nbondj == 1).OR.(nbondj == 2)) THEN
633	DO ji=1,jpi
634	zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
635	END DO
636	ENDIF
637	ENDIF
638	#endif
639	IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
640	!
641	! Sum over sub-time-steps to compute advective velocities
642	za2 = wgtbtp2(jn)
643	zu_sum (:,:) = zu_sum (:,:) + za2 * zwx (:,:) / e2u (:,:)
644	zv_sum (:,:) = zv_sum (:,:) + za2 * zwy (:,:) / e1v (:,:)
645	!
646	! Set next sea level:
647	DO jj = 2, jpjm1
648	DO ji = fs_2, fs_jpim1 ! vector opt.
649	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
650	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e12t(ji,jj)
651	END DO
652	END DO
653	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)
654	IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))
655	CALL lbc_lnk( ssha_e, 'T', 1._wp )
656
657	#if defined key_bdy
658	! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
659	IF (lk_bdy) CALL bdy_ssh( ssha_e )
660	#endif
661	#if defined key_agrif
662	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
663	#endif
664	!
665	! Sea Surface Height at u-,v-points (vvl case only)
666	IF ( lk_vvl ) THEN
667	DO jj = 2, jpjm1
668	DO ji = 2, jpim1 ! NO Vector Opt.
669	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
670	& * ( e12t(ji ,jj ) * ssha_e(ji ,jj ) &
671	& + e12t(ji+1,jj ) * ssha_e(ji+1,jj ) )
672	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
673	& * ( e12t(ji ,jj ) * ssha_e(ji ,jj ) &
674	& + e12t(ji ,jj+1) * ssha_e(ji ,jj+1) )
675	END DO
676	END DO
677	CALL lbc_lnk( zsshu_a, 'U', 1._wp ) ; CALL lbc_lnk( zsshv_a, 'V', 1._wp )
678	ENDIF
679	!
680	! Half-step back interpolation of SSH for surface pressure computation:
681	!----------------------------------------------------------------------
682	IF ((jn==1).AND.ll_init) THEN
683	za0=1._wp ! Forward-backward
684	za1=0._wp
685	za2=0._wp
686	za3=0._wp
687	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
688	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
689	za1=-0.1666666666666_wp ! za1 = gam
690	za2= 0.0833333333333_wp ! za2 = eps
691	za3= 0._wp
692	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
693	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
694	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
695	za2=0.088_wp ! za2 = gam
696	za3=0.013_wp ! za3 = eps
697	ENDIF
698
699	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
700	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
701
702	!
703	! Compute associated depths at U and V points:
704	IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN
705	!
706	DO jj = 2, jpjm1
707	DO ji = 2, jpim1
708	zx1 = z1_2 * umask(ji ,jj,1) * r1_e12u(ji ,jj) &
709	& * ( e12t(ji ,jj ) * zsshp2_e(ji ,jj) &
710	& + e12t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
711	zy1 = z1_2 * vmask(ji ,jj,1) * r1_e12v(ji ,jj ) &
712	& * ( e12t(ji ,jj ) * zsshp2_e(ji ,jj ) &
713	& + e12t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
714	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
715	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
716	END DO
717	END DO
718
719	IF(ln_wd) THEN
720	zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 )
721	zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 )
722	END IF
723	!shall we call lbc_lnk for zhu(v)st_e() here?
724
725	ENDIF
726	!
727	! Add Coriolis trend:
728	! zwz array below or triads normally depend on sea level with key_vvl and should be updated
729	! at each time step. We however keep them constant here for optimization.
730	! Recall that zwx and zwy arrays hold fluxes at this stage:
731	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
732	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
733	!
734	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
735	DO jj = 2, jpjm1
736	DO ji = fs_2, fs_jpim1 ! vector opt.
737	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
738	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
739	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
740	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
741	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
742	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
743	END DO
744	END DO
745	!
746	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
747	DO jj = 2, jpjm1
748	DO ji = fs_2, fs_jpim1 ! vector opt.
749	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
750	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
751	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
752	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
753	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
754	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
755	END DO
756	END DO
757	!
758	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!
759	DO jj = 2, jpjm1
760	DO ji = fs_2, fs_jpim1 ! vector opt.
761	zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
762	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
763	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
764	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
765	zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
766	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
767	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
768	& + ftne(ji,jj ) * zwx(ji ,jj ) )
769	END DO
770	END DO
771	!
772	ENDIF
773	!
774	! Add tidal astronomical forcing if defined
775	IF ( lk_tide.AND.ln_tide_pot ) THEN
776	DO jj = 2, jpjm1
777	DO ji = fs_2, fs_jpim1 ! vector opt.
778	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
779	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
780	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
781	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
782	END DO
783	END DO
784	ENDIF
785	!
786	! Add bottom stresses:
787	zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
788	zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
789	!
790	IF(ln_wd) THEN ! Calculating and applying W/D gravity filters
791	DO jj = 2, jpjm1
792	DO ji = 2, jpim1
793	ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))
794	ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + &
795	& rn_wdmin1 + rn_wdmin2
796	IF(ll_tmp1) THEN
797	zcpx(ji,jj) = 1.0_wp
798	ELSE IF(ll_tmp2) THEN
799	! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here
800	zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /&
801	& (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)))
802	ELSE
803	zcpx(ji,jj) = 0._wp
804	END IF
805
806	ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))
807	ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + &
808	& rn_wdmin1 + rn_wdmin2
809	IF(ll_tmp1) THEN
810	zcpy(ji,jj) = 1.0_wp
811	ELSE IF(ll_tmp2) THEN
812	! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here
813	zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /&
814	& (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)))
815	ELSE
816	zcpy(ji,jj) = 0._wp
817	END IF
818	END DO
819	END DO
820	CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp )
821	ENDIF
822
823	IF(ln_wd) THEN
824	DO jj = 2, jpjm1
825	DO ji = 2, jpim1
826	! Add surface pressure gradient
827	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
828	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
829	zwx(ji,jj) = zu_spg * zcpx(ji,jj)
830	zwy(ji,jj) = zv_spg * zcpy(ji,jj)
831	END DO
832	END DO
833	ELSE
834	DO jj = 2, jpjm1
835	DO ji = fs_2, fs_jpim1 ! vector opt.
836	! Add surface pressure gradient
837	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
838	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
839	zwx(ji,jj) = zu_spg
840	zwy(ji,jj) = zv_spg
841	END DO
842	END DO
843	END IF
844	!
845	! Set next velocities:
846	IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN ! Vector form
847	DO jj = 2, jpjm1
848	DO ji = fs_2, fs_jpim1 ! vector opt.
849	ua_e(ji,jj) = ( zun_e(ji,jj) &
850	& + rdtbt * ( zwx(ji,jj) &
851	& + zu_trd(ji,jj) &
852	& + zu_frc(ji,jj) ) &
853	& ) * umask(ji,jj,1)
854
855	va_e(ji,jj) = ( zvn_e(ji,jj) &
856	& + rdtbt * ( zwy(ji,jj) &
857	& + zv_trd(ji,jj) &
858	& + zv_frc(ji,jj) ) &
859	& ) * vmask(ji,jj,1)
860	END DO
861	END DO
862
863	ELSE ! Flux form
864	DO jj = 2, jpjm1
865	DO ji = fs_2, fs_jpim1 ! vector opt.
866
867	IF(ln_wd) THEN
868	zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1)
869	zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1)
870	ELSE
871	zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
872	zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
873	END IF
874
875	zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1))
876	zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1))
877
878	ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) &
879	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
880	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
881	& + hu(ji,jj) * zu_frc(ji,jj) ) &
882	& ) * zhura
883
884	va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &
885	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
886	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
887	& + hv(ji,jj) * zv_frc(ji,jj) ) &
888	& ) * zhvra
889	END DO
890	END DO
891	ENDIF
892	!
893	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
894	! ! ----------------------------------------------
895	IF(ln_wd) THEN
896	hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1)
897	hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1)
898	ELSE
899	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
900	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
901	END IF
902
903	hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
904	hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
905	!
906	ENDIF
907	! !* domain lateral boundary
908	! ! -----------------------
909	!
910	CALL lbc_lnk( ua_e , 'U', -1._wp ) ! local domain boundaries
911	CALL lbc_lnk( va_e , 'V', -1._wp )
912
913	#if defined key_bdy
914	! open boundaries
915	IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )
916	#endif
917	#if defined key_agrif
918	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
919	#endif
920	! !* Swap
921	! ! ----
922	ubb_e (:,:) = ub_e (:,:)
923	ub_e (:,:) = zun_e (:,:)
924	zun_e (:,:) = ua_e (:,:)
925	!
926	vbb_e (:,:) = vb_e (:,:)
927	vb_e (:,:) = zvn_e (:,:)
928	zvn_e (:,:) = va_e (:,:)
929	!
930	sshbb_e(:,:) = sshb_e(:,:)
931	sshb_e (:,:) = sshn_e(:,:)
932	sshn_e (:,:) = ssha_e(:,:)
933
934	! !* Sum over whole bt loop
935	! ! ----------------------
936	za1 = wgtbtp1(jn)
937	IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities
938	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
939	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
940	ELSE ! Sum transports
941	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
942	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
943	ENDIF
944	! ! Sum sea level
945	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
946	! ! ==================== !
947	END DO ! end loop !
948	! ! ==================== !
949	! -----------------------------------------------------------------------------
950	! Phase 3. update the general trend with the barotropic trend
951	! -----------------------------------------------------------------------------
952	!
953	! At this stage ssha holds a time averaged value
954	! ! Sea Surface Height at u-,v- and f-points
955	IF( lk_vvl ) THEN ! (required only in key_vvl case)
956	DO jj = 1, jpjm1
957	DO ji = 1, jpim1 ! NO Vector Opt.
958	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
959	& * ( e12t(ji ,jj) * ssha(ji ,jj) &
960	& + e12t(ji+1,jj) * ssha(ji+1,jj) )
961	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
962	& * ( e12t(ji,jj ) * ssha(ji,jj ) &
963	& + e12t(ji,jj+1) * ssha(ji,jj+1) )
964	END DO
965	END DO
966	CALL lbc_lnk( zsshu_a, 'U', 1._wp ) ; CALL lbc_lnk( zsshv_a, 'V', 1._wp ) ! Boundary conditions
967	ENDIF
968	!
969	! Set advection velocity correction:
970	IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN
971	un_adv(:,:) = zu_sum(:,:)*hur(:,:)
972	vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
973	ELSE
974	un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
975	vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
976	END IF
977
978	IF (ln_bt_fw) THEN ! Save integrated transport for next computation
979	ub2_b(:,:) = zu_sum(:,:)
980	vb2_b(:,:) = zv_sum(:,:)
981	ENDIF
982	!
983	! Update barotropic trend:
984	IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
985	DO jk=1,jpkm1
986	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
987	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
988	END DO
989	ELSE
990	DO jk=1,jpkm1
991	ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
992	va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
993	END DO
994	! Save barotropic velocities not transport:
995	ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
996	va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
997	ENDIF
998	!
999	DO jk = 1, jpkm1
1000	! Correct velocities:
1001	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
1002	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
1003	!
1004	END DO
1005	!
1006	#if defined key_agrif
1007	! Save time integrated fluxes during child grid integration
1008	! (used to update coarse grid transports)
1009	! Useless with 2nd order momentum schemes
1010	!
1011	IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN
1012	IF ( Agrif_NbStepint().EQ.0 ) THEN
1013	ub2_i_b(:,:) = 0.e0
1014	vb2_i_b(:,:) = 0.e0
1015	END IF
1016	!
1017	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1018	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1019	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1020	ENDIF
1021	!
1022	!
1023	#endif
1024	!
1025	! !* write time-spliting arrays in the restart
1026	IF(lrst_oce .AND.ln_bt_fw) CALL ts_rst( kt, 'WRITE' )
1027	!
1028	CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv )
1029	CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
1030	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
1031	CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
1032	CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a )
1033	CALL wrk_dealloc( jpi, jpj, zhf )
1034	IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy )
1035	!
1036	IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts')
1037	!
1038	END SUBROUTINE dyn_spg_ts
1039
1040	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1041	!!---------------------------------------------------------------------
1042	!! * ROUTINE ts_wgt *
1043	!!
1044	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1045	!!----------------------------------------------------------------------
1046	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
1047	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
1048	INTEGER, INTENT(inout) :: jpit ! cycle length
1049	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
1050	zwgt2 ! Secondary weights
1051
1052	INTEGER :: jic, jn, ji ! temporary integers
1053	REAL(wp) :: za1, za2
1054	!!----------------------------------------------------------------------
1055
1056	zwgt1(:) = 0._wp
1057	zwgt2(:) = 0._wp
1058
1059	! Set time index when averaged value is requested
1060	IF (ll_fw) THEN
1061	jic = nn_baro
1062	ELSE
1063	jic = 2 * nn_baro
1064	ENDIF
1065
1066	! Set primary weights:
1067	IF (ll_av) THEN
1068	! Define simple boxcar window for primary weights
1069	! (width = nn_baro, centered around jic)
1070	SELECT CASE ( nn_bt_flt )
1071	CASE( 0 ) ! No averaging
1072	zwgt1(jic) = 1._wp
1073	jpit = jic
1074
1075	CASE( 1 ) ! Boxcar, width = nn_baro
1076	DO jn = 1, 3*nn_baro
1077	za1 = ABS(float(jn-jic))/float(nn_baro)
1078	IF (za1 < 0.5_wp) THEN
1079	zwgt1(jn) = 1._wp
1080	jpit = jn
1081	ENDIF
1082	ENDDO
1083
1084	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
1085	DO jn = 1, 3*nn_baro
1086	za1 = ABS(float(jn-jic))/float(nn_baro)
1087	IF (za1 < 1._wp) THEN
1088	zwgt1(jn) = 1._wp
1089	jpit = jn
1090	ENDIF
1091	ENDDO
1092	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1093	END SELECT
1094
1095	ELSE ! No time averaging
1096	zwgt1(jic) = 1._wp
1097	jpit = jic
1098	ENDIF
1099
1100	! Set secondary weights
1101	DO jn = 1, jpit
1102	DO ji = jn, jpit
1103	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1104	END DO
1105	END DO
1106
1107	! Normalize weigths:
1108	za1 = 1._wp / SUM(zwgt1(1:jpit))
1109	za2 = 1._wp / SUM(zwgt2(1:jpit))
1110	DO jn = 1, jpit
1111	zwgt1(jn) = zwgt1(jn) * za1
1112	zwgt2(jn) = zwgt2(jn) * za2
1113	END DO
1114	!
1115	END SUBROUTINE ts_wgt
1116
1117	SUBROUTINE ts_rst( kt, cdrw )
1118	!!---------------------------------------------------------------------
1119	!! * ROUTINE ts_rst *
1120	!!
1121	!! ** Purpose : Read or write time-splitting arrays in restart file
1122	!!----------------------------------------------------------------------
1123	INTEGER , INTENT(in) :: kt ! ocean time-step
1124	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1125	!
1126	!!----------------------------------------------------------------------
1127	!
1128	IF( TRIM(cdrw) == 'READ' ) THEN
1129	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) )
1130	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) )
1131	IF( .NOT.ln_bt_av ) THEN
1132	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) )
1133	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) )
1134	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) )
1135	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) )
1136	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) )
1137	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) )
1138	ENDIF
1139	#if defined key_agrif
1140	! Read time integrated fluxes
1141	IF ( .NOT.Agrif_Root() ) THEN
1142	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:) )
1143	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:) )
1144	ENDIF
1145	#endif
1146	!
1147	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1148	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
1149	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
1150	!
1151	IF (.NOT.ln_bt_av) THEN
1152	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
1153	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
1154	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
1155	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
1156	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
1157	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
1158	ENDIF
1159	#if defined key_agrif
1160	! Save time integrated fluxes
1161	IF ( .NOT.Agrif_Root() ) THEN
1162	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
1163	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
1164	ENDIF
1165	#endif
1166	ENDIF
1167	!
1168	END SUBROUTINE ts_rst
1169
1170	SUBROUTINE dyn_spg_ts_init( kt )
1171	!!---------------------------------------------------------------------
1172	!! * ROUTINE dyn_spg_ts_init *
1173	!!
1174	!! ** Purpose : Set time splitting options
1175	!!----------------------------------------------------------------------
1176	INTEGER , INTENT(in) :: kt ! ocean time-step
1177	!
1178	INTEGER :: ji ,jj
1179	INTEGER :: ios ! Local integer output status for namelist read
1180	REAL(wp) :: zxr2, zyr2, zcmax
1181	REAL(wp), POINTER, DIMENSION(:,:) :: zcu
1182	!!
1183	NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
1184	& nn_baro, rn_bt_cmax, nn_bt_flt
1185	!!----------------------------------------------------------------------
1186	!
1187	REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters
1188	READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901)
1189	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp )
1190
1191	REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters
1192	READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 )
1193	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp )
1194	IF(lwm) WRITE ( numond, namsplit )
1195	!
1196	! ! Max courant number for ext. grav. waves
1197	!
1198	CALL wrk_alloc( jpi, jpj, zcu )
1199	!
1200	IF (lk_vvl) THEN
1201	DO jj = 1, jpj
1202	DO ji =1, jpi
1203	zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1204	zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1205	zcu(ji,jj) = sqrt(gravht_0(ji,jj)(zxr2 + zyr2) )
1206	END DO
1207	END DO
1208	ELSE
1209	DO jj = 1, jpj
1210	DO ji =1, jpi
1211	zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1212	zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1213	zcu(ji,jj) = sqrt(gravht(ji,jj)(zxr2 + zyr2) )
1214	END DO
1215	END DO
1216	ENDIF
1217
1218	zcmax = MAXVAL(zcu(:,:))
1219	IF( lk_mpp ) CALL mpp_max( zcmax )
1220
1221	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1222	IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1223
1224	rdtbt = rdt / FLOAT(nn_baro)
1225	zcmax = zcmax * rdtbt
1226	! Print results
1227	IF(lwp) WRITE(numout,*)
1228	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1229	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1230	IF( ln_bt_nn_auto ) THEN
1231	IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro '
1232	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1233	ELSE
1234	IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist '
1235	ENDIF
1236
1237	IF(ln_bt_av) THEN
1238	IF(lwp) WRITE(numout,*) ' ln_bt_av=.true. => Time averaging over nn_baro time steps is on '
1239	ELSE
1240	IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
1241	ENDIF
1242	!
1243	!
1244	IF(ln_bt_fw) THEN
1245	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1246	ELSE
1247	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1248	ENDIF
1249	!
1250	#if defined key_agrif
1251	! Restrict the use of Agrif to the forward case only
1252	IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root())) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1253	#endif
1254	!
1255	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1256	SELECT CASE ( nn_bt_flt )
1257	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1258	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1259	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1260	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1261	END SELECT
1262	!
1263	IF(lwp) WRITE(numout,*) ' '
1264	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1265	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1266	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1267	!
1268	IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
1269	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1270	ENDIF
1271	IF ( zcmax>0.9_wp ) THEN
1272	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1273	ENDIF
1274	!
1275	CALL wrk_dealloc( jpi, jpj, zcu )
1276	!
1277	END SUBROUTINE dyn_spg_ts_init
1278
1279	#else
1280	!!---------------------------------------------------------------------------
1281	!! Default case : Empty module No split explicit free surface
1282	!!---------------------------------------------------------------------------
1283	CONTAINS
1284	INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function
1285	dyn_spg_ts_alloc = 0
1286	END FUNCTION dyn_spg_ts_alloc
1287	SUBROUTINE dyn_spg_ts( kt ) ! Empty routine
1288	INTEGER, INTENT(in) :: kt
1289	WRITE(,) 'dyn_spg_ts: You should not have seen this print! error?', kt
1290	END SUBROUTINE dyn_spg_ts
1291	SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine
1292	INTEGER , INTENT(in) :: kt ! ocean time-step
1293	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1294	WRITE(,) 'ts_rst : You should not have seen this print! error?', kt, cdrw
1295	END SUBROUTINE ts_rst
1296	SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine
1297	INTEGER , INTENT(in) :: kt ! ocean time-step
1298	WRITE(,) 'dyn_spg_ts_init : You should not have seen this print! error?', kt
1299	END SUBROUTINE dyn_spg_ts_init
1300	#endif
1301
1302	!!======================================================================
1303	END MODULE dynspg_ts
1304
1305
1306

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

Download in other formats: