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

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

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

Download in other formats: