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

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4487

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

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

Download in other formats: