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

Context Navigation

source: branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4960

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

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

Download in other formats: