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

Context Navigation

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5866

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

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

Download in other formats: