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

Context Navigation

source: branches/2015/dev_5894_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 6876

Last change on this file since 6876 was 5899, checked in by emanuelaclementi, 9 years ago
ticket 1632 Add in changes from the 2014/dev_4822_INGV_WAVE branch
Property svn:keywords set to `Id`
File size: 55.3 KB

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

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

Download in other formats: