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

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4496

Last change on this file since 4496 was 4496, checked in by jchanut, 10 years ago

Fix misplaced call to ts_rst
o-This line, and those below, will be ignored--

M NEMO/OPA_SRC/DYN/dynspg.F90
M NEMO/OPA_SRC/DYN/dynspg_ts.F90

Property svn:keywords set to Id

File size: 54.7 KB

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

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

Download in other formats: