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/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5014

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

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

Download in other formats: