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

Context Navigation

source: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5944

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

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

Download in other formats: