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

Context Navigation

source: branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4223

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

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

Download in other formats: