Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynspg_ts.F90 @ 8243

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

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

Download in other formats: