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 @ 7037

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

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

Download in other formats: