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

Context Navigation

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 6748

Last change on this file since 6748 was 6748, checked in by mocavero, 8 years ago
GYRE hybrid parallelization
Property svn:keywords set to `Id`
File size: 66.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_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	!!----------------------------------------------------------------------
160	!
161	IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts')
162	!
163	! !* Allocate temporary arrays
164	CALL wrk_alloc( jpi,jpj, zsshp2_e, zhdiv )
165	CALL wrk_alloc( jpi,jpj, zu_trd, zv_trd)
166	CALL wrk_alloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc)
167	CALL wrk_alloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e)
168	CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a )
169	CALL wrk_alloc( jpi,jpj, zhf )
170	IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 )
171	!
172	zmdi=1.e+20 ! missing data indicator for masking
173	! !* Local constant initialization
174	z1_12 = 1._wp / 12._wp
175	z1_8 = 0.125_wp
176	z1_4 = 0.25_wp
177	z1_2 = 0.5_wp
178	zraur = 1._wp / rau0
179	! ! reciprocal of baroclinic time step
180	IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt
181	ELSE ; z2dt_bf = 2.0_wp * rdt
182	ENDIF
183	z1_2dt_b = 1.0_wp / z2dt_bf
184	!
185	ll_init = ln_bt_av ! if no time averaging, then no specific restart
186	ll_fw_start = .FALSE.
187	! ! time offset in steps for bdy data update
188	IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro
189	ELSE ; noffset = 0
190	ENDIF
191	!
192	IF( kt == nit000 ) THEN !* initialisation
193	!
194	IF(lwp) WRITE(numout,*)
195	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
196	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
197	IF(lwp) WRITE(numout,*)
198	!
199	IF( neuler == 0 ) ll_init=.TRUE.
200	!
201	IF( ln_bt_fw .OR. neuler == 0 ) THEN
202	ll_fw_start =.TRUE.
203	noffset = 0
204	ELSE
205	ll_fw_start =.FALSE.
206	ENDIF
207	!
208	! Set averaging weights and cycle length:
209	CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
210	!
211	ENDIF
212	!
213	! Set arrays to remove/compute coriolis trend.
214	! Do it once at kt=nit000 if volume is fixed, else at each long time step.
215	! Note that these arrays are also used during barotropic loop. These are however frozen
216	! although they should be updated in the variable volume case. Not a big approximation.
217	! To remove this approximation, copy lines below inside barotropic loop
218	! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
219	!
220	IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN
221	IF( ln_dynvor_een ) THEN !== EEN scheme ==!
222	SELECT CASE( nn_een_e3f ) !* ff/e3 at F-point
223	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
224	!$OMP PARALLEL DO schedule(static) private(jj, ji)
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	!$OMP PARALLEL DO schedule(static) private(jj, ji)
234	DO jj = 1, jpjm1
235	DO ji = 1, jpim1
236	zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &
237	& ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) &
238	& / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + &
239	& tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) )
240	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj)
241	END DO
242	END DO
243	END SELECT
244	CALL lbc_lnk( zwz, 'F', 1._wp )
245	!
246	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
247
248	!$OMP PARALLEL DO schedule(static) private(jj, ji)
249	DO jj = 2, jpj
250	DO ji = 2, jpi
251	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
252	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
253	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
254	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
255	END DO
256	END DO
257	!
258	ELSE !== all other schemes (ENE, ENS, MIX)
259	!$OMP PARALLEL WORKSHARE
260	zwz(:,:) = 0._wp
261	zhf(:,:) = 0._wp
262	!$OMP END PARALLEL WORKSHARE
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	END IF
276
277	!$OMP PARALLEL DO schedule(static) private(jj, ji)
278	DO jj = 1, jpjm1
279	DO ji = 1, jpi
280	zhf(ji,jj) = zhf(ji,jj) * (1._wp- umask(ji,jj,1) * umask(ji,jj+1,1))
281	END DO
282	END DO
283
284	DO jk = 1, jpkm1
285	!$OMP PARALLEL DO schedule(static) private(jj, ji)
286	DO jj = 1, jpjm1
287	DO ji = 1, jpi
288	zhf(ji,jj) = zhf(ji,jj) + e3f_n(ji,jj,jk) * umask(ji,jj,jk) * umask(ji,jj+1,jk)
289	END DO
290	END DO
291	END DO
292	CALL lbc_lnk( zhf, 'F', 1._wp )
293	! JC: TBC. hf should be greater than 0
294	!$OMP PARALLEL DO schedule(static) private(jj, ji)
295	DO jj = 1, jpj
296	DO ji = 1, jpi
297	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
298	END DO
299	END DO
300	!$OMP PARALLEL WORKSHARE
301	zwz(:,:) = ff(:,:) * zwz(:,:)
302	!$OMP END PARALLEL WORKSHARE
303	ENDIF
304	ENDIF
305	!
306	! If forward start at previous time step, and centered integration,
307	! then update averaging weights:
308	IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
309	ll_fw_start=.FALSE.
310	CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
311	ENDIF
312
313	! -----------------------------------------------------------------------------
314	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
315	! -----------------------------------------------------------------------------
316	!
317	!
318	! !* e3*d/dt(Ua) (Vertically integrated)
319	! ! --------------------------------------------------
320	!$OMP PARALLEL WORKSHARE
321	zu_frc(:,:) = 0._wp
322	zv_frc(:,:) = 0._wp
323	!$OMP END PARALLEL WORKSHARE
324	!
325	DO jk = 1, jpkm1
326	!$OMP PARALLEL DO schedule(static) private(jj,ji)
327	DO jj=1,jpj
328	DO ji=1,jpi
329	zu_frc(ji,jj) = zu_frc(ji,jj) + e3u_n(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
330	zv_frc(ji,jj) = zv_frc(ji,jj) + e3v_n(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
331	END DO
332	END DO
333	END DO
334	!
335	!$OMP PARALLEL WORKSHARE
336	zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:)
337	zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:)
338	!$OMP END PARALLEL WORKSHARE
339	!
340	!
341	! !* baroclinic momentum trend (remove the vertical mean trend)
342	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
343	DO jk = 1, jpkm1 ! -----------------------------------------------------------
344	DO jj = 2, jpjm1
345	DO ji = fs_2, fs_jpim1 ! vector opt.
346	ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
347	va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
348	END DO
349	END DO
350	END DO
351	! !* barotropic Coriolis trends (vorticity scheme dependent)
352	! ! --------------------------------------------------------
353	!$OMP PARALLEL WORKSHARE
354	zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes
355	zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)
356	!$OMP END PARALLEL WORKSHARE
357	!
358	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
359	!$OMP PARALLEL DO schedule(static) private(jj,ji,zy1,zy2,zx1,zx2)
360	DO jj = 2, jpjm1
361	DO ji = fs_2, fs_jpim1 ! vector opt.
362	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
363	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
364	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
365	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
366	! energy conserving formulation for planetary vorticity term
367	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
368	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
369	END DO
370	END DO
371	!
372	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
373	!$OMP PARALLEL DO schedule(static) private(jj,ji,zy1,zx1)
374	DO jj = 2, jpjm1
375	DO ji = fs_2, fs_jpim1 ! vector opt.
376	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
377	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
378	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
379	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
380	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
381	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
382	END DO
383	END DO
384	!
385	ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme
386	!$OMP PARALLEL DO schedule(static) private(jj,ji)
387	DO jj = 2, jpjm1
388	DO ji = fs_2, fs_jpim1 ! vector opt.
389	zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
390	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
391	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
392	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
393	zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
394	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
395	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
396	& + ftne(ji,jj ) * zwx(ji ,jj ) )
397	END DO
398	END DO
399	!
400	ENDIF
401	!
402	! !* Right-Hand-Side of the barotropic momentum equation
403	! ! ----------------------------------------------------
404	IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient
405	IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters
406	wduflt1(:,:) = 1.0_wp
407	wdvflt1(:,:) = 1.0_wp
408	!$OMP PARALLEL DO schedule(static) private(jj,ji,ll_tmp1,ll_tmp2)
409	DO jj = 2, jpjm1
410	DO ji = 2, jpim1
411	ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) &
412	& .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) &
413	& > rn_wdmin1 + rn_wdmin2
414	ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) &
415	& + rn_wdmin1 + rn_wdmin2
416	IF(ll_tmp1) THEN
417	zcpx(ji,jj) = 1.0_wp
418	ELSEIF(ll_tmp2) THEN
419	! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen here
420	zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) &
421	& /(sshn(ji+1,jj) - sshn(ji,jj)))
422	ELSE
423	zcpx(ji,jj) = 0._wp
424	wduflt1(ji,jj) = 0.0_wp
425	END IF
426
427	ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) &
428	& .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) &
429	& > rn_wdmin1 + rn_wdmin2
430	ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) &
431	& + rn_wdmin1 + rn_wdmin2
432	IF(ll_tmp1) THEN
433	zcpy(ji,jj) = 1.0_wp
434	ELSEIF(ll_tmp2) THEN
435	! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen here
436	zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) &
437	& /(sshn(ji,jj+1) - sshn(ji,jj)))
438	ELSE
439	zcpy(ji,jj) = 0._wp
440	wdvflt1(ji,jj) = 0.0_wp
441	ENDIF
442
443	END DO
444	END DO
445
446	CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp )
447
448	!$OMP PARALLEL DO schedule(static) private(jj,ji)
449	DO jj = 2, jpjm1
450	DO ji = 2, jpim1
451	zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) &
452	& * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj)
453	zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) &
454	& * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj)
455	END DO
456	END DO
457
458	ELSE
459
460	!$OMP PARALLEL DO schedule(static) private(jj,ji)
461	DO jj = 2, jpjm1
462	DO ji = fs_2, fs_jpim1 ! vector opt.
463	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj)
464	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj)
465	END DO
466	END DO
467	ENDIF
468
469	ENDIF
470
471	!$OMP PARALLEL DO schedule(static) private(jj,ji)
472	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
473	DO ji = fs_2, fs_jpim1
474	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
475	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
476	END DO
477	END DO
478	!
479	! ! Add bottom stress contribution from baroclinic velocities:
480	IF (ln_bt_fw) THEN
481	!$OMP PARALLEL DO schedule(static) private(jj,ji,ikbu,ikbv)
482	DO jj = 2, jpjm1
483	DO ji = fs_2, fs_jpim1 ! vector opt.
484	ikbu = mbku(ji,jj)
485	ikbv = mbkv(ji,jj)
486	zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
487	zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
488	END DO
489	END DO
490	ELSE
491	!$OMP PARALLEL DO schedule(static) private(jj,ji,ikbu,ikbv)
492	DO jj = 2, jpjm1
493	DO ji = fs_2, fs_jpim1 ! vector opt.
494	ikbu = mbku(ji,jj)
495	ikbv = mbkv(ji,jj)
496	zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
497	zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
498	END DO
499	END DO
500	ENDIF
501	!
502	! Note that the "unclipped" bottom friction parameter is used even with explicit drag
503	IF( ln_wd ) THEN
504	zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:)
505	zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:)
506	ELSE
507	zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:)
508	zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:)
509	END IF
510	!
511	! ! Add top stress contribution from baroclinic velocities:
512	IF (ln_bt_fw) THEN
513	!$OMP PARALLEL DO schedule(static) private(jj,ji,ikbu,ikbv)
514	DO jj = 2, jpjm1
515	DO ji = fs_2, fs_jpim1 ! vector opt.
516	iktu = miku(ji,jj)
517	iktv = mikv(ji,jj)
518	zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
519	zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
520	END DO
521	END DO
522	ELSE
523	!$OMP PARALLEL DO schedule(static) private(jj,ji,ikbu,ikbv)
524	DO jj = 2, jpjm1
525	DO ji = fs_2, fs_jpim1 ! vector opt.
526	iktu = miku(ji,jj)
527	iktv = mikv(ji,jj)
528	zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
529	zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
530	END DO
531	END DO
532	ENDIF
533	!
534	! Note that the "unclipped" top friction parameter is used even with explicit drag
535	!$OMP PARALLEL WORKSHARE
536	zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * tfrua(:,:) * zwx(:,:)
537	zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * tfrva(:,:) * zwy(:,:)
538	!$OMP END PARALLEL WORKSHARE
539	!
540	IF (ln_bt_fw) THEN ! Add wind forcing
541	!$OMP PARALLEL WORKSHARE
542	zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:)
543	zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:)
544	!$OMP END PARALLEL WORKSHARE
545	ELSE
546	!$OMP PARALLEL WORKSHARE
547	zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:)
548	zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:)
549	!$OMP END PARALLEL WORKSHARE
550	ENDIF
551	!
552	IF ( ln_apr_dyn ) THEN ! Add atm pressure forcing
553	IF (ln_bt_fw) THEN
554	!$OMP PARALLEL DO schedule(static) private(jj,ji,zu_spg,zv_spg)
555	DO jj = 2, jpjm1
556	DO ji = fs_2, fs_jpim1 ! vector opt.
557	zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
558	zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
559	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
560	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
561	END DO
562	END DO
563	ELSE
564	!$OMP PARALLEL DO schedule(static) private(jj,ji,zu_spg,zv_spg)
565	DO jj = 2, jpjm1
566	DO ji = fs_2, fs_jpim1 ! vector opt.
567	zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
568	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
569	zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
570	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
571	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
572	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
573	END DO
574	END DO
575	ENDIF
576	ENDIF
577	! !* Right-Hand-Side of the barotropic ssh equation
578	! ! -----------------------------------------------
579	! ! Surface net water flux and rivers
580	IF (ln_bt_fw) THEN
581	zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
582	ELSE
583	zssh_frc(:,:) = zraur * z1_2 * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) &
584	& + fwfisf(:,:) + fwfisf_b(:,:) )
585	ENDIF
586	#if defined key_asminc
587	! ! Include the IAU weighted SSH increment
588	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
589	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
590	ENDIF
591	#endif
592	! !* Fill boundary data arrays for AGRIF
593	! ! ------------------------------------
594	#if defined key_agrif
595	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
596	#endif
597	!
598	! -----------------------------------------------------------------------
599	! Phase 2 : Integration of the barotropic equations
600	! -----------------------------------------------------------------------
601	!
602	! ! ==================== !
603	! ! Initialisations !
604	! ! ==================== !
605	! Initialize barotropic variables:
606	IF( ll_init )THEN
607	!$OMP PARALLEL WORKSHARE
608	sshbb_e(:,:) = 0._wp
609	ubb_e (:,:) = 0._wp
610	vbb_e (:,:) = 0._wp
611	sshb_e (:,:) = 0._wp
612	ub_e (:,:) = 0._wp
613	vb_e (:,:) = 0._wp
614	!$OMP END PARALLEL WORKSHARE
615	ENDIF
616
617	IF( ln_wd ) THEN !preserve the positivity of water depth
618	!ssh[b,n,a] should have already been processed for this
619	!$OMP PARALLEL WORKSHARE
620	sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:))
621	sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:))
622	!$OMP END PARALLEL WORKSHARE
623	ENDIF
624	!
625	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
626	!$OMP PARALLEL WORKSHARE
627	sshn_e(:,:) = sshn(:,:)
628	un_e (:,:) = un_b(:,:)
629	vn_e (:,:) = vn_b(:,:)
630	!
631	hu_e (:,:) = hu_n(:,:)
632	hv_e (:,:) = hv_n(:,:)
633	hur_e (:,:) = r1_hu_n(:,:)
634	hvr_e (:,:) = r1_hv_n(:,:)
635	!$OMP END PARALLEL WORKSHARE
636	ELSE ! CENTRED integration: start from BEFORE fields
637	!$OMP PARALLEL WORKSHARE
638	sshn_e(:,:) = sshb(:,:)
639	un_e (:,:) = ub_b(:,:)
640	vn_e (:,:) = vb_b(:,:)
641	!
642	hu_e (:,:) = hu_b(:,:)
643	hv_e (:,:) = hv_b(:,:)
644	hur_e (:,:) = r1_hu_b(:,:)
645	hvr_e (:,:) = r1_hv_b(:,:)
646	!$OMP END PARALLEL WORKSHARE
647	ENDIF
648	!
649	!
650	!
651	! Initialize sums:
652	!$OMP PARALLEL WORKSHARE
653	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
654	va_b (:,:) = 0._wp
655	ssha (:,:) = 0._wp ! Sum for after averaged sea level
656	un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop
657	vn_adv(:,:) = 0._wp
658	!$OMP END PARALLEL WORKSHARE
659	! ! ==================== !
660	DO jn = 1, icycle ! sub-time-step loop !
661	! ! ==================== !
662	! !* Update the forcing (BDY and tides)
663	! ! ------------------
664	! Update only tidal forcing at open boundaries
665	#if defined key_tide
666	IF( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 )
667	IF( ln_tide_pot .AND. lk_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset )
668	#endif
669	!
670	! Set extrapolation coefficients for predictor step:
671	IF ((jn<3).AND.ll_init) THEN ! Forward
672	za1 = 1._wp
673	za2 = 0._wp
674	za3 = 0._wp
675	ELSE ! AB3-AM4 Coefficients: bet=0.281105
676	za1 = 1.781105_wp ! za1 = 3/2 + bet
677	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
678	za3 = 0.281105_wp ! za3 = bet
679	ENDIF
680
681	! Extrapolate barotropic velocities at step jit+0.5:
682	!$OMP PARALLEL WORKSHARE
683	ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
684	va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
685	!$OMP END PARALLEL WORKSHARE
686
687	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
688	! ! ------------------
689	! Extrapolate Sea Level at step jit+0.5:
690	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
691	!
692	!$OMP PARALLEL DO schedule(static) private(jj,ji)
693	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
694	DO ji = 2, fs_jpim1 ! Vector opt.
695	zwx(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
696	& * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) &
697	& + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
698	zwy(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
699	& * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) &
700	& + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
701	END DO
702	END DO
703	CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
704	!
705	!$OMP PARALLEL WORKSHARE
706	zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
707	zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
708	!$OMP END PARALLEL WORKSHARE
709	IF( ln_wd ) THEN
710	zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1)
711	zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1)
712	END IF
713	ELSE
714	!$OMP PARALLEL WORKSHARE
715	zhup2_e (:,:) = hu_n(:,:)
716	zhvp2_e (:,:) = hv_n(:,:)
717	!$OMP END PARALLEL WORKSHARE
718	ENDIF
719	! !* after ssh
720	! ! -----------
721	! One should enforce volume conservation at open boundaries here
722	! considering fluxes below:
723	!
724	!$OMP PARALLEL WORKSHARE
725	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
726	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
727	!$OMP END PARALLEL WORKSHARE
728	!
729	#if defined key_agrif
730	! Set fluxes during predictor step to ensure volume conservation
731	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
732	IF((nbondi == -1).OR.(nbondi == 2)) THEN
733	DO jj=1,jpj
734	zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
735	END DO
736	ENDIF
737	IF((nbondi == 1).OR.(nbondi == 2)) THEN
738	DO jj=1,jpj
739	zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
740	END DO
741	ENDIF
742	IF((nbondj == -1).OR.(nbondj == 2)) THEN
743	DO ji=1,jpi
744	zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
745	END DO
746	ENDIF
747	IF((nbondj == 1).OR.(nbondj == 2)) THEN
748	DO ji=1,jpi
749	zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
750	END DO
751	ENDIF
752	ENDIF
753	#endif
754	IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
755	!
756	! Sum over sub-time-steps to compute advective velocities
757	za2 = wgtbtp2(jn)
758	!$OMP PARALLEL WORKSHARE
759	un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
760	vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
761	!$OMP END PARALLEL WORKSHARE
762	!
763	! Set next sea level:
764	!$OMP PARALLEL DO schedule(static) private(jj,ji)
765	DO jj = 2, jpjm1
766	DO ji = fs_2, fs_jpim1 ! vector opt.
767	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
768	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj)
769	END DO
770	END DO
771	!$OMP PARALLEL WORKSHARE
772	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)
773	!$OMP END PARALLEL WORKSHARE
774	IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))
775	CALL lbc_lnk( ssha_e, 'T', 1._wp )
776
777	#if defined key_bdy
778	! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)
779	IF( lk_bdy ) CALL bdy_ssh( ssha_e )
780	#endif
781	#if defined key_agrif
782	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
783	#endif
784	!
785	! Sea Surface Height at u-,v-points (vvl case only)
786	IF( .NOT.ln_linssh ) THEN
787	!$OMP PARALLEL DO schedule(static) private(jj,ji)
788	DO jj = 2, jpjm1
789	DO ji = 2, jpim1 ! NO Vector Opt.
790	zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &
791	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
792	& + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )
793	zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &
794	& * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &
795	& + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )
796	END DO
797	END DO
798	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
799	ENDIF
800	!
801	! Half-step back interpolation of SSH for surface pressure computation:
802	!----------------------------------------------------------------------
803	IF ((jn==1).AND.ll_init) THEN
804	za0=1._wp ! Forward-backward
805	za1=0._wp
806	za2=0._wp
807	za3=0._wp
808	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
809	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
810	za1=-0.1666666666666_wp ! za1 = gam
811	za2= 0.0833333333333_wp ! za2 = eps
812	za3= 0._wp
813	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
814	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
815	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
816	za2=0.088_wp ! za2 = gam
817	za3=0.013_wp ! za3 = eps
818	ENDIF
819	!
820	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
821	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
822	IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters
823	wduflt1(:,:) = 1._wp
824	wdvflt1(:,:) = 1._wp
825	!$OMP PARALLEL DO schedule(static) private(jj,ji,ll_tmp1,ll_tmp2)
826	DO jj = 2, jpjm1
827	DO ji = 2, jpim1
828	ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) &
829	& .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) ) &
830	& > rn_wdmin1 + rn_wdmin2
831	ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) &
832	& + rn_wdmin1 + rn_wdmin2
833	IF(ll_tmp1) THEN
834	zcpx(ji,jj) = 1._wp
835	ELSE IF(ll_tmp2) THEN
836	! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen here
837	zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) &
838	& / (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)) )
839	ELSE
840	zcpx(ji,jj) = 0._wp
841	wduflt1(ji,jj) = 0._wp
842	END IF
843
844	ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) &
845	& .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) ) &
846	& > rn_wdmin1 + rn_wdmin2
847	ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) &
848	& + rn_wdmin1 + rn_wdmin2
849	IF(ll_tmp1) THEN
850	zcpy(ji,jj) = 1._wp
851	ELSE IF(ll_tmp2) THEN
852	! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen here
853	zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) &
854	& / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)) )
855	ELSE
856	zcpy(ji,jj) = 0._wp
857	wdvflt1(ji,jj) = 0._wp
858	END IF
859	END DO
860	END DO
861	CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp )
862	ENDIF
863	!
864	! Compute associated depths at U and V points:
865	IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form
866	!
867	!$OMP PARALLEL DO schedule(static) private(jj,ji,zx1,zy1)
868	DO jj = 2, jpjm1
869	DO ji = 2, jpim1
870	zx1 = z1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) &
871	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) &
872	& + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
873	zy1 = z1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) &
874	& * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) &
875	& + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
876	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
877	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
878	END DO
879	END DO
880
881	IF( ln_wd ) THEN
882	zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 )
883	zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 )
884	END IF
885
886	ENDIF
887	!
888	! Add Coriolis trend:
889	! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
890	! at each time step. We however keep them constant here for optimization.
891	! Recall that zwx and zwy arrays hold fluxes at this stage:
892	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
893	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
894	!
895	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
896	!$OMP PARALLEL DO schedule(static) private(jj,ji,zy1,zy2,zx1,zx2)
897	DO jj = 2, jpjm1
898	DO ji = fs_2, fs_jpim1 ! vector opt.
899	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
900	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
901	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
902	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
903	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
904	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
905	END DO
906	END DO
907	!
908	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
909	!$OMP PARALLEL DO schedule(static) private(jj,ji,zx1,zy1)
910	DO jj = 2, jpjm1
911	DO ji = fs_2, fs_jpim1 ! vector opt.
912	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
913	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj)
914	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
915	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj)
916	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
917	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
918	END DO
919	END DO
920	!
921	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!
922	!$OMP PARALLEL DO schedule(static) private(jj,ji)
923	DO jj = 2, jpjm1
924	DO ji = fs_2, fs_jpim1 ! vector opt.
925	zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
926	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
927	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
928	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
929	zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
930	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
931	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
932	& + ftne(ji,jj ) * zwx(ji ,jj ) )
933	END DO
934	END DO
935	!
936	ENDIF
937	!
938	! Add tidal astronomical forcing if defined
939	IF ( lk_tide.AND.ln_tide_pot ) THEN
940	!$OMP PARALLEL DO schedule(static) private(jj,ji,zu_spg,zv_spg)
941	DO jj = 2, jpjm1
942	DO ji = fs_2, fs_jpim1 ! vector opt.
943	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
944	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
945	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
946	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
947	END DO
948	END DO
949	ENDIF
950	!
951	! Add bottom stresses:
952	!$OMP PARALLEL WORKSHARE
953	zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:)
954	zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:)
955	!$OMP END PARALLEL WORKSHARE
956	!
957	! Add top stresses:
958	!$OMP PARALLEL WORKSHARE
959	zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:)
960	zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:)
961	!$OMP END PARALLEL WORKSHARE
962	!
963	! Surface pressure trend:
964
965	IF( ln_wd ) THEN
966	!$OMP PARALLEL DO schedule(static) private(jj,ji,zu_spg,zv_spg)
967	DO jj = 2, jpjm1
968	DO ji = 2, jpim1
969	! Add surface pressure gradient
970	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
971	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
972	zwx(ji,jj) = zu_spg * zcpx(ji,jj)
973	zwy(ji,jj) = zv_spg * zcpy(ji,jj)
974	END DO
975	END DO
976	ELSE
977	!$OMP PARALLEL DO schedule(static) private(jj,ji,zu_spg,zv_spg)
978	DO jj = 2, jpjm1
979	DO ji = fs_2, fs_jpim1 ! vector opt.
980	! Add surface pressure gradient
981	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
982	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
983	zwx(ji,jj) = zu_spg
984	zwy(ji,jj) = zv_spg
985	END DO
986	END DO
987	END IF
988
989	!
990	! Set next velocities:
991	IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form
992	!$OMP PARALLEL DO schedule(static) private(jj,ji)
993	DO jj = 2, jpjm1
994	DO ji = fs_2, fs_jpim1 ! vector opt.
995	ua_e(ji,jj) = ( un_e(ji,jj) &
996	& + rdtbt * ( zwx(ji,jj) &
997	& + zu_trd(ji,jj) &
998	& + zu_frc(ji,jj) ) &
999	& ) * ssumask(ji,jj)
1000
1001	va_e(ji,jj) = ( vn_e(ji,jj) &
1002	& + rdtbt * ( zwy(ji,jj) &
1003	& + zv_trd(ji,jj) &
1004	& + zv_frc(ji,jj) ) &
1005	& ) * ssvmask(ji,jj)
1006	END DO
1007	END DO
1008	!
1009	ELSE !* Flux form
1010	!$OMP PARALLEL DO schedule(static) private(jj,ji,zhura,zhvra)
1011	DO jj = 2, jpjm1
1012	DO ji = fs_2, fs_jpim1 ! vector opt.
1013
1014	IF( ln_wd ) THEN
1015	zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1)
1016	zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1)
1017	ELSE
1018	zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
1019	zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
1020	END IF
1021	zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
1022	zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
1023
1024	ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) &
1025	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
1026	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
1027	& + hu_n(ji,jj) * zu_frc(ji,jj) ) &
1028	& ) * zhura
1029
1030	va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) &
1031	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
1032	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
1033	& + hv_n(ji,jj) * zv_frc(ji,jj) ) &
1034	& ) * zhvra
1035	END DO
1036	END DO
1037	ENDIF
1038	!
1039	IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
1040	IF( ln_wd ) THEN
1041	hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1)
1042	hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1)
1043	ELSE
1044	!$OMP PARALLEL WORKSHARE
1045	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
1046	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
1047	!$OMP END PARALLEL WORKSHARE
1048	END IF
1049	!$OMP PARALLEL WORKSHARE
1050	hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) )
1051	hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) )
1052	!$OMP END PARALLEL WORKSHARE
1053	!
1054	ENDIF
1055	! !* domain lateral boundary
1056	CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
1057	!
1058	#if defined key_bdy
1059	! ! open boundaries
1060	IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )
1061	#endif
1062	#if defined key_agrif
1063	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
1064	#endif
1065	! !* Swap
1066	! ! ----
1067	!$OMP PARALLEL WORKSHARE
1068	ubb_e (:,:) = ub_e (:,:)
1069	ub_e (:,:) = un_e (:,:)
1070	un_e (:,:) = ua_e (:,:)
1071	!
1072	vbb_e (:,:) = vb_e (:,:)
1073	vb_e (:,:) = vn_e (:,:)
1074	vn_e (:,:) = va_e (:,:)
1075	!
1076	sshbb_e(:,:) = sshb_e(:,:)
1077	sshb_e (:,:) = sshn_e(:,:)
1078	sshn_e (:,:) = ssha_e(:,:)
1079	!$OMP END PARALLEL WORKSHARE
1080
1081	! !* Sum over whole bt loop
1082	! ! ----------------------
1083	za1 = wgtbtp1(jn)
1084	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
1085	!$OMP PARALLEL WORKSHARE
1086	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
1087	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
1088	!$OMP END PARALLEL WORKSHARE
1089	ELSE ! Sum transports
1090	!$OMP PARALLEL WORKSHARE
1091	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
1092	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
1093	!$OMP END PARALLEL WORKSHARE
1094	ENDIF
1095	! ! Sum sea level
1096	!$OMP PARALLEL WORKSHARE
1097	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
1098	!$OMP END PARALLEL WORKSHARE
1099	! ! ==================== !
1100	END DO ! end loop !
1101	! ! ==================== !
1102	! -----------------------------------------------------------------------------
1103	! Phase 3. update the general trend with the barotropic trend
1104	! -----------------------------------------------------------------------------
1105	!
1106	! Set advection velocity correction:
1107	!$OMP PARALLEL WORKSHARE
1108	zwx(:,:) = un_adv(:,:)
1109	zwy(:,:) = vn_adv(:,:)
1110	!$OMP END PARALLEL WORKSHARE
1111	IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN
1112	!$OMP PARALLEL WORKSHARE
1113	un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:)
1114	vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:)
1115	!$OMP END PARALLEL WORKSHARE
1116	ELSE
1117	!$OMP PARALLEL WORKSHARE
1118	un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:)
1119	vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:)
1120	!$OMP END PARALLEL WORKSHARE
1121	END IF
1122
1123	IF( ln_bt_fw ) THEN ! Save integrated transport for next computation
1124	!$OMP PARALLEL WORKSHARE
1125	ub2_b(:,:) = zwx(:,:)
1126	vb2_b(:,:) = zwy(:,:)
1127	!$OMP END PARALLEL WORKSHARE
1128	ENDIF
1129	!
1130	! Update barotropic trend:
1131	IF( ln_dynadv_vec .OR. ln_linssh ) THEN
1132	!$OMP PARALLEL DO schedule(static) private(jk)
1133	DO jk=1,jpkm1
1134	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
1135	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
1136	END DO
1137	ELSE
1138	! At this stage, ssha has been corrected: compute new depths at velocity points
1139	!$OMP PARALLEL DO schedule(static) private(jj,ji)
1140	DO jj = 1, jpjm1
1141	DO ji = 1, jpim1 ! NO Vector Opt.
1142	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &
1143	& * ( e1e2t(ji ,jj) * ssha(ji ,jj) &
1144	& + e1e2t(ji+1,jj) * ssha(ji+1,jj) )
1145	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &
1146	& * ( e1e2t(ji,jj ) * ssha(ji,jj ) &
1147	& + e1e2t(ji,jj+1) * ssha(ji,jj+1) )
1148	END DO
1149	END DO
1150	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
1151	!
1152	!$OMP PARALLEL DO schedule(static) private(jk)
1153	DO jk=1,jpkm1
1154	ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
1155	va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
1156	END DO
1157	! Save barotropic velocities not transport:
1158	!$OMP PARALLEL WORKSHARE
1159	ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
1160	va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
1161	!$OMP END PARALLEL WORKSHARE
1162	ENDIF
1163	!
1164	!$OMP PARALLEL DO schedule(static) private(jk)
1165	DO jk = 1, jpkm1
1166	! Correct velocities:
1167	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk)
1168	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk)
1169	!
1170	END DO
1171	!
1172	CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current
1173	CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic i-current
1174	!
1175	#if defined key_agrif
1176	! Save time integrated fluxes during child grid integration
1177	! (used to update coarse grid transports at next time step)
1178	!
1179	IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN
1180	IF( Agrif_NbStepint() == 0 ) THEN
1181	ub2_i_b(:,:) = 0._wp
1182	vb2_i_b(:,:) = 0._wp
1183	END IF
1184	!
1185	za1 = 1._wp / REAL(Agrif_rhot(), wp)
1186	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
1187	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
1188	ENDIF
1189	#endif
1190	! !* write time-spliting arrays in the restart
1191	IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
1192	!
1193	CALL wrk_dealloc( jpi,jpj, zsshp2_e, zhdiv )
1194	CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd )
1195	CALL wrk_dealloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc )
1196	CALL wrk_dealloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
1197	CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a )
1198	CALL wrk_dealloc( jpi,jpj, zhf )
1199	IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 )
1200	!
1201	IF ( ln_diatmb ) THEN
1202	CALL iom_put( "baro_u" , un_bumask(:,:,1)+zmdi(1-umask(:,:,1 ) ) ) ! Barotropic U Velocity
1203	CALL iom_put( "baro_v" , vn_bvmask(:,:,1)+zmdi(1-vmask(:,:,1 ) ) ) ! Barotropic V Velocity
1204	ENDIF
1205	IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts')
1206	!
1207	END SUBROUTINE dyn_spg_ts
1208
1209
1210	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
1211	!!---------------------------------------------------------------------
1212	!! * ROUTINE ts_wgt *
1213	!!
1214	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
1215	!!----------------------------------------------------------------------
1216	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
1217	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
1218	INTEGER, INTENT(inout) :: jpit ! cycle length
1219	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
1220	zwgt2 ! Secondary weights
1221
1222	INTEGER :: jic, jn, ji ! temporary integers
1223	REAL(wp) :: za1, za2
1224	!!----------------------------------------------------------------------
1225
1226	zwgt1(:) = 0._wp
1227	zwgt2(:) = 0._wp
1228
1229	! Set time index when averaged value is requested
1230	IF (ll_fw) THEN
1231	jic = nn_baro
1232	ELSE
1233	jic = 2 * nn_baro
1234	ENDIF
1235
1236	! Set primary weights:
1237	IF (ll_av) THEN
1238	! Define simple boxcar window for primary weights
1239	! (width = nn_baro, centered around jic)
1240	SELECT CASE ( nn_bt_flt )
1241	CASE( 0 ) ! No averaging
1242	zwgt1(jic) = 1._wp
1243	jpit = jic
1244
1245	CASE( 1 ) ! Boxcar, width = nn_baro
1246	DO jn = 1, 3*nn_baro
1247	za1 = ABS(float(jn-jic))/float(nn_baro)
1248	IF (za1 < 0.5_wp) THEN
1249	zwgt1(jn) = 1._wp
1250	jpit = jn
1251	ENDIF
1252	ENDDO
1253
1254	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
1255	DO jn = 1, 3*nn_baro
1256	za1 = ABS(float(jn-jic))/float(nn_baro)
1257	IF (za1 < 1._wp) THEN
1258	zwgt1(jn) = 1._wp
1259	jpit = jn
1260	ENDIF
1261	ENDDO
1262	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
1263	END SELECT
1264
1265	ELSE ! No time averaging
1266	zwgt1(jic) = 1._wp
1267	jpit = jic
1268	ENDIF
1269
1270	! Set secondary weights
1271	DO jn = 1, jpit
1272	DO ji = jn, jpit
1273	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1274	END DO
1275	END DO
1276
1277	! Normalize weigths:
1278	za1 = 1._wp / SUM(zwgt1(1:jpit))
1279	za2 = 1._wp / SUM(zwgt2(1:jpit))
1280	DO jn = 1, jpit
1281	zwgt1(jn) = zwgt1(jn) * za1
1282	zwgt2(jn) = zwgt2(jn) * za2
1283	END DO
1284	!
1285	END SUBROUTINE ts_wgt
1286
1287
1288	SUBROUTINE ts_rst( kt, cdrw )
1289	!!---------------------------------------------------------------------
1290	!! * ROUTINE ts_rst *
1291	!!
1292	!! ** Purpose : Read or write time-splitting arrays in restart file
1293	!!----------------------------------------------------------------------
1294	INTEGER , INTENT(in) :: kt ! ocean time-step
1295	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1296	!
1297	!!----------------------------------------------------------------------
1298	!
1299	IF( TRIM(cdrw) == 'READ' ) THEN
1300	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) )
1301	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) )
1302	IF( .NOT.ln_bt_av ) THEN
1303	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) )
1304	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) )
1305	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) )
1306	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) )
1307	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) )
1308	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) )
1309	ENDIF
1310	#if defined key_agrif
1311	! Read time integrated fluxes
1312	IF ( .NOT.Agrif_Root() ) THEN
1313	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:) )
1314	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:) )
1315	ENDIF
1316	#endif
1317	!
1318	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1319	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
1320	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
1321	!
1322	IF (.NOT.ln_bt_av) THEN
1323	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
1324	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
1325	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
1326	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
1327	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
1328	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
1329	ENDIF
1330	#if defined key_agrif
1331	! Save time integrated fluxes
1332	IF ( .NOT.Agrif_Root() ) THEN
1333	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
1334	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
1335	ENDIF
1336	#endif
1337	ENDIF
1338	!
1339	END SUBROUTINE ts_rst
1340
1341
1342	SUBROUTINE dyn_spg_ts_init
1343	!!---------------------------------------------------------------------
1344	!! * ROUTINE dyn_spg_ts_init *
1345	!!
1346	!! ** Purpose : Set time splitting options
1347	!!----------------------------------------------------------------------
1348	INTEGER :: ji ,jj ! dummy loop indices
1349	REAL(wp) :: zxr2, zyr2, zcmax ! local scalar
1350	REAL(wp), POINTER, DIMENSION(:,:) :: zcu
1351	!!----------------------------------------------------------------------
1352	!
1353	! Max courant number for ext. grav. waves
1354	!
1355	CALL wrk_alloc( jpi,jpj, zcu )
1356	!
1357	!$OMP PARALLEL DO schedule(static) private(jj, ji, zxr2, zyr2)
1358	DO jj = 1, jpj
1359	DO ji =1, jpi
1360	zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
1361	zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
1362	zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) )
1363	END DO
1364	END DO
1365	!
1366	zcmax = MAXVAL( zcu(:,:) )
1367	IF( lk_mpp ) CALL mpp_max( zcmax )
1368
1369	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1370	IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1371
1372	rdtbt = rdt / REAL( nn_baro , wp )
1373	zcmax = zcmax * rdtbt
1374	! Print results
1375	IF(lwp) WRITE(numout,*)
1376	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1377	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1378	IF( ln_bt_auto ) THEN
1379	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro '
1380	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1381	ELSE
1382	IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist '
1383	ENDIF
1384
1385	IF(ln_bt_av) THEN
1386	IF(lwp) WRITE(numout,*) ' ln_bt_av=.true. => Time averaging over nn_baro time steps is on '
1387	ELSE
1388	IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
1389	ENDIF
1390	!
1391	!
1392	IF(ln_bt_fw) THEN
1393	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1394	ELSE
1395	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1396	ENDIF
1397	!
1398	#if defined key_agrif
1399	! Restrict the use of Agrif to the forward case only
1400	IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1401	#endif
1402	!
1403	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1404	SELECT CASE ( nn_bt_flt )
1405	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1406	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1407	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1408	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1409	END SELECT
1410	!
1411	IF(lwp) WRITE(numout,*) ' '
1412	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1413	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1414	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1415	!
1416	IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
1417	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1418	ENDIF
1419	IF( zcmax>0.9_wp ) THEN
1420	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1421	ENDIF
1422	!
1423	CALL wrk_dealloc( jpi,jpj, zcu )
1424	!
1425	END SUBROUTINE dyn_spg_ts_init
1426
1427	!!======================================================================
1428	END MODULE dynspg_ts

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

Download in other formats: