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

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

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

Context Navigation

source: branches/2014/dev_r5134_UKMO4_CF_compliance/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5350

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

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

Download in other formats: