Context Navigation

dynspg_ts.F90 @ 9466

Last change on this file since 9466 was 9466, checked in by jcastill, 6 years ago
Move again some lines to allow merging with operational branches
File size: 56.3 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	!! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence
14	!!---------------------------------------------------------------------
15	#if defined key_dynspg_ts \|\| defined key_esopa
16	!!----------------------------------------------------------------------
17	!! 'key_dynspg_ts' split explicit free surface
18	!!----------------------------------------------------------------------
19	!! dyn_spg_ts : compute surface pressure gradient trend using a time-
20	!! splitting scheme and add to the general trend
21	!!----------------------------------------------------------------------
22	USE oce ! ocean dynamics and tracers
23	USE dom_oce ! ocean space and time domain
24	USE sbc_oce ! surface boundary condition: ocean
25	USE sbcisf ! ice shelf variable (fwfisf)
26	USE dynspg_oce ! surface pressure gradient variables
27	USE phycst ! physical constants
28	USE dynvor ! vorticity term
29	USE bdy_par ! for lk_bdy
30	USE bdytides ! open boundary condition data
31	USE bdydyn2d ! open boundary conditions on barotropic variables
32	USE sbctide ! tides
33	USE updtide ! tide potential
34	USE sbcwave ! surface wave
35	!
36	USE lib_mpp ! distributed memory computing library
37	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
38	USE prtctl ! Print control
39	USE in_out_manager ! I/O manager
40	USE iom ! IOM library
41	USE restart ! only for lrst_oce
42	USE zdf_oce ! Bottom friction coefts
43	USE wrk_nemo ! Memory Allocation
44	USE timing ! Timing
45	USE sbcapr ! surface boundary condition: atmospheric pressure
46	USE dynadv, ONLY: ln_dynadv_vec
47	#if defined key_agrif
48	USE agrif_opa_interp ! agrif
49	#endif
50	#if defined key_asminc
51	USE asminc ! Assimilation increment
52	#endif
53
54	IMPLICIT NONE
55	PRIVATE
56
57	PUBLIC dyn_spg_ts ! routine called in dynspg.F90
58	PUBLIC dyn_spg_ts_alloc ! " " " "
59	PUBLIC dyn_spg_ts_init ! " " " "
60	PUBLIC ts_rst ! " " " "
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
336	!!gm Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum....
337	!!gm Is it correct to do so ? I think so...
338
339
340	! !* barotropic Coriolis trends (vorticity scheme dependent)
341	! ! --------------------------------------------------------
342	zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:) ! now fluxes
343	zwy(:,:) = vn_b(:,:) * hv(:,:) * e1v(:,:)
344	!
345	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
346	DO jj = 2, jpjm1
347	DO ji = fs_2, fs_jpim1 ! vector opt.
348	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
349	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
350	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
351	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
352	! energy conserving formulation for planetary vorticity term
353	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
354	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
355	END DO
356	END DO
357	!
358	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
359	DO jj = 2, jpjm1
360	DO ji = fs_2, fs_jpim1 ! vector opt.
361	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
362	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
363	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
364	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
365	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
366	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
367	END DO
368	END DO
369	!
370	ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN ! enstrophy and energy conserving scheme
371	DO jj = 2, jpjm1
372	DO ji = fs_2, fs_jpim1 ! vector opt.
373	zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
374	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
375	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
376	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
377	zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
378	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
379	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
380	& + ftne(ji,jj ) * zwx(ji ,jj ) )
381	END DO
382	END DO
383	!
384	ENDIF
385	!
386	! !* Right-Hand-Side of the barotropic momentum equation
387	! ! ----------------------------------------------------
388	IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient
389	DO jj = 2, jpjm1
390	DO ji = fs_2, fs_jpim1 ! vector opt.
391	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj)
392	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj)
393	END DO
394	END DO
395	ENDIF
396
397	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
398	DO ji = fs_2, fs_jpim1
399	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)
400	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
401	END DO
402	END DO
403	!
404	! ! Add bottom stress contribution from baroclinic velocities:
405	IF (ln_bt_fw) THEN
406	DO jj = 2, jpjm1
407	DO ji = fs_2, fs_jpim1 ! vector opt.
408	ikbu = mbku(ji,jj)
409	ikbv = mbkv(ji,jj)
410	zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
411	zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
412	END DO
413	END DO
414	ELSE
415	DO jj = 2, jpjm1
416	DO ji = fs_2, fs_jpim1 ! vector opt.
417	ikbu = mbku(ji,jj)
418	ikbv = mbkv(ji,jj)
419	zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
420	zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
421	END DO
422	END DO
423	ENDIF
424	!
425	! Note that the "unclipped" bottom friction parameter is used even with explicit drag
426	zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
427	zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
428	!
429	IF (ln_bt_fw) THEN ! Add wind forcing
430	zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
431	zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)
432	ELSE
433	zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)
434	zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
435	ENDIF
436	!
437	IF ( ln_apr_dyn ) THEN ! Add atm pressure forcing
438	IF (ln_bt_fw) THEN
439	DO jj = 2, jpjm1
440	DO ji = fs_2, fs_jpim1 ! vector opt.
441	zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
442	zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)
443	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
444	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
445	END DO
446	END DO
447	ELSE
448	DO jj = 2, jpjm1
449	DO ji = fs_2, fs_jpim1 ! vector opt.
450	zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
451	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj)
452	zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
453	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)
454	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
455	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
456	END DO
457	END DO
458	ENDIF
459	ENDIF
460	! !* Right-Hand-Side of the barotropic ssh equation
461	! ! -----------------------------------------------
462	! ! Surface net water flux and rivers
463	IF (ln_bt_fw) THEN
464	zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
465	ELSE
466	zssh_frc(:,:) = zraur * z1_2 * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) &
467	& + fwfisf(:,:) + fwfisf_b(:,:) )
468	ENDIF
469	#if defined key_asminc
470	! ! Include the IAU weighted SSH increment
471	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
472	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
473	ENDIF
474	#endif
475	! !* Fill boundary data arrays for AGRIF
476	! ! ------------------------------------
477	#if defined key_agrif
478	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
479	#endif
480	!
481	! -----------------------------------------------------------------------
482	! Phase 2 : Integration of the barotropic equations
483	! -----------------------------------------------------------------------
484	!
485	! ! ==================== !
486	! ! Initialisations !
487	! ! ==================== !
488	! Initialize barotropic variables:
489	IF( ll_init )THEN
490	sshbb_e(:,:) = 0._wp
491	ubb_e (:,:) = 0._wp
492	vbb_e (:,:) = 0._wp
493	sshb_e (:,:) = 0._wp
494	ub_e (:,:) = 0._wp
495	vb_e (:,:) = 0._wp
496	ENDIF
497	!
498	! Moved here to allow merging with other branch
499	IF( ln_sdw ) THEN ! Stokes drift divergence added if necessary
500	zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
501	ENDIF
502	!
503	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
504	sshn_e(:,:) = sshn (:,:)
505	zun_e (:,:) = un_b (:,:)
506	zvn_e (:,:) = vn_b (:,:)
507	!
508	hu_e (:,:) = hu (:,:)
509	hv_e (:,:) = hv (:,:)
510	hur_e (:,:) = hur (:,:)
511	hvr_e (:,:) = hvr (:,:)
512	ELSE ! CENTRED integration: start from BEFORE fields
513	sshn_e(:,:) = sshb (:,:)
514	zun_e (:,:) = ub_b (:,:)
515	zvn_e (:,:) = vb_b (:,:)
516	!
517	hu_e (:,:) = hu_b (:,:)
518	hv_e (:,:) = hv_b (:,:)
519	hur_e (:,:) = hur_b(:,:)
520	hvr_e (:,:) = hvr_b(:,:)
521	ENDIF
522	!
523	!
524	!
525	! Initialize sums:
526	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
527	va_b (:,:) = 0._wp
528	ssha (:,:) = 0._wp ! Sum for after averaged sea level
529	zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop
530	zv_sum(:,:) = 0._wp
531	! ! ==================== !
532	DO jn = 1, icycle ! sub-time-step loop !
533	! ! ==================== !
534	! !* Update the forcing (BDY and tides)
535	! ! ------------------
536	! Update only tidal forcing at open boundaries
537	#if defined key_tide
538	IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
539	IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset )
540	#endif
541	!
542	! Set extrapolation coefficients for predictor step:
543	IF ((jn<3).AND.ll_init) THEN ! Forward
544	za1 = 1._wp
545	za2 = 0._wp
546	za3 = 0._wp
547	ELSE ! AB3-AM4 Coefficients: bet=0.281105
548	za1 = 1.781105_wp ! za1 = 3/2 + bet
549	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
550	za3 = 0.281105_wp ! za3 = bet
551	ENDIF
552
553	! Extrapolate barotropic velocities at step jit+0.5:
554	ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
555	va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
556
557	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
558	! ! ------------------
559	! Extrapolate Sea Level at step jit+0.5:
560	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
561	!
562	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
563	DO ji = 2, fs_jpim1 ! Vector opt.
564	zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
565	& * ( e12t(ji ,jj) * zsshp2_e(ji ,jj) &
566	& + e12t(ji+1,jj) * zsshp2_e(ji+1,jj) )
567	zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
568	& * ( e12t(ji,jj ) * zsshp2_e(ji,jj ) &
569	& + e12t(ji,jj+1) * zsshp2_e(ji,jj+1) )
570	END DO
571	END DO
572	CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
573	!
574	zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
575	zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
576	ELSE
577	zhup2_e (:,:) = hu(:,:)
578	zhvp2_e (:,:) = hv(:,:)
579	ENDIF
580	! !* after ssh
581	! ! -----------
582	! One should enforce volume conservation at open boundaries here
583	! considering fluxes below:
584	!
585	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
586	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
587	!
588	#if defined key_agrif
589	! Set fluxes during predictor step to ensure
590	! volume conservation
591	IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN
592	IF((nbondi == -1).OR.(nbondi == 2)) THEN
593	DO jj=1,jpj
594	zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
595	END DO
596	ENDIF
597	IF((nbondi == 1).OR.(nbondi == 2)) THEN
598	DO jj=1,jpj
599	zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
600	END DO
601	ENDIF
602	IF((nbondj == -1).OR.(nbondj == 2)) THEN
603	DO ji=1,jpi
604	zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
605	END DO
606	ENDIF
607	IF((nbondj == 1).OR.(nbondj == 2)) THEN
608	DO ji=1,jpi
609	zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
610	END DO
611	ENDIF
612	ENDIF
613	#endif
614	!
615	! Sum over sub-time-steps to compute advective velocities
616	za2 = wgtbtp2(jn)
617	zu_sum (:,:) = zu_sum (:,:) + za2 * zwx (:,:) / e2u (:,:)
618	zv_sum (:,:) = zv_sum (:,:) + za2 * zwy (:,:) / e1v (:,:)
619	!
620	! Set next sea level:
621	DO jj = 2, jpjm1
622	DO ji = fs_2, fs_jpim1 ! vector opt.
623	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
624	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e12t(ji,jj)
625	END DO
626	END DO
627	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)
628	CALL lbc_lnk( ssha_e, 'T', 1._wp )
629
630	#if defined key_bdy
631	! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
632	IF (lk_bdy) CALL bdy_ssh( ssha_e )
633	#endif
634	#if defined key_agrif
635	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
636	#endif
637	!
638	! Sea Surface Height at u-,v-points (vvl case only)
639	IF ( lk_vvl ) THEN
640	DO jj = 2, jpjm1
641	DO ji = 2, jpim1 ! NO Vector Opt.
642	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
643	& * ( e12t(ji ,jj ) * ssha_e(ji ,jj ) &
644	& + e12t(ji+1,jj ) * ssha_e(ji+1,jj ) )
645	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
646	& * ( e12t(ji ,jj ) * ssha_e(ji ,jj ) &
647	& + e12t(ji ,jj+1) * ssha_e(ji ,jj+1) )
648	END DO
649	END DO
650	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
651	ENDIF
652	!
653	! Half-step back interpolation of SSH for surface pressure computation:
654	!----------------------------------------------------------------------
655	IF ((jn==1).AND.ll_init) THEN
656	za0=1._wp ! Forward-backward
657	za1=0._wp
658	za2=0._wp
659	za3=0._wp
660	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
661	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
662	za1=-0.1666666666666_wp ! za1 = gam
663	za2= 0.0833333333333_wp ! za2 = eps
664	za3= 0._wp
665	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
666	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
667	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
668	za2=0.088_wp ! za2 = gam
669	za3=0.013_wp ! za3 = eps
670	ENDIF
671
672	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
673	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
674
675	!
676	! Compute associated depths at U and V points:
677	IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN
678	!
679	DO jj = 2, jpjm1
680	DO ji = 2, jpim1
681	zx1 = z1_2 * umask(ji ,jj,1) * r1_e12u(ji ,jj) &
682	& * ( e12t(ji ,jj ) * zsshp2_e(ji ,jj) &
683	& + e12t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
684	zy1 = z1_2 * vmask(ji ,jj,1) * r1_e12v(ji ,jj ) &
685	& * ( e12t(ji ,jj ) * zsshp2_e(ji ,jj ) &
686	& + e12t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
687	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
688	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
689	END DO
690	END DO
691	ENDIF
692	!
693	! Add Coriolis trend:
694	! zwz array below or triads normally depend on sea level with key_vvl and should be updated
695	! at each time step. We however keep them constant here for optimization.
696	! Recall that zwx and zwy arrays hold fluxes at this stage:
697	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
698	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
699	!
700	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
701	DO jj = 2, jpjm1
702	DO ji = fs_2, fs_jpim1 ! vector opt.
703	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
704	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
705	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
706	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
707	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
708	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
709	END DO
710	END DO
711	!
712	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
713	DO jj = 2, jpjm1
714	DO ji = fs_2, fs_jpim1 ! vector opt.
715	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
716	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
717	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
718	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
719	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
720	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
721	END DO
722	END DO
723	!
724	ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !== energy and enstrophy conserving scheme ==!
725	DO jj = 2, jpjm1
726	DO ji = fs_2, fs_jpim1 ! vector opt.
727	zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
728	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
729	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
730	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
731	zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
732	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
733	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
734	& + ftne(ji,jj ) * zwx(ji ,jj ) )
735	END DO
736	END DO
737	!
738	ENDIF
739	!
740	! Add tidal astronomical forcing if defined
741	IF ( lk_tide.AND.ln_tide_pot ) THEN
742	DO jj = 2, jpjm1
743	DO ji = fs_2, fs_jpim1 ! vector opt.
744	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
745	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
746	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
747	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
748	END DO
749	END DO
750	ENDIF
751	!
752	! Add bottom stresses:
753	zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
754	zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
755	!
756	! Surface pressure trend:
757	DO jj = 2, jpjm1
758	DO ji = fs_2, fs_jpim1 ! vector opt.
759	! Add surface pressure gradient
760	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
761	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
762	zwx(ji,jj) = zu_spg
763	zwy(ji,jj) = zv_spg
764	END DO
765	END DO
766	!
767	! Set next velocities:
768	IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN ! Vector form
769	DO jj = 2, jpjm1
770	DO ji = fs_2, fs_jpim1 ! vector opt.
771	ua_e(ji,jj) = ( zun_e(ji,jj) &
772	& + rdtbt * ( zwx(ji,jj) &
773	& + zu_trd(ji,jj) &
774	& + zu_frc(ji,jj) ) &
775	& ) * umask(ji,jj,1)
776
777	va_e(ji,jj) = ( zvn_e(ji,jj) &
778	& + rdtbt * ( zwy(ji,jj) &
779	& + zv_trd(ji,jj) &
780	& + zv_frc(ji,jj) ) &
781	& ) * vmask(ji,jj,1)
782	END DO
783	END DO
784
785	ELSE ! Flux form
786	DO jj = 2, jpjm1
787	DO ji = fs_2, fs_jpim1 ! vector opt.
788
789	zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
790	zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
791
792	ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) &
793	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
794	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
795	& + hu(ji,jj) * zu_frc(ji,jj) ) &
796	& ) * zhura
797
798	va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &
799	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
800	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
801	& + hv(ji,jj) * zv_frc(ji,jj) ) &
802	& ) * zhvra
803	END DO
804	END DO
805	ENDIF
806	!
807	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
808	! ! ----------------------------------------------
809	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
810	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
811	hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
812	hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
813	!
814	ENDIF
815	! !* domain lateral boundary
816	! ! -----------------------
817	!
818	CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
819
820	#if defined key_bdy
821	! open boundaries
822	IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )
823	#endif
824	#if defined key_agrif
825	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
826	#endif
827	! !* Swap
828	! ! ----
829	ubb_e (:,:) = ub_e (:,:)
830	ub_e (:,:) = zun_e (:,:)
831	zun_e (:,:) = ua_e (:,:)
832	!
833	vbb_e (:,:) = vb_e (:,:)
834	vb_e (:,:) = zvn_e (:,:)
835	zvn_e (:,:) = va_e (:,:)
836	!
837	sshbb_e(:,:) = sshb_e(:,:)
838	sshb_e (:,:) = sshn_e(:,:)
839	sshn_e (:,:) = ssha_e(:,:)
840
841	! !* Sum over whole bt loop
842	! ! ----------------------
843	za1 = wgtbtp1(jn)
844	IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities
845	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
846	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
847	ELSE ! Sum transports
848	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
849	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
850	ENDIF
851	! ! Sum sea level
852	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
853	! ! ==================== !
854	END DO ! end loop !
855	! ! ==================== !
856	! -----------------------------------------------------------------------------
857	! Phase 3. update the general trend with the barotropic trend
858	! -----------------------------------------------------------------------------
859	!
860	! At this stage ssha holds a time averaged value
861	! ! Sea Surface Height at u-,v- and f-points
862	IF( lk_vvl ) THEN ! (required only in key_vvl case)
863	DO jj = 1, jpjm1
864	DO ji = 1, jpim1 ! NO Vector Opt.
865	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
866	& * ( e12t(ji ,jj) * ssha(ji ,jj) &
867	& + e12t(ji+1,jj) * ssha(ji+1,jj) )
868	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
869	& * ( e12t(ji,jj ) * ssha(ji,jj ) &
870	& + e12t(ji,jj+1) * ssha(ji,jj+1) )
871	END DO
872	END DO
873	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
874	ENDIF
875	!
876	! Set advection velocity correction:
877	IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN
878	un_adv(:,:) = zu_sum(:,:)*hur(:,:)
879	vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
880	ELSE
881	un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
882	vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
883	END IF
884
885	IF (ln_bt_fw) THEN ! Save integrated transport for next computation
886	ub2_b(:,:) = zu_sum(:,:)
887	vb2_b(:,:) = zv_sum(:,:)
888	ENDIF
889	!
890	! Update barotropic trend:
891	IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
892	DO jk=1,jpkm1
893	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
894	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
895	END DO
896	ELSE
897	DO jk=1,jpkm1
898	ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
899	va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
900	END DO
901	! Save barotropic velocities not transport:
902	ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
903	va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
904	ENDIF
905	!
906	DO jk = 1, jpkm1
907	! Correct velocities:
908	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
909	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
910	!
911	END DO
912	!
913	#if defined key_agrif
914	! Save time integrated fluxes during child grid integration
915	! (used to update coarse grid transports at next time step)
916	!
917	IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN
918	IF ( Agrif_NbStepint().EQ.0 ) THEN
919	ub2_i_b(:,:) = 0.e0
920	vb2_i_b(:,:) = 0.e0
921	END IF
922	!
923	za1 = 1._wp / REAL(Agrif_rhot(), wp)
924	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
925	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
926	ENDIF
927	!
928	!
929	#endif
930	!
931	! !* write time-spliting arrays in the restart
932	IF(lrst_oce .AND.ln_bt_fw) CALL ts_rst( kt, 'WRITE' )
933	!
934	CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv )
935	CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
936	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
937	CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
938	CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a )
939	CALL wrk_dealloc( jpi, jpj, zhf )
940	!
941	IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts')
942	!
943	END SUBROUTINE dyn_spg_ts
944
945	SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
946	!!---------------------------------------------------------------------
947	!! * ROUTINE ts_wgt *
948	!!
949	!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
950	!!----------------------------------------------------------------------
951	LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
952	LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
953	INTEGER, INTENT(inout) :: jpit ! cycle length
954	REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights
955	zwgt2 ! Secondary weights
956
957	INTEGER :: jic, jn, ji ! temporary integers
958	REAL(wp) :: za1, za2
959	!!----------------------------------------------------------------------
960
961	zwgt1(:) = 0._wp
962	zwgt2(:) = 0._wp
963
964	! Set time index when averaged value is requested
965	IF (ll_fw) THEN
966	jic = nn_baro
967	ELSE
968	jic = 2 * nn_baro
969	ENDIF
970
971	! Set primary weights:
972	IF (ll_av) THEN
973	! Define simple boxcar window for primary weights
974	! (width = nn_baro, centered around jic)
975	SELECT CASE ( nn_bt_flt )
976	CASE( 0 ) ! No averaging
977	zwgt1(jic) = 1._wp
978	jpit = jic
979
980	CASE( 1 ) ! Boxcar, width = nn_baro
981	DO jn = 1, 3*nn_baro
982	za1 = ABS(float(jn-jic))/float(nn_baro)
983	IF (za1 < 0.5_wp) THEN
984	zwgt1(jn) = 1._wp
985	jpit = jn
986	ENDIF
987	ENDDO
988
989	CASE( 2 ) ! Boxcar, width = 2 * nn_baro
990	DO jn = 1, 3*nn_baro
991	za1 = ABS(float(jn-jic))/float(nn_baro)
992	IF (za1 < 1._wp) THEN
993	zwgt1(jn) = 1._wp
994	jpit = jn
995	ENDIF
996	ENDDO
997	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
998	END SELECT
999
1000	ELSE ! No time averaging
1001	zwgt1(jic) = 1._wp
1002	jpit = jic
1003	ENDIF
1004
1005	! Set secondary weights
1006	DO jn = 1, jpit
1007	DO ji = jn, jpit
1008	zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
1009	END DO
1010	END DO
1011
1012	! Normalize weigths:
1013	za1 = 1._wp / SUM(zwgt1(1:jpit))
1014	za2 = 1._wp / SUM(zwgt2(1:jpit))
1015	DO jn = 1, jpit
1016	zwgt1(jn) = zwgt1(jn) * za1
1017	zwgt2(jn) = zwgt2(jn) * za2
1018	END DO
1019	!
1020	END SUBROUTINE ts_wgt
1021
1022	SUBROUTINE ts_rst( kt, cdrw )
1023	!!---------------------------------------------------------------------
1024	!! * ROUTINE ts_rst *
1025	!!
1026	!! ** Purpose : Read or write time-splitting arrays in restart file
1027	!!----------------------------------------------------------------------
1028	INTEGER , INTENT(in) :: kt ! ocean time-step
1029	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1030	!
1031	!!----------------------------------------------------------------------
1032	!
1033	IF( TRIM(cdrw) == 'READ' ) THEN
1034	CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) )
1035	CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) )
1036	IF( .NOT.ln_bt_av ) THEN
1037	CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) )
1038	CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) )
1039	CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) )
1040	CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) )
1041	CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) )
1042	CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) )
1043	ENDIF
1044	#if defined key_agrif
1045	! Read time integrated fluxes
1046	IF ( .NOT.Agrif_Root() ) THEN
1047	CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:) )
1048	CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:) )
1049	ENDIF
1050	#endif
1051	!
1052	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
1053	CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
1054	CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
1055	!
1056	IF (.NOT.ln_bt_av) THEN
1057	CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
1058	CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) )
1059	CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) )
1060	CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) )
1061	CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) )
1062	CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) )
1063	ENDIF
1064	#if defined key_agrif
1065	! Save time integrated fluxes
1066	IF ( .NOT.Agrif_Root() ) THEN
1067	CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
1068	CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
1069	ENDIF
1070	#endif
1071	ENDIF
1072	!
1073	END SUBROUTINE ts_rst
1074
1075	SUBROUTINE dyn_spg_ts_init( kt )
1076	!!---------------------------------------------------------------------
1077	!! * ROUTINE dyn_spg_ts_init *
1078	!!
1079	!! ** Purpose : Set time splitting options
1080	!!----------------------------------------------------------------------
1081	INTEGER , INTENT(in) :: kt ! ocean time-step
1082	!
1083	INTEGER :: ji ,jj
1084	INTEGER :: ios ! Local integer output status for namelist read
1085	REAL(wp) :: zxr2, zyr2, zcmax
1086	REAL(wp), POINTER, DIMENSION(:,:) :: zcu
1087	!!
1088	NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
1089	& nn_baro, rn_bt_cmax, nn_bt_flt
1090	!!----------------------------------------------------------------------
1091	!
1092	REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters
1093	READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901)
1094	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp )
1095
1096	REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters
1097	READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 )
1098	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp )
1099	IF(lwm) WRITE ( numond, namsplit )
1100	!
1101	! ! Max courant number for ext. grav. waves
1102	!
1103	CALL wrk_alloc( jpi, jpj, zcu )
1104	!
1105	IF (lk_vvl) THEN
1106	DO jj = 1, jpj
1107	DO ji =1, jpi
1108	zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1109	zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1110	zcu(ji,jj) = sqrt(gravht_0(ji,jj)(zxr2 + zyr2) )
1111	END DO
1112	END DO
1113	ELSE
1114	DO jj = 1, jpj
1115	DO ji =1, jpi
1116	zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
1117	zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
1118	zcu(ji,jj) = sqrt(gravht(ji,jj)(zxr2 + zyr2) )
1119	END DO
1120	END DO
1121	ENDIF
1122
1123	zcmax = MAXVAL(zcu(:,:))
1124	IF( lk_mpp ) CALL mpp_max( zcmax )
1125
1126	! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
1127	IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
1128
1129	rdtbt = rdt / FLOAT(nn_baro)
1130	zcmax = zcmax * rdtbt
1131	! Print results
1132	IF(lwp) WRITE(numout,*)
1133	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
1134	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
1135	IF( ln_bt_nn_auto ) THEN
1136	IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro '
1137	IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax
1138	ELSE
1139	IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist '
1140	ENDIF
1141
1142	IF(ln_bt_av) THEN
1143	IF(lwp) WRITE(numout,*) ' ln_bt_av=.true. => Time averaging over nn_baro time steps is on '
1144	ELSE
1145	IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables '
1146	ENDIF
1147	!
1148	!
1149	IF(ln_bt_fw) THEN
1150	IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
1151	ELSE
1152	IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
1153	ENDIF
1154	!
1155	#if defined key_agrif
1156	! Restrict the use of Agrif to the forward case only
1157	IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root())) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
1158	#endif
1159	!
1160	IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt
1161	SELECT CASE ( nn_bt_flt )
1162	CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac'
1163	CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'
1164	CASE( 2 ) ; IF(lwp) WRITE(numout,) ' Boxcar: width = 2nn_baro'
1165	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
1166	END SELECT
1167	!
1168	IF(lwp) WRITE(numout,*) ' '
1169	IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro
1170	IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt
1171	IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax
1172	!
1173	IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN
1174	CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
1175	ENDIF
1176	IF ( zcmax>0.9_wp ) THEN
1177	CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
1178	ENDIF
1179	!
1180	CALL wrk_dealloc( jpi, jpj, zcu )
1181	!
1182	END SUBROUTINE dyn_spg_ts_init
1183
1184	#else
1185	!!---------------------------------------------------------------------------
1186	!! Default case : Empty module No split explicit free surface
1187	!!---------------------------------------------------------------------------
1188	CONTAINS
1189	INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function
1190	dyn_spg_ts_alloc = 0
1191	END FUNCTION dyn_spg_ts_alloc
1192	SUBROUTINE dyn_spg_ts( kt ) ! Empty routine
1193	INTEGER, INTENT(in) :: kt
1194	WRITE(,) 'dyn_spg_ts: You should not have seen this print! error?', kt
1195	END SUBROUTINE dyn_spg_ts
1196	SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine
1197	INTEGER , INTENT(in) :: kt ! ocean time-step
1198	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1199	WRITE(,) 'ts_rst : You should not have seen this print! error?', kt, cdrw
1200	END SUBROUTINE ts_rst
1201	SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine
1202	INTEGER , INTENT(in) :: kt ! ocean time-step
1203	WRITE(,) 'dyn_spg_ts_init : You should not have seen this print! error?', kt
1204	END SUBROUTINE dyn_spg_ts_init
1205	#endif
1206
1207	!!======================================================================
1208	END MODULE dynspg_ts
1209
1210
1211

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

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

Context Navigation

source: branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 9466

Download in other formats: