Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90 @ 12359

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

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/dev_r5518_obs_oper_update_PS44/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 12359

Download in other formats: