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

Context Navigation

source: branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 6603

Last change on this file since 6603 was 6603, checked in by hliu, 8 years ago
wet/dry for dynspg_ts.F90
Property svn:keywords set to `Id`
File size: 64.4 KB

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

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

Download in other formats: