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/r5518_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 7099

Last change on this file since 7099 was 7099, checked in by jcastill, 7 years ago
Changes as in r7078 of the original ING branch: branches/2015/dev_r5936_INGV1_WAVE
Property svn:keywords set to `Id`
File size: 55.8 KB

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

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

Download in other formats: