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

Context Navigation

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 6035

Last change on this file since 6035 was 6035, checked in by timgraham, 8 years ago
Merged branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics into branch
Property svn:keywords set to `Id`
File size: 55.0 KB

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

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

Download in other formats: