New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynspg_ts.F90 in branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5989

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

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

Download in other formats: