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

Context Navigation

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago
#2050 first version. Compiled OK in moci test suite
File size: 57.0 KB

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

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

Download in other formats: