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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7698

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

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

Download in other formats: