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/2016/dev_merge_2016/NEMOGCM/CONFIG/WAD_TEST_CASES/MY_SRC – NEMO

Context Navigation

source: branches/2016/dev_merge_2016/NEMOGCM/CONFIG/WAD_TEST_CASES/MY_SRC/dynspg_ts.F90 @ 7485

Last change on this file since 7485 was 7467, checked in by acc, 8 years ago
Changes to WAD_TEST_CASES/MY_SRC files and namelist_cfg to successfully run with the merge branch. Some changes will be moved into the main source once verified
File size: 64.3 KB

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

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

Download in other formats: