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

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

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

Context Navigation

source: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 5208

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

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

Download in other formats: