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/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7167

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

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

Download in other formats: