Context Navigation

dynspg_ts.F90 @ 5666

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

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/2015_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5666

Download in other formats: