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/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7845

Last change on this file since 7845 was 7845, checked in by glong, 7 years ago
Split dyn_vrt_diag in divcur.F90 into two parts - and changed calls in dynhpg.F90, dynkeg.F90, dynldf_bilap.F90, dynvor.F90, dynzad.F90. Added calls to dyn_vrt_diag in dynspg_ts.F90. Also added call in dynzdf.F90 but this is more of a place holder as it is not currently correct.
File size: 56.3 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 divcur ! for dyn_vrt_dia_2d
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	IMPLICIT NONE
53	PRIVATE
54
55	PUBLIC dyn_spg_ts ! routine called in dynspg.F90
56	PUBLIC dyn_spg_ts_alloc ! " " " "
57	PUBLIC dyn_spg_ts_init ! " " " "
58	PUBLIC ts_rst ! " " " "
59
60	INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
61	REAL(wp),SAVE :: rdtbt ! Barotropic time step
62
63	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: &
64	wgtbtp1, & ! Primary weights used for time filtering of barotropic variables
65	wgtbtp2 ! Secondary weights used for time filtering of barotropic variables
66
67	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff/h at F points
68	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
69	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
70
71	! Arrays below are saved to allow testing of the "no time averaging" option
72	! If this option is not retained, these could be replaced by temporary arrays
73	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, & ! Instantaneous barotropic arrays
74	ubb_e, ub_e, &
75	vbb_e, vb_e
76
77	!! * Substitutions
78	# include "domzgr_substitute.h90"
79	# include "vectopt_loop_substitute.h90"
80	!!----------------------------------------------------------------------
81	!! NEMO/OPA 3.5 , NEMO Consortium (2013)
82	!! $Id$
83	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
84	!!----------------------------------------------------------------------
85	CONTAINS
86
87	INTEGER FUNCTION dyn_spg_ts_alloc()
88	!!----------------------------------------------------------------------
89	!! * routine dyn_spg_ts_alloc *
90	!!----------------------------------------------------------------------
91	INTEGER :: ierr(3)
92	!!----------------------------------------------------------------------
93	ierr(:) = 0
94
95	ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), &
96	& ub_e(jpi,jpj) , vb_e(jpi,jpj) , &
97	& ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , STAT= ierr(1) )
98
99	ALLOCATE( wgtbtp1(3nn_baro), wgtbtp2(3nn_baro), zwz(jpi,jpj), STAT= ierr(2) )
100
101	IF( ln_dynvor_een .or. ln_dynvor_een_old ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &
102	& ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )
103
104	dyn_spg_ts_alloc = MAXVAL(ierr(:))
105
106	IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc )
107	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays')
108	!
109	END FUNCTION dyn_spg_ts_alloc
110
111	SUBROUTINE dyn_spg_ts( kt )
112	!!----------------------------------------------------------------------
113	!!
114	!! ** Purpose :
115	!! -Compute the now trend due to the explicit time stepping
116	!! of the quasi-linear barotropic system.
117	!!
118	!! ** Method :
119	!! Barotropic variables are advanced from internal time steps
120	!! "n" to "n+1" if ln_bt_fw=T
121	!! or from
122	!! "n-1" to "n+1" if ln_bt_fw=F
123	!! thanks to a generalized forward-backward time stepping (see ref. below).
124	!!
125	!! ** Action :
126	!! -Update the filtered free surface at step "n+1" : ssha
127	!! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b
128	!! -Compute barotropic advective velocities at step "n" : un_adv, vn_adv
129	!! These are used to advect tracers and are compliant with discrete
130	!! continuity equation taken at the baroclinic time steps. This
131	!! ensures tracers conservation.
132	!! -Update 3d trend (ua, va) with barotropic component.
133	!!
134	!! References : Shchepetkin, A.F. and J.C. McWilliams, 2005:
135	!! The regional oceanic modeling system (ROMS):
136	!! a split-explicit, free-surface,
137	!! topography-following-coordinate oceanic model.
138	!! Ocean Modelling, 9, 347-404.
139	!!---------------------------------------------------------------------
140	!
141	INTEGER, INTENT(in) :: kt ! ocean time-step index
142	!
143	LOGICAL :: ll_fw_start ! if true, forward integration
144	LOGICAL :: ll_init ! if true, special startup of 2d equations
145	INTEGER :: ji, jj, jk, jn ! dummy loop indices
146	INTEGER :: ikbu, ikbv, noffset ! local integers
147	INTEGER :: id_dia_vrt_spg_int = 1 ! TODO remove once flags set properly
148	INTEGER :: id_dia_vrt_spg_mean = 1 ! TODO remove once flags set properly
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	! !* Local constant initialization
175	z1_12 = 1._wp / 12._wp
176	z1_8 = 0.125_wp
177	z1_4 = 0.25_wp
178	z1_2 = 0.5_wp
179	zraur = 1._wp / rau0
180	!
181	IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step
182	z2dt_bf = rdt
183	ELSE
184	z2dt_bf = 2.0_wp * rdt
185	ENDIF
186	z1_2dt_b = 1.0_wp / z2dt_bf
187	!
188	ll_init = ln_bt_av ! if no time averaging, then no specific restart
189	ll_fw_start = .FALSE.
190	!
191	! time offset in steps for bdy data update
192	IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ; noffset = 0 ; ENDIF
193	!
194	IF( kt == nit000 ) THEN !* initialisation
195	!
196	IF(lwp) WRITE(numout,*)
197	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
198	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
199	IF(lwp) WRITE(numout,*)
200	!
201	IF (neuler==0) ll_init=.TRUE.
202	!
203	IF (ln_bt_fw.OR.(neuler==0)) THEN
204	ll_fw_start=.TRUE.
205	noffset = 0
206	ELSE
207	ll_fw_start=.FALSE.
208	ENDIF
209	!
210	! Set averaging weights and cycle length:
211	CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
212	!
213	!
214	ENDIF
215	!
216	! Set arrays to remove/compute coriolis trend.
217	! Do it once at kt=nit000 if volume is fixed, else at each long time step.
218	! Note that these arrays are also used during barotropic loop. These are however frozen
219	! although they should be updated in the variable volume case. Not a big approximation.
220	! To remove this approximation, copy lines below inside barotropic loop
221	! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
222	!
223	IF ( kt == nit000 .OR. lk_vvl ) THEN
224	IF ( ln_dynvor_een_old ) THEN
225	DO jj = 1, jpjm1
226	DO ji = 1, jpim1
227	zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + &
228	& ht(ji ,jj ) + ht(ji+1,jj ) ) / 4._wp
229	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj)
230	END DO
231	END DO
232	CALL lbc_lnk( zwz, 'F', 1._wp )
233	zwz(:,:) = ff(:,:) * zwz(:,:)
234
235	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
236	DO jj = 2, jpj
237	DO ji = fs_2, jpi ! vector opt.
238	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
239	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
240	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
241	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
242	END DO
243	END DO
244	ELSE IF ( ln_dynvor_een ) THEN
245	DO jj = 1, jpjm1
246	DO ji = 1, jpim1
247	zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + &
248	& ht(ji ,jj ) + ht(ji+1,jj ) ) &
249	& / ( MAX( 1.0_wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + &
250	& tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) )
251	IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj)
252	END DO
253	END DO
254	CALL lbc_lnk( zwz, 'F', 1._wp )
255	zwz(:,:) = ff(:,:) * zwz(:,:)
256
257	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
258	DO jj = 2, jpj
259	DO ji = fs_2, jpi ! vector opt.
260	ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
261	ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
262	ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
263	ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
264	END DO
265	END DO
266	ELSE
267	zwz(:,:) = 0._wp
268	zhf(:,:) = 0.
269	IF ( .not. ln_sco ) THEN
270	! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
271	! Set it to zero for the time being
272	! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level
273	! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
274	! ENDIF
275	! zhf(:,:) = gdepw_0(:,:,jk+1)
276	ELSE
277	zhf(:,:) = hbatf(:,:)
278	END IF
279
280	DO jj = 1, jpjm1
281	zhf(:,jj) = zhf(:,jj)(1._wp- umask(:,jj,1) umask(:,jj+1,1))
282	END DO
283
284	DO jk = 1, jpkm1
285	DO jj = 1, jpjm1
286	zhf(:,jj) = zhf(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
287	END DO
288	END DO
289	CALL lbc_lnk( zhf, 'F', 1._wp )
290	! JC: TBC. hf should be greater than 0
291	DO jj = 1, jpj
292	DO ji = 1, jpi
293	IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
294	END DO
295	END DO
296	zwz(:,:) = ff(:,:) * zwz(:,:)
297	ENDIF
298	ENDIF
299	!
300	! If forward start at previous time step, and centered integration,
301	! then update averaging weights:
302	IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN
303	ll_fw_start=.FALSE.
304	CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
305	ENDIF
306
307	! -----------------------------------------------------------------------------
308	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
309	! -----------------------------------------------------------------------------
310	!
311	!
312	! !* e3*d/dt(Ua) (Vertically integrated)
313	! ! --------------------------------------------------
314	zu_frc(:,:) = 0._wp
315	zv_frc(:,:) = 0._wp
316	!
317	DO jk = 1, jpkm1
318	zu_frc(:,:) = zu_frc(:,:) + fse3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
319	zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
320	END DO
321	!
322	zu_frc(:,:) = zu_frc(:,:) * hur(:,:)
323	zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)
324	!
325	!
326	! !* baroclinic momentum trend (remove the vertical mean trend)
327	DO jk = 1, jpkm1 ! -----------------------------------------------------------
328	DO jj = 2, jpjm1
329	DO ji = fs_2, fs_jpim1 ! vector opt.
330	ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
331	va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
332	END DO
333	END DO
334	END DO
335	! !* barotropic Coriolis trends (vorticity scheme dependent)
336	! ! --------------------------------------------------------
337	zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:) ! now fluxes
338	zwy(:,:) = vn_b(:,:) * hv(:,:) * e1v(:,:)
339	!
340	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
341	DO jj = 2, jpjm1
342	DO ji = fs_2, fs_jpim1 ! vector opt.
343	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
344	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
345	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
346	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
347	! energy conserving formulation for planetary vorticity term
348	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
349	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
350	END DO
351	END DO
352	!
353	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
354	DO jj = 2, jpjm1
355	DO ji = fs_2, fs_jpim1 ! vector opt.
356	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
357	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
358	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
359	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
360	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
361	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
362	END DO
363	END DO
364	!
365	ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN ! enstrophy and energy conserving scheme
366	DO jj = 2, jpjm1
367	DO ji = fs_2, fs_jpim1 ! vector opt.
368	zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
369	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
370	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
371	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
372	zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
373	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
374	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
375	& + ftne(ji,jj ) * zwx(ji ,jj ) )
376	END DO
377	END DO
378	!
379	ENDIF
380	!
381	! !* Right-Hand-Side of the barotropic momentum equation
382	! ! ----------------------------------------------------
383	IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient
384	DO jj = 2, jpjm1
385	DO ji = fs_2, fs_jpim1 ! vector opt.
386	zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj)
387	zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj)
388	END DO
389	END DO
390	ENDIF
391
392	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
393	DO ji = fs_2, fs_jpim1
394	zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)
395	zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
396	END DO
397	END DO
398	!
399	! ! Add bottom stress contribution from baroclinic velocities:
400	IF (ln_bt_fw) THEN
401	DO jj = 2, jpjm1
402	DO ji = fs_2, fs_jpim1 ! vector opt.
403	ikbu = mbku(ji,jj)
404	ikbv = mbkv(ji,jj)
405	zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
406	zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
407	END DO
408	END DO
409	ELSE
410	DO jj = 2, jpjm1
411	DO ji = fs_2, fs_jpim1 ! vector opt.
412	ikbu = mbku(ji,jj)
413	ikbv = mbkv(ji,jj)
414	zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
415	zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
416	END DO
417	END DO
418	ENDIF
419	!
420	! Note that the "unclipped" bottom friction parameter is used even with explicit drag
421	zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
422	zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
423	!
424	IF (ln_bt_fw) THEN ! Add wind forcing
425	zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
426	zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)
427	ELSE
428	zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)
429	zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
430	ENDIF
431	!
432	IF ( ln_apr_dyn ) THEN ! Add atm pressure forcing
433	IF (ln_bt_fw) THEN
434	DO jj = 2, jpjm1
435	DO ji = fs_2, fs_jpim1 ! vector opt.
436	zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
437	zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)
438	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
439	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
440	END DO
441	END DO
442	ELSE
443	DO jj = 2, jpjm1
444	DO ji = fs_2, fs_jpim1 ! vector opt.
445	zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
446	& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj)
447	zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
448	& + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)
449	zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
450	zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
451	END DO
452	END DO
453	ENDIF
454	ENDIF
455	!
456	! Calculate dia_vor_int &
457	! dia_vor_mn if required
458	IF ( ( id_dia_vrt_spg_int == 1 ) .or. (id_dia_vrt_spg_mean == 1 ) ) THEN
459	CALL dyn_vrt_dia_2d(zu_frc, zv_frc, di_dia_vrt_spg_int, id_dia_vrt_spg_mean)
460	END IF
461	! !* Right-Hand-Side of the barotropic ssh equation
462	! ! -----------------------------------------------
463	! ! Surface net water flux and rivers
464	IF (ln_bt_fw) THEN
465	zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
466	ELSE
467	zssh_frc(:,:) = zraur * z1_2 * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) &
468	& + fwfisf(:,:) + fwfisf_b(:,:) )
469	ENDIF
470	#if defined key_asminc
471	! ! Include the IAU weighted SSH increment
472	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
473	zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
474	ENDIF
475	#endif
476	! !* Fill boundary data arrays for AGRIF
477	! ! ------------------------------------
478	#if defined key_agrif
479	IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
480	#endif
481	!
482	! -----------------------------------------------------------------------
483	! Phase 2 : Integration of the barotropic equations
484	! -----------------------------------------------------------------------
485	!
486	! ! ==================== !
487	! ! Initialisations !
488	! ! ==================== !
489	! Initialize barotropic variables:
490	IF( ll_init )THEN
491	sshbb_e(:,:) = 0._wp
492	ubb_e (:,:) = 0._wp
493	vbb_e (:,:) = 0._wp
494	sshb_e (:,:) = 0._wp
495	ub_e (:,:) = 0._wp
496	vb_e (:,:) = 0._wp
497	ENDIF
498	!
499	IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields
500	sshn_e(:,:) = sshn (:,:)
501	zun_e (:,:) = un_b (:,:)
502	zvn_e (:,:) = vn_b (:,:)
503	!
504	hu_e (:,:) = hu (:,:)
505	hv_e (:,:) = hv (:,:)
506	hur_e (:,:) = hur (:,:)
507	hvr_e (:,:) = hvr (:,:)
508	ELSE ! CENTRED integration: start from BEFORE fields
509	sshn_e(:,:) = sshb (:,:)
510	zun_e (:,:) = ub_b (:,:)
511	zvn_e (:,:) = vb_b (:,:)
512	!
513	hu_e (:,:) = hu_b (:,:)
514	hv_e (:,:) = hv_b (:,:)
515	hur_e (:,:) = hur_b(:,:)
516	hvr_e (:,:) = hvr_b(:,:)
517	ENDIF
518	!
519	!
520	!
521	! Initialize sums:
522	ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)
523	va_b (:,:) = 0._wp
524	ssha (:,:) = 0._wp ! Sum for after averaged sea level
525	zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop
526	zv_sum(:,:) = 0._wp
527	! ! ==================== !
528	DO jn = 1, icycle ! sub-time-step loop !
529	! ! ==================== !
530	! !* Update the forcing (BDY and tides)
531	! ! ------------------
532	! Update only tidal forcing at open boundaries
533	#if defined key_tide
534	IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )
535	IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset )
536	#endif
537	!
538	! Set extrapolation coefficients for predictor step:
539	IF ((jn<3).AND.ll_init) THEN ! Forward
540	za1 = 1._wp
541	za2 = 0._wp
542	za3 = 0._wp
543	ELSE ! AB3-AM4 Coefficients: bet=0.281105
544	za1 = 1.781105_wp ! za1 = 3/2 + bet
545	za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet)
546	za3 = 0.281105_wp ! za3 = bet
547	ENDIF
548
549	! Extrapolate barotropic velocities at step jit+0.5:
550	ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
551	va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
552
553	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
554	! ! ------------------
555	! Extrapolate Sea Level at step jit+0.5:
556	zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
557	!
558	DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points
559	DO ji = 2, fs_jpim1 ! Vector opt.
560	zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
561	& * ( e12t(ji ,jj) * zsshp2_e(ji ,jj) &
562	& + e12t(ji+1,jj) * zsshp2_e(ji+1,jj) )
563	zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
564	& * ( e12t(ji,jj ) * zsshp2_e(ji,jj ) &
565	& + e12t(ji,jj+1) * zsshp2_e(ji,jj+1) )
566	END DO
567	END DO
568	CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
569	!
570	zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points
571	zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
572	ELSE
573	zhup2_e (:,:) = hu(:,:)
574	zhvp2_e (:,:) = hv(:,:)
575	ENDIF
576	! !* after ssh
577	! ! -----------
578	! One should enforce volume conservation at open boundaries here
579	! considering fluxes below:
580	!
581	zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
582	zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
583	!
584	#if defined key_agrif
585	! Set fluxes during predictor step to ensure
586	! volume conservation
587	IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN
588	IF((nbondi == -1).OR.(nbondi == 2)) THEN
589	DO jj=1,jpj
590	zwx(2,jj) = ubdy_w(jj) * e2u(2,jj)
591	END DO
592	ENDIF
593	IF((nbondi == 1).OR.(nbondi == 2)) THEN
594	DO jj=1,jpj
595	zwx(nlci-2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)
596	END DO
597	ENDIF
598	IF((nbondj == -1).OR.(nbondj == 2)) THEN
599	DO ji=1,jpi
600	zwy(ji,2) = vbdy_s(ji) * e1v(ji,2)
601	END DO
602	ENDIF
603	IF((nbondj == 1).OR.(nbondj == 2)) THEN
604	DO ji=1,jpi
605	zwy(ji,nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-2)
606	END DO
607	ENDIF
608	ENDIF
609	#endif
610	!
611	! Sum over sub-time-steps to compute advective velocities
612	za2 = wgtbtp2(jn)
613	zu_sum (:,:) = zu_sum (:,:) + za2 * zwx (:,:) / e2u (:,:)
614	zv_sum (:,:) = zv_sum (:,:) + za2 * zwy (:,:) / e1v (:,:)
615	!
616	! Set next sea level:
617	DO jj = 2, jpjm1
618	DO ji = fs_2, fs_jpim1 ! vector opt.
619	zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) &
620	& + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e12t(ji,jj)
621	END DO
622	END DO
623	ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)
624	CALL lbc_lnk( ssha_e, 'T', 1._wp )
625
626	#if defined key_bdy
627	! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.)
628	IF (lk_bdy) CALL bdy_ssh( ssha_e )
629	#endif
630	#if defined key_agrif
631	IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
632	#endif
633	!
634	! Sea Surface Height at u-,v-points (vvl case only)
635	IF ( lk_vvl ) THEN
636	DO jj = 2, jpjm1
637	DO ji = 2, jpim1 ! NO Vector Opt.
638	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
639	& * ( e12t(ji ,jj ) * ssha_e(ji ,jj ) &
640	& + e12t(ji+1,jj ) * ssha_e(ji+1,jj ) )
641	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
642	& * ( e12t(ji ,jj ) * ssha_e(ji ,jj ) &
643	& + e12t(ji ,jj+1) * ssha_e(ji ,jj+1) )
644	END DO
645	END DO
646	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
647	ENDIF
648	!
649	! Half-step back interpolation of SSH for surface pressure computation:
650	!----------------------------------------------------------------------
651	IF ((jn==1).AND.ll_init) THEN
652	za0=1._wp ! Forward-backward
653	za1=0._wp
654	za2=0._wp
655	za3=0._wp
656	ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
657	za0= 1.0833333333333_wp ! za0 = 1-gam-eps
658	za1=-0.1666666666666_wp ! za1 = gam
659	za2= 0.0833333333333_wp ! za2 = eps
660	za3= 0._wp
661	ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
662	za0=0.614_wp ! za0 = 1/2 + gam + 2*eps
663	za1=0.285_wp ! za1 = 1/2 - 2gam - 3eps
664	za2=0.088_wp ! za2 = gam
665	za3=0.013_wp ! za3 = eps
666	ENDIF
667
668	zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) &
669	& + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
670
671	!
672	! Compute associated depths at U and V points:
673	IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN
674	!
675	DO jj = 2, jpjm1
676	DO ji = 2, jpim1
677	zx1 = z1_2 * umask(ji ,jj,1) * r1_e12u(ji ,jj) &
678	& * ( e12t(ji ,jj ) * zsshp2_e(ji ,jj) &
679	& + e12t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )
680	zy1 = z1_2 * vmask(ji ,jj,1) * r1_e12v(ji ,jj ) &
681	& * ( e12t(ji ,jj ) * zsshp2_e(ji ,jj ) &
682	& + e12t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )
683	zhust_e(ji,jj) = hu_0(ji,jj) + zx1
684	zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
685	END DO
686	END DO
687	ENDIF
688	!
689	! Add Coriolis trend:
690	! zwz array below or triads normally depend on sea level with key_vvl and should be updated
691	! at each time step. We however keep them constant here for optimization.
692	! Recall that zwx and zwy arrays hold fluxes at this stage:
693	! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5
694	! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
695	!
696	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
697	DO jj = 2, jpjm1
698	DO ji = fs_2, fs_jpim1 ! vector opt.
699	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
700	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
701	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
702	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
703	zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
704	zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
705	END DO
706	END DO
707	!
708	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
709	DO jj = 2, jpjm1
710	DO ji = fs_2, fs_jpim1 ! vector opt.
711	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
712	& + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
713	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
714	& + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
715	zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )
716	zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )
717	END DO
718	END DO
719	!
720	ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !== energy and enstrophy conserving scheme ==!
721	DO jj = 2, jpjm1
722	DO ji = fs_2, fs_jpim1 ! vector opt.
723	zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &
724	& + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
725	& + ftse(ji,jj ) * zwy(ji ,jj-1) &
726	& + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
727	zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
728	& + ftse(ji,jj+1) * zwx(ji ,jj+1) &
729	& + ftnw(ji,jj ) * zwx(ji-1,jj ) &
730	& + ftne(ji,jj ) * zwx(ji ,jj ) )
731	END DO
732	END DO
733	!
734	ENDIF
735	!
736	! Add tidal astronomical forcing if defined
737	IF ( lk_tide.AND.ln_tide_pot ) THEN
738	DO jj = 2, jpjm1
739	DO ji = fs_2, fs_jpim1 ! vector opt.
740	zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
741	zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
742	zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
743	zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
744	END DO
745	END DO
746	ENDIF
747	!
748	! Add bottom stresses:
749	zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)
750	zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)
751	!
752	! Surface pressure trend:
753	DO jj = 2, jpjm1
754	DO ji = fs_2, fs_jpim1 ! vector opt.
755	! Add surface pressure gradient
756	zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
757	zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
758	zwx(ji,jj) = zu_spg
759	zwy(ji,jj) = zv_spg
760	END DO
761	END DO
762	!
763	! Set next velocities:
764	IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN ! Vector form
765	DO jj = 2, jpjm1
766	DO ji = fs_2, fs_jpim1 ! vector opt.
767	ua_e(ji,jj) = ( zun_e(ji,jj) &
768	& + rdtbt * ( zwx(ji,jj) &
769	& + zu_trd(ji,jj) &
770	& + zu_frc(ji,jj) ) &
771	& ) * umask(ji,jj,1)
772
773	va_e(ji,jj) = ( zvn_e(ji,jj) &
774	& + rdtbt * ( zwy(ji,jj) &
775	& + zv_trd(ji,jj) &
776	& + zv_frc(ji,jj) ) &
777	& ) * vmask(ji,jj,1)
778	END DO
779	END DO
780
781	ELSE ! Flux form
782	DO jj = 2, jpjm1
783	DO ji = fs_2, fs_jpim1 ! vector opt.
784
785	zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
786	zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
787
788	ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) &
789	& + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &
790	& + zhup2_e(ji,jj) * zu_trd(ji,jj) &
791	& + hu(ji,jj) * zu_frc(ji,jj) ) &
792	& ) * zhura
793
794	va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &
795	& + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) &
796	& + zhvp2_e(ji,jj) * zv_trd(ji,jj) &
797	& + hv(ji,jj) * zv_frc(ji,jj) ) &
798	& ) * zhvra
799	END DO
800	END DO
801	ENDIF
802	!
803	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
804	! ! ----------------------------------------------
805	hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
806	hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
807	hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
808	hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
809	!
810	ENDIF
811	! !* domain lateral boundary
812	! ! -----------------------
813	!
814	CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp )
815
816	#if defined key_bdy
817	! open boundaries
818	IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )
819	#endif
820	#if defined key_agrif
821	IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif
822	#endif
823	! !* Swap
824	! ! ----
825	ubb_e (:,:) = ub_e (:,:)
826	ub_e (:,:) = zun_e (:,:)
827	zun_e (:,:) = ua_e (:,:)
828	!
829	vbb_e (:,:) = vb_e (:,:)
830	vb_e (:,:) = zvn_e (:,:)
831	zvn_e (:,:) = va_e (:,:)
832	!
833	sshbb_e(:,:) = sshb_e(:,:)
834	sshb_e (:,:) = sshn_e(:,:)
835	sshn_e (:,:) = ssha_e(:,:)
836
837	! !* Sum over whole bt loop
838	! ! ----------------------
839	za1 = wgtbtp1(jn)
840	IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities
841	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)
842	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)
843	ELSE ! Sum transports
844	ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)
845	va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)
846	ENDIF
847	! ! Sum sea level
848	ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
849	! ! ==================== !
850	END DO ! end loop !
851	! ! ==================== !
852	! -----------------------------------------------------------------------------
853	! Phase 3. update the general trend with the barotropic trend
854	! -----------------------------------------------------------------------------
855	!
856	! At this stage ssha holds a time averaged value
857	! ! Sea Surface Height at u-,v- and f-points
858	IF( lk_vvl ) THEN ! (required only in key_vvl case)
859	DO jj = 1, jpjm1
860	DO ji = 1, jpim1 ! NO Vector Opt.
861	zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &
862	& * ( e12t(ji ,jj) * ssha(ji ,jj) &
863	& + e12t(ji+1,jj) * ssha(ji+1,jj) )
864	zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &
865	& * ( e12t(ji,jj ) * ssha(ji,jj ) &
866	& + e12t(ji,jj+1) * ssha(ji,jj+1) )
867	END DO
868	END DO
869	CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
870	ENDIF
871	!
872	! Set advection velocity correction:
873	IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN
874	un_adv(:,:) = zu_sum(:,:)*hur(:,:)
875	vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
876	ELSE
877	un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:)
878	vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:)
879	END IF
880
881	IF (ln_bt_fw) THEN ! Save integrated transport for next computation
882	ub2_b(:,:) = zu_sum(:,:)
883	vb2_b(:,:) = zv_sum(:,:)
884	ENDIF
885	!
886	! Update barotropic trend:
887	IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN
888	DO jk=1,jpkm1
889	ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
890	va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
891	END DO
892	ELSE
893	DO jk=1,jpkm1
894	ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
895	va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
896	END DO
897	! Save barotropic velocities not transport:
898	ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )
899	va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )
900	ENDIF
901	!
902	DO jk = 1, jpkm1
903	! Correct velocities:
904	un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk)
905	vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk)
906	!
907	END DO
908	!
909	#if defined key_agrif
910	! Save time integrated fluxes during child grid integration
911	! (used to update coarse grid transports at next time step)
912	!
913	IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN
914	IF ( Agrif_NbStepint().EQ.0 ) THEN
915	ub2_i_b(:,:) = 0.e0
916	vb2_i_b(:,:) = 0.e0
917	END IF
918	!
919	za1 = 1._wp / REAL(Agrif_rhot(), wp)
920	ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
921	vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
922	ENDIF
923	!
924	!
925	#endif
926	!
927	! !* write time-spliting arrays in the restart
928	IF(lrst_oce .AND.ln_bt_fw) CALL ts_rst( kt, 'WRITE' )
929	!
930	CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv )
931	CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e )
932	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )
933	CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e )
934	CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a )
935	CALL wrk_dealloc( jpi, jpj, zhf )
936	!
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.

Download in other formats: