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 @ 7900

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

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

Download in other formats: