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 @ 5956

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

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

Download in other formats: