Home

Context Navigation

source: palm/trunk/SOURCE/advec_ws.f90 @ 2341

Last change on this file since 2341 was 2329, checked in by knoop, 8 years ago
Bugfix for topography usage with anelastic approximation and boussinesq approximation with air density != 1
Property svn:keywords set to `Id`
File size: 278.1 KB

Line
1	!> @file advec_ws.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: advec_ws.f90 2329 2017-08-03 14:24:56Z Giersch $
27	! Bugfix concerning density in divergence correction close to buildings
28	!
29	! 2292 2017-06-20 09:51:42Z schwenkel
30	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
31	! includes two more prognostic equations for cloud drop concentration (nc)
32	! and cloud water content (qc).
33	!
34	! 2233 2017-05-30 18:08:54Z suehring
35	!
36	! 2232 2017-05-30 17:47:52Z suehring
37	! Rename wall_flags_0 and wall_flags_00 into advc_flags_1 and advc_flags_2,
38	! respectively.
39	! Set advc_flags_1/2 on basis of wall_flags_0/00 instead of nzb_s/u/v/w_inner.
40	! Setting advc_flags_1/2 also for downward-facing walls
41	!
42	! 2200 2017-04-11 11:37:51Z suehring
43	! monotonic_adjustment removed
44	!
45	! 2118 2017-01-17 16:38:49Z raasch
46	! OpenACC version of subroutines removed
47	!
48	! 2037 2016-10-26 11:15:40Z knoop
49	! Anelastic approximation implemented
50	!
51	! 2000 2016-08-20 18:09:15Z knoop
52	! Forced header and separation lines into 80 columns
53	!
54	! 1996 2016-08-18 11:42:29Z suehring
55	! Bugfix concerning calculation of turbulent of turbulent fluxes
56	!
57	! 1960 2016-07-12 16:34:24Z suehring
58	! Separate humidity and passive scalar
59	!
60	! 1942 2016-06-14 12:18:18Z suehring
61	! Initialization of flags for ws-scheme moved from init_grid.
62	!
63	! 1873 2016-04-18 14:50:06Z maronga
64	! Module renamed (removed _mod)
65	!
66	!
67	! 1850 2016-04-08 13:29:27Z maronga
68	! Module renamed
69	!
70	!
71	! 1822 2016-04-07 07:49:42Z hoffmann
72	! icloud_scheme removed, microphysics_seifert added
73	!
74	! 1682 2015-10-07 23:56:08Z knoop
75	! Code annotations made doxygen readable
76	!
77	! 1630 2015-08-26 16:57:23Z suehring
78	!
79	!
80	! 1629 2015-08-26 16:56:11Z suehring
81	! Bugfix concerning wall_flags at left and south PE boundaries
82	!
83	! 1581 2015-04-10 13:45:59Z suehring
84	!
85	!
86	! 1580 2015-04-10 13:43:49Z suehring
87	! Bugfix: statistical evaluation of scalar fluxes in case of monotonic limiter
88	!
89	! 1567 2015-03-10 17:57:55Z suehring
90	! Bugfixes in monotonic limiter.
91	!
92	! 2015-03-09 13:10:37Z heinze
93	! Bugfix: REAL constants provided with KIND-attribute in call of
94	! intrinsic functions like MAX and MIN
95	!
96	! 1557 2015-03-05 16:43:04Z suehring
97	! Enable monotone advection for scalars using monotonic limiter
98	!
99	! 1374 2014-04-25 12:55:07Z raasch
100	! missing variables added to ONLY list
101	!
102	! 1361 2014-04-16 15:17:48Z hoffmann
103	! accelerator and vector version for qr and nr added
104	!
105	! 1353 2014-04-08 15:21:23Z heinze
106	! REAL constants provided with KIND-attribute,
107	! module kinds added
108	! some formatting adjustments
109	!
110	! 1322 2014-03-20 16:38:49Z raasch
111	! REAL constants defined as wp-kind
112	!
113	! 1320 2014-03-20 08:40:49Z raasch
114	! ONLY-attribute added to USE-statements,
115	! kind-parameters added to all INTEGER and REAL declaration statements,
116	! kinds are defined in new module kinds,
117	! old module precision_kind is removed,
118	! revision history before 2012 removed,
119	! comment fields (!:) to be used for variable explanations added to
120	! all variable declaration statements
121	!
122	! 1257 2013-11-08 15:18:40Z raasch
123	! accelerator loop directives removed
124	!
125	! 1221 2013-09-10 08:59:13Z raasch
126	! wall_flags_00 introduced, which holds bits 32-...
127	!
128	! 1128 2013-04-12 06:19:32Z raasch
129	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
130	! j_north
131	!
132	! 1115 2013-03-26 18:16:16Z hoffmann
133	! calculation of qr and nr is restricted to precipitation
134	!
135	! 1053 2012-11-13 17:11:03Z hoffmann
136	! necessary expansions according to the two new prognostic equations (nr, qr)
137	! of the two-moment cloud physics scheme:
138	! +flux_l_, flux_s_, diss_l_, diss_s_, sums_ws*s_ws_l
139	!
140	! 1036 2012-10-22 13:43:42Z raasch
141	! code put under GPL (PALM 3.9)
142	!
143	! 1027 2012-10-15 17:18:39Z suehring
144	! Bugfix in calculation indices k_mm, k_pp in accelerator version
145	!
146	! 1019 2012-09-28 06:46:45Z raasch
147	! small change in comment lines
148	!
149	! 1015 2012-09-27 09:23:24Z raasch
150	! accelerator versions (*_acc) added
151	!
152	! 1010 2012-09-20 07:59:54Z raasch
153	! cpp switch __nopointer added for pointer free version
154	!
155	! 888 2012-04-20 15:03:46Z suehring
156	! Number of IBITS() calls with identical arguments is reduced.
157	!
158	! 862 2012-03-26 14:21:38Z suehring
159	! ws-scheme also work with topography in combination with vector version.
160	! ws-scheme also work with outflow boundaries in combination with
161	! vector version.
162	! Degradation of the applied order of scheme is now steered by multiplying with
163	! Integer advc_flags_1. 2nd order scheme, WS3 and WS5 are calculated on each
164	! grid point and mulitplied with the appropriate flag.
165	! 2nd order numerical dissipation term changed. Now the appropriate 2nd order
166	! term derived according to the 4th and 6th order terms is applied. It turns
167	! out that diss_2nd does not provide sufficient dissipation near walls.
168	! Therefore, the function diss_2nd is removed.
169	! Near walls a divergence correction is necessary to overcome numerical
170	! instabilities due to too less divergence reduction of the velocity field.
171	! boundary_flags and logicals steering the degradation are removed.
172	! Empty SUBROUTINE local_diss removed.
173	! Further formatting adjustments.
174	!
175	! 801 2012-01-10 17:30:36Z suehring
176	! Bugfix concerning OpenMP parallelization. Summation of sums_wsus_ws_l,
177	! sums_wsvs_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wspts_ws_l,
178	! sums_wsqs_ws_l, sums_wssas_ws_l is now thread-safe by adding an additional
179	! dimension.
180	!
181	! Initial revision
182	!
183	! 411 2009-12-11 12:31:43 Z suehring
184	!
185	! Description:
186	! ------------
187	!> Advection scheme for scalars and momentum using the flux formulation of
188	!> Wicker and Skamarock 5th order. Additionally the module contains of a
189	!> routine using for initialisation and steering of the statical evaluation.
190	!> The computation of turbulent fluxes takes place inside the advection
191	!> routines.
192	!> Near non-cyclic boundaries the order of the applied advection scheme is
193	!> degraded.
194	!> A divergence correction is applied. It is necessary for topography, since
195	!> the divergence is not sufficiently reduced, resulting in erroneous fluxes and
196	!> partly numerical instabilities.
197	!-----------------------------------------------------------------------------!
198	MODULE advec_ws
199
200
201
202	PRIVATE
203	PUBLIC advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, ws_init, &
204	ws_init_flags, ws_statistics
205
206	INTERFACE ws_init
207	MODULE PROCEDURE ws_init
208	END INTERFACE ws_init
209
210	INTERFACE ws_init_flags
211	MODULE PROCEDURE ws_init_flags
212	END INTERFACE ws_init_flags
213
214	INTERFACE ws_statistics
215	MODULE PROCEDURE ws_statistics
216	END INTERFACE ws_statistics
217
218	INTERFACE advec_s_ws
219	MODULE PROCEDURE advec_s_ws
220	MODULE PROCEDURE advec_s_ws_ij
221	END INTERFACE advec_s_ws
222
223	INTERFACE advec_u_ws
224	MODULE PROCEDURE advec_u_ws
225	MODULE PROCEDURE advec_u_ws_ij
226	END INTERFACE advec_u_ws
227
228	INTERFACE advec_v_ws
229	MODULE PROCEDURE advec_v_ws
230	MODULE PROCEDURE advec_v_ws_ij
231	END INTERFACE advec_v_ws
232
233	INTERFACE advec_w_ws
234	MODULE PROCEDURE advec_w_ws
235	MODULE PROCEDURE advec_w_ws_ij
236	END INTERFACE advec_w_ws
237
238	CONTAINS
239
240
241	!------------------------------------------------------------------------------!
242	! Description:
243	! ------------
244	!> Initialization of WS-scheme
245	!------------------------------------------------------------------------------!
246	SUBROUTINE ws_init
247
248	USE arrays_3d, &
249	ONLY: diss_l_e, diss_l_nc, diss_l_nr, diss_l_pt, diss_l_q, &
250	diss_l_qc, diss_l_qr, diss_l_s, diss_l_sa, diss_l_u, &
251	diss_l_v, diss_l_w, flux_l_e, flux_l_nc, flux_l_nr, &
252	flux_l_pt, flux_l_q, flux_l_qc, flux_l_qr, flux_l_s, &
253	flux_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_s_e, &
254	diss_s_nc, diss_s_nr, diss_s_pt, diss_s_q, diss_s_qc, &
255	diss_s_qr, diss_s_s, diss_s_sa, diss_s_u, diss_s_v, &
256	diss_s_w, flux_s_e, flux_s_nc, flux_s_nr, flux_s_pt, &
257	flux_s_q, flux_s_qc, flux_s_qr, flux_s_s, flux_s_sa, &
258	flux_s_u, flux_s_v, flux_s_w
259
260	USE constants, &
261	ONLY: adv_mom_1, adv_mom_3, adv_mom_5, adv_sca_1, adv_sca_3, &
262	adv_sca_5
263
264	USE control_parameters, &
265	ONLY: cloud_physics, humidity, loop_optimization, &
266	passive_scalar, microphysics_morrison, microphysics_seifert, &
267	ocean, ws_scheme_mom, ws_scheme_sca
268
269	USE indices, &
270	ONLY: nyn, nys, nzb, nzt
271
272	USE kinds
273
274	USE pegrid
275
276	USE statistics, &
277	ONLY: sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsncs_ws_l,&
278	sums_wsnrs_ws_l,sums_wspts_ws_l, sums_wsqcs_ws_l, &
279	sums_wsqrs_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l, &
280	sums_wssas_ws_l, sums_wsss_ws_l, sums_wsus_ws_l, &
281	sums_wsvs_ws_l
282
283
284	!
285	!-- Set the appropriate factors for scalar and momentum advection.
286	adv_sca_5 = 1.0_wp / 60.0_wp
287	adv_sca_3 = 1.0_wp / 12.0_wp
288	adv_sca_1 = 1.0_wp / 2.0_wp
289	adv_mom_5 = 1.0_wp / 120.0_wp
290	adv_mom_3 = 1.0_wp / 24.0_wp
291	adv_mom_1 = 1.0_wp / 4.0_wp
292	!
293	!-- Arrays needed for statical evaluation of fluxes.
294	IF ( ws_scheme_mom ) THEN
295
296	ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1), &
297	sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1), &
298	sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
299	sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
300	sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
301
302	sums_wsus_ws_l = 0.0_wp
303	sums_wsvs_ws_l = 0.0_wp
304	sums_us2_ws_l = 0.0_wp
305	sums_vs2_ws_l = 0.0_wp
306	sums_ws2_ws_l = 0.0_wp
307
308	ENDIF
309
310	IF ( ws_scheme_sca ) THEN
311
312	ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
313	sums_wspts_ws_l = 0.0_wp
314
315	IF ( humidity ) THEN
316	ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
317	sums_wsqs_ws_l = 0.0_wp
318	ENDIF
319
320	IF ( passive_scalar ) THEN
321	ALLOCATE( sums_wsss_ws_l(nzb:nzt+1,0:threads_per_task-1) )
322	sums_wsss_ws_l = 0.0_wp
323	ENDIF
324
325	IF ( cloud_physics .AND. microphysics_morrison ) THEN
326	ALLOCATE( sums_wsqcs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
327	ALLOCATE( sums_wsncs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
328	sums_wsqcs_ws_l = 0.0_wp
329	sums_wsncs_ws_l = 0.0_wp
330	ENDIF
331
332	IF ( cloud_physics .AND. microphysics_seifert ) THEN
333	ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
334	ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
335	sums_wsqrs_ws_l = 0.0_wp
336	sums_wsnrs_ws_l = 0.0_wp
337	ENDIF
338
339	IF ( ocean ) THEN
340	ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
341	sums_wssas_ws_l = 0.0_wp
342	ENDIF
343
344	ENDIF
345
346	!
347	!-- Arrays needed for reasons of speed optimization for cache version.
348	!-- For the vector version the buffer arrays are not necessary,
349	!-- because the the fluxes can swapped directly inside the loops of the
350	!-- advection routines.
351	IF ( loop_optimization /= 'vector' ) THEN
352
353	IF ( ws_scheme_mom ) THEN
354
355	ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1), &
356	flux_s_v(nzb+1:nzt,0:threads_per_task-1), &
357	flux_s_w(nzb+1:nzt,0:threads_per_task-1), &
358	diss_s_u(nzb+1:nzt,0:threads_per_task-1), &
359	diss_s_v(nzb+1:nzt,0:threads_per_task-1), &
360	diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
361	ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
362	flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
363	flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
364	diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
365	diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
366	diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
367
368	ENDIF
369
370	IF ( ws_scheme_sca ) THEN
371
372	ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1), &
373	flux_s_e(nzb+1:nzt,0:threads_per_task-1), &
374	diss_s_pt(nzb+1:nzt,0:threads_per_task-1), &
375	diss_s_e(nzb+1:nzt,0:threads_per_task-1) )
376	ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
377	flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
378	diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
379	diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
380
381	IF ( humidity ) THEN
382	ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1), &
383	diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
384	ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
385	diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
386	ENDIF
387
388	IF ( passive_scalar ) THEN
389	ALLOCATE( flux_s_s(nzb+1:nzt,0:threads_per_task-1), &
390	diss_s_s(nzb+1:nzt,0:threads_per_task-1) )
391	ALLOCATE( flux_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
392	diss_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
393	ENDIF
394
395	IF ( cloud_physics .AND. microphysics_morrison ) THEN
396	ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1), &
397	diss_s_qc(nzb+1:nzt,0:threads_per_task-1), &
398	flux_s_nc(nzb+1:nzt,0:threads_per_task-1), &
399	diss_s_nc(nzb+1:nzt,0:threads_per_task-1) )
400	ALLOCATE( flux_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
401	diss_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
402	flux_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
403	diss_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
404	ENDIF
405
406	IF ( cloud_physics .AND. microphysics_seifert ) THEN
407	ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1), &
408	diss_s_qr(nzb+1:nzt,0:threads_per_task-1), &
409	flux_s_nr(nzb+1:nzt,0:threads_per_task-1), &
410	diss_s_nr(nzb+1:nzt,0:threads_per_task-1) )
411	ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
412	diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
413	flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
414	diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
415	ENDIF
416
417	IF ( ocean ) THEN
418	ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1), &
419	diss_s_sa(nzb+1:nzt,0:threads_per_task-1) )
420	ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
421	diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
422	ENDIF
423
424	ENDIF
425
426	ENDIF
427
428	END SUBROUTINE ws_init
429
430	!------------------------------------------------------------------------------!
431	! Description:
432	! ------------
433	!> Initialization of flags for WS-scheme used to degrade the order of the scheme
434	!> near walls.
435	!------------------------------------------------------------------------------!
436	SUBROUTINE ws_init_flags
437
438	USE control_parameters, &
439	ONLY: inflow_l, inflow_n, inflow_r, inflow_s, momentum_advec, &
440	nest_bound_l, nest_bound_n, nest_bound_r, nest_bound_s, &
441	outflow_l, outflow_n, outflow_r, outflow_s, scalar_advec
442
443	USE indices, &
444	ONLY: advc_flags_1, advc_flags_2, nbgp, nxl, nxlu, nxr, nyn, nys, &
445	nysv, nzb, nzt, wall_flags_0
446
447	USE kinds
448
449	IMPLICIT NONE
450
451	INTEGER(iwp) :: i !< index variable along x
452	INTEGER(iwp) :: j !< index variable along y
453	INTEGER(iwp) :: k !< index variable along z
454	INTEGER(iwp) :: k_mm !< dummy index along z
455	INTEGER(iwp) :: k_pp !< dummy index along z
456	INTEGER(iwp) :: k_ppp !< dummy index along z
457
458	LOGICAL :: flag_set !< steering variable for advection flags
459
460
461	IF ( scalar_advec == 'ws-scheme' ) THEN
462	!
463	!-- Set flags to steer the degradation of the advection scheme in advec_ws
464	!-- near topography, inflow- and outflow boundaries as well as bottom and
465	!-- top of model domain. advc_flags_1 remains zero for all non-prognostic
466	!-- grid points.
467	DO i = nxl, nxr
468	DO j = nys, nyn
469	DO k = nzb+1, nzt
470	!
471	!-- scalar - x-direction
472	!-- WS1 (0), WS3 (1), WS5 (2)
473	IF ( ( .NOT. BTEST(wall_flags_0(k,j,i+1),0) &
474	.OR. .NOT. BTEST(wall_flags_0(k,j,i+2),0) &
475	.OR. .NOT. BTEST(wall_flags_0(k,j,i-1),0) ) &
476	.OR. ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
477	.AND. i == nxl ) &
478	.OR. ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
479	.AND. i == nxr ) ) &
480	THEN
481	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 0 )
482	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+3),0) &
483	.AND. BTEST(wall_flags_0(k,j,i+1),0) &
484	.AND. BTEST(wall_flags_0(k,j,i+2),0) &
485	.AND. BTEST(wall_flags_0(k,j,i-1),0) &
486	) .OR. &
487	( .NOT. BTEST(wall_flags_0(k,j,i-2),0) &
488	.AND. BTEST(wall_flags_0(k,j,i+1),0) &
489	.AND. BTEST(wall_flags_0(k,j,i+2),0) &
490	.AND. BTEST(wall_flags_0(k,j,i-1),0) &
491	) &
492	.OR. &
493	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
494	.AND. i == nxr-1 ) .OR. &
495	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
496	.AND. i == nxlu ) ) &
497	THEN
498	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 1 )
499	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),0) &
500	.AND. BTEST(wall_flags_0(k,j,i+2),0) &
501	.AND. BTEST(wall_flags_0(k,j,i+3),0) &
502	.AND. BTEST(wall_flags_0(k,j,i-1),0) &
503	.AND. BTEST(wall_flags_0(k,j,i-2),0) ) &
504	THEN
505	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 2 )
506	ENDIF
507	!
508	!-- scalar - y-direction
509	!-- WS1 (3), WS3 (4), WS5 (5)
510	IF ( ( .NOT. BTEST(wall_flags_0(k,j+1,i),0) &
511	.OR. .NOT. BTEST(wall_flags_0(k,j+2,i),0) &
512	.OR. .NOT. BTEST(wall_flags_0(k,j-1,i),0)) &
513	.OR. ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
514	.AND. j == nys ) &
515	.OR. ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
516	.AND. j == nyn ) ) &
517	THEN
518	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 3 )
519	!
520	!-- WS3
521	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+3,i),0) &
522	.AND. BTEST(wall_flags_0(k,j+1,i),0) &
523	.AND. BTEST(wall_flags_0(k,j+2,i),0) &
524	.AND. BTEST(wall_flags_0(k,j-1,i),0) &
525	) .OR. &
526	( .NOT. BTEST(wall_flags_0(k,j-2,i),0) &
527	.AND. BTEST(wall_flags_0(k,j+1,i),0) &
528	.AND. BTEST(wall_flags_0(k,j+2,i),0) &
529	.AND. BTEST(wall_flags_0(k,j-1,i),0) &
530	) &
531	.OR. &
532	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
533	.AND. j == nysv ) .OR. &
534	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
535	.AND. j == nyn-1 ) ) &
536	THEN
537	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 4 )
538	!
539	!-- WS5
540	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),0) &
541	.AND. BTEST(wall_flags_0(k,j+2,i),0) &
542	.AND. BTEST(wall_flags_0(k,j+3,i),0) &
543	.AND. BTEST(wall_flags_0(k,j-1,i),0) &
544	.AND. BTEST(wall_flags_0(k,j-2,i),0) ) &
545	THEN
546	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 5 )
547	ENDIF
548	!
549	!-- scalar - z-direction
550	!-- WS1 (6), WS3 (7), WS5 (8)
551	IF ( k == nzb+1 ) THEN
552	k_mm = nzb
553	ELSE
554	k_mm = k - 2
555	ENDIF
556	IF ( k > nzt-1 ) THEN
557	k_pp = nzt+1
558	ELSE
559	k_pp = k + 2
560	ENDIF
561	IF ( k > nzt-2 ) THEN
562	k_ppp = nzt+1
563	ELSE
564	k_ppp = k + 3
565	ENDIF
566
567	flag_set = .FALSE.
568	IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),0) .AND. &
569	BTEST(wall_flags_0(k,j,i),0) .OR. &
570	.NOT. BTEST(wall_flags_0(k_pp,j,i),0) .AND. &
571	BTEST(wall_flags_0(k,j,i),0) .OR. &
572	k == nzt ) &
573	THEN
574	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 6 )
575	flag_set = .TRUE.
576	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),0) .OR. &
577	.NOT. BTEST(wall_flags_0(k_ppp,j,i),0) ) .AND. &
578	BTEST(wall_flags_0(k-1,j,i),0) .AND. &
579	BTEST(wall_flags_0(k,j,i),0) .AND. &
580	BTEST(wall_flags_0(k+1,j,i),0) .AND. &
581	BTEST(wall_flags_0(k_pp,j,i),0) .OR. &
582	k == nzt - 1 ) &
583	THEN
584	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 7 )
585	flag_set = .TRUE.
586	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),0) &
587	.AND. BTEST(wall_flags_0(k-1,j,i),0) &
588	.AND. BTEST(wall_flags_0(k,j,i),0) &
589	.AND. BTEST(wall_flags_0(k+1,j,i),0) &
590	.AND. BTEST(wall_flags_0(k_pp,j,i),0) &
591	.AND. BTEST(wall_flags_0(k_ppp,j,i),0) &
592	.AND. .NOT. flag_set ) &
593	THEN
594	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 8 )
595	ENDIF
596
597	ENDDO
598	ENDDO
599	ENDDO
600	ENDIF
601
602	IF ( momentum_advec == 'ws-scheme' ) THEN
603	!
604	!-- Set advc_flags_1 to steer the degradation of the advection scheme in advec_ws
605	!-- near topography, inflow- and outflow boundaries as well as bottom and
606	!-- top of model domain. advc_flags_1 remains zero for all non-prognostic
607	!-- grid points.
608	DO i = nxl, nxr
609	DO j = nys, nyn
610	DO k = nzb+1, nzt
611	!
612	!-- At first, set flags to WS1.
613	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
614	!-- in order to handle the left/south flux.
615	!-- near vertical walls.
616	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 9 )
617	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 12 )
618	!
619	!-- u component - x-direction
620	!-- WS1 (9), WS3 (10), WS5 (11)
621	IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),1) .OR. &
622	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
623	.AND. i <= nxlu ) .OR. &
624	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
625	.AND. i == nxr ) ) &
626	THEN
627	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 9 )
628	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),1) .AND. &
629	BTEST(wall_flags_0(k,j,i+1),1) .OR. &
630	.NOT. BTEST(wall_flags_0(k,j,i-1),1) ) &
631	.OR. &
632	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
633	.AND. i == nxr-1 ) .OR. &
634	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
635	.AND. i == nxlu+1) ) &
636	THEN
637	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 10 )
638	!
639	!-- Clear flag for WS1
640	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 9 )
641	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),1) .AND. &
642	BTEST(wall_flags_0(k,j,i+2),1) .AND. &
643	BTEST(wall_flags_0(k,j,i-1),1) ) &
644	THEN
645	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 11 )
646	!
647	!-- Clear flag for WS1
648	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 9 )
649	ENDIF
650	!
651	!-- u component - y-direction
652	!-- WS1 (12), WS3 (13), WS5 (14)
653	IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),1) .OR. &
654	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
655	.AND. j == nys ) .OR. &
656	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
657	.AND. j == nyn ) ) &
658	THEN
659	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 12 )
660	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),1) .AND. &
661	BTEST(wall_flags_0(k,j+1,i),1) .OR. &
662	.NOT. BTEST(wall_flags_0(k,j-1,i),1) ) &
663	.OR. &
664	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
665	.AND. j == nysv ) .OR. &
666	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
667	.AND. j == nyn-1 ) ) &
668	THEN
669	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 13 )
670	!
671	!-- Clear flag for WS1
672	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 12 )
673	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),1) .AND. &
674	BTEST(wall_flags_0(k,j+2,i),1) .AND. &
675	BTEST(wall_flags_0(k,j-1,i),1) ) &
676	THEN
677	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 14 )
678	!
679	!-- Clear flag for WS1
680	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 12 )
681	ENDIF
682	!
683	!-- u component - z-direction
684	!-- WS1 (15), WS3 (16), WS5 (17)
685	IF ( k == nzb+1 ) THEN
686	k_mm = nzb
687	ELSE
688	k_mm = k - 2
689	ENDIF
690	IF ( k > nzt-1 ) THEN
691	k_pp = nzt+1
692	ELSE
693	k_pp = k + 2
694	ENDIF
695	IF ( k > nzt-2 ) THEN
696	k_ppp = nzt+1
697	ELSE
698	k_ppp = k + 3
699	ENDIF
700
701	flag_set = .FALSE.
702	IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),1) .AND. &
703	BTEST(wall_flags_0(k,j,i),1) .OR. &
704	.NOT. BTEST(wall_flags_0(k_pp,j,i),1) .AND. &
705	BTEST(wall_flags_0(k,j,i),1) .OR. &
706	k == nzt ) &
707	THEN
708	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 15 )
709	flag_set = .TRUE.
710	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),1) .OR. &
711	.NOT. BTEST(wall_flags_0(k_ppp,j,i),1) ) .AND. &
712	BTEST(wall_flags_0(k-1,j,i),1) .AND. &
713	BTEST(wall_flags_0(k,j,i),1) .AND. &
714	BTEST(wall_flags_0(k+1,j,i),1) .AND. &
715	BTEST(wall_flags_0(k_pp,j,i),1) .OR. &
716	k == nzt - 1 ) &
717	THEN
718	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 16 )
719	flag_set = .TRUE.
720	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),1) .AND. &
721	BTEST(wall_flags_0(k-1,j,i),1) .AND. &
722	BTEST(wall_flags_0(k,j,i),1) .AND. &
723	BTEST(wall_flags_0(k+1,j,i),1) .AND. &
724	BTEST(wall_flags_0(k_pp,j,i),1) .AND. &
725	BTEST(wall_flags_0(k_ppp,j,i),1) .AND. &
726	.NOT. flag_set ) &
727	THEN
728	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 17 )
729	ENDIF
730
731	ENDDO
732	ENDDO
733	ENDDO
734
735	DO i = nxl, nxr
736	DO j = nys, nyn
737	DO k = nzb+1, nzt
738	!
739	!-- At first, set flags to WS1.
740	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
741	!-- in order to handle the left/south flux.
742	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 18 )
743	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 21 )
744	!
745	!-- v component - x-direction
746	!-- WS1 (18), WS3 (19), WS5 (20)
747	IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),2) .OR. &
748	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
749	.AND. i == nxl ) .OR. &
750	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
751	.AND. i == nxr ) ) &
752	THEN
753	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 18 )
754	!
755	!-- WS3
756	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),2) .AND. &
757	BTEST(wall_flags_0(k,j,i+1),2) ) .OR. &
758	.NOT. BTEST(wall_flags_0(k,j,i-1),2) &
759	.OR. &
760	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
761	.AND. i == nxr-1 ) .OR. &
762	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
763	.AND. i == nxlu ) ) &
764	THEN
765	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 19 )
766	!
767	!-- Clear flag for WS1
768	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 18 )
769	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),2) .AND. &
770	BTEST(wall_flags_0(k,j,i+2),2) .AND. &
771	BTEST(wall_flags_0(k,j,i-1),2) ) &
772	THEN
773	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 20 )
774	!
775	!-- Clear flag for WS1
776	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 18 )
777	ENDIF
778	!
779	!-- v component - y-direction
780	!-- WS1 (21), WS3 (22), WS5 (23)
781	IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),2) .OR. &
782	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
783	.AND. j <= nysv ) .OR. &
784	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
785	.AND. j == nyn ) ) &
786	THEN
787	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 21 )
788	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),2) .AND. &
789	BTEST(wall_flags_0(k,j+1,i),2) .OR. &
790	.NOT. BTEST(wall_flags_0(k,j-1,i),2) ) &
791	.OR. &
792	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
793	.AND. j == nysv+1) .OR. &
794	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
795	.AND. j == nyn-1 ) ) &
796	THEN
797	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 22 )
798	!
799	!-- Clear flag for WS1
800	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 21 )
801	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),2) .AND. &
802	BTEST(wall_flags_0(k,j+2,i),2) .AND. &
803	BTEST(wall_flags_0(k,j-1,i),2) ) &
804	THEN
805	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 23 )
806	!
807	!-- Clear flag for WS1
808	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 21 )
809	ENDIF
810	!
811	!-- v component - z-direction
812	!-- WS1 (24), WS3 (25), WS5 (26)
813	IF ( k == nzb+1 ) THEN
814	k_mm = nzb
815	ELSE
816	k_mm = k - 2
817	ENDIF
818	IF ( k > nzt-1 ) THEN
819	k_pp = nzt+1
820	ELSE
821	k_pp = k + 2
822	ENDIF
823	IF ( k > nzt-2 ) THEN
824	k_ppp = nzt+1
825	ELSE
826	k_ppp = k + 3
827	ENDIF
828
829	flag_set = .FALSE.
830	IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),2) .AND. &
831	BTEST(wall_flags_0(k,j,i),2) .OR. &
832	.NOT. BTEST(wall_flags_0(k_pp,j,i),2) .AND. &
833	BTEST(wall_flags_0(k,j,i),2) .OR. &
834	k == nzt ) &
835	THEN
836	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 24 )
837	flag_set = .TRUE.
838	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),2) .OR. &
839	.NOT. BTEST(wall_flags_0(k_ppp,j,i),2) ) .AND. &
840	BTEST(wall_flags_0(k-1,j,i),2) .AND. &
841	BTEST(wall_flags_0(k,j,i),2) .AND. &
842	BTEST(wall_flags_0(k+1,j,i),2) .AND. &
843	BTEST(wall_flags_0(k_pp,j,i),2) .OR. &
844	k == nzt - 1 ) &
845	THEN
846	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 25 )
847	flag_set = .TRUE.
848	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),2) .AND. &
849	BTEST(wall_flags_0(k-1,j,i),2) .AND. &
850	BTEST(wall_flags_0(k,j,i),2) .AND. &
851	BTEST(wall_flags_0(k+1,j,i),2) .AND. &
852	BTEST(wall_flags_0(k_pp,j,i),2) .AND. &
853	BTEST(wall_flags_0(k_ppp,j,i),2) .AND. &
854	.NOT. flag_set ) &
855	THEN
856	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 26 )
857	ENDIF
858
859	ENDDO
860	ENDDO
861	ENDDO
862	DO i = nxl, nxr
863	DO j = nys, nyn
864	DO k = nzb+1, nzt
865	!
866	!-- At first, set flags to WS1.
867	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
868	!-- in order to handle the left/south flux.
869	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 27 )
870	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 30 )
871	!
872	!-- w component - x-direction
873	!-- WS1 (27), WS3 (28), WS5 (29)
874	IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),3) .OR. &
875	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
876	.AND. i == nxl ) .OR. &
877	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
878	.AND. i == nxr ) ) &
879	THEN
880	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 27 )
881	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),3) .AND. &
882	BTEST(wall_flags_0(k,j,i+1),3) .OR. &
883	.NOT. BTEST(wall_flags_0(k,j,i-1),3) ) &
884	.OR. &
885	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
886	.AND. i == nxr-1 ) .OR. &
887	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
888	.AND. i == nxlu ) ) &
889	THEN
890	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 28 )
891	!
892	!-- Clear flag for WS1
893	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 27 )
894	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),3) .AND. &
895	BTEST(wall_flags_0(k,j,i+2),3) .AND. &
896	BTEST(wall_flags_0(k,j,i-1),3) ) &
897	THEN
898	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i),29 )
899	!
900	!-- Clear flag for WS1
901	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 27 )
902	ENDIF
903	!
904	!-- w component - y-direction
905	!-- WS1 (30), WS3 (31), WS5 (32)
906	IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),3) .OR. &
907	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
908	.AND. j == nys ) .OR. &
909	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
910	.AND. j == nyn ) ) &
911	THEN
912	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 30 )
913	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),3) .AND. &
914	BTEST(wall_flags_0(k,j+1,i),3) .OR. &
915	.NOT. BTEST(wall_flags_0(k,j-1,i),3) ) &
916	.OR. &
917	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
918	.AND. j == nysv ) .OR. &
919	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
920	.AND. j == nyn-1 ) ) &
921	THEN
922	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 31 )
923	!
924	!-- Clear flag for WS1
925	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 30 )
926	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),3) .AND. &
927	BTEST(wall_flags_0(k,j+2,i),3) .AND. &
928	BTEST(wall_flags_0(k,j-1,i),3) ) &
929	THEN
930	advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 0 )
931	!
932	!-- Clear flag for WS1
933	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 30 )
934	ENDIF
935	!
936	!-- w component - z-direction
937	!-- WS1 (33), WS3 (34), WS5 (35)
938	flag_set = .FALSE.
939	IF ( k == nzb+1 ) THEN
940	k_mm = nzb
941	ELSE
942	k_mm = k - 2
943	ENDIF
944	IF ( k > nzt-1 ) THEN
945	k_pp = nzt+1
946	ELSE
947	k_pp = k + 2
948	ENDIF
949	IF ( k > nzt-2 ) THEN
950	k_ppp = nzt+1
951	ELSE
952	k_ppp = k + 3
953	ENDIF
954
955	IF ( ( .NOT. BTEST(wall_flags_0(k-1,j,i),3) .AND. &
956	.NOT. BTEST(wall_flags_0(k,j,i),3) .AND. &
957	BTEST(wall_flags_0(k+1,j,i),3) ) .OR. &
958	( .NOT. BTEST(wall_flags_0(k-1,j,i),3) .AND. &
959	BTEST(wall_flags_0(k,j,i),3) ) .OR. &
960	( .NOT. BTEST(wall_flags_0(k+1,j,i),3) .AND. &
961	BTEST(wall_flags_0(k,j,i),3) ) .OR. &
962	k == nzt ) &
963	THEN
964	!
965	!-- Please note, at k == nzb_w_inner(j,i) a flag is explictely
966	!-- set, although this is not a prognostic level. However,
967	!-- contrary to the advection of u,v and s this is necessary
968	!-- because flux_t(nzb_w_inner(j,i)) is used for the tendency
969	!-- at k == nzb_w_inner(j,i)+1.
970	advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 1 )
971	flag_set = .TRUE.
972	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),3) .OR. &
973	.NOT. BTEST(wall_flags_0(k_ppp,j,i),3) ) .AND. &
974	BTEST(wall_flags_0(k-1,j,i),3) .AND. &
975	BTEST(wall_flags_0(k,j,i),3) .AND. &
976	BTEST(wall_flags_0(k+1,j,i),3) .OR. &
977	k == nzt - 1 ) &
978	THEN
979	advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 2 )
980	flag_set = .TRUE.
981	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),3) .AND. &
982	BTEST(wall_flags_0(k-1,j,i),3) .AND. &
983	BTEST(wall_flags_0(k,j,i),3) .AND. &
984	BTEST(wall_flags_0(k+1,j,i),3) .AND. &
985	BTEST(wall_flags_0(k_pp,j,i),3) .AND. &
986	BTEST(wall_flags_0(k_ppp,j,i),3) .AND. &
987	.NOT. flag_set ) &
988	THEN
989	advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 3 )
990	ENDIF
991
992	ENDDO
993	ENDDO
994	ENDDO
995
996	ENDIF
997
998
999	!
1000	!-- Exchange 3D integer wall_flags.
1001	IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' &
1002	) THEN
1003	!
1004	!-- Exchange ghost points for advection flags
1005	CALL exchange_horiz_int( advc_flags_1, nbgp )
1006	CALL exchange_horiz_int( advc_flags_2, nbgp )
1007	!
1008	!-- Set boundary flags at inflow and outflow boundary in case of
1009	!-- non-cyclic boundary conditions.
1010	IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) THEN
1011	advc_flags_1(:,:,nxl-1) = advc_flags_1(:,:,nxl)
1012	advc_flags_2(:,:,nxl-1) = advc_flags_2(:,:,nxl)
1013	ENDIF
1014
1015	IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) THEN
1016	advc_flags_1(:,:,nxr+1) = advc_flags_1(:,:,nxr)
1017	advc_flags_2(:,:,nxr+1) = advc_flags_2(:,:,nxr)
1018	ENDIF
1019
1020	IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) THEN
1021	advc_flags_1(:,nyn+1,:) = advc_flags_1(:,nyn,:)
1022	advc_flags_2(:,nyn+1,:) = advc_flags_2(:,nyn,:)
1023	ENDIF
1024
1025	IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) THEN
1026	advc_flags_1(:,nys-1,:) = advc_flags_1(:,nys,:)
1027	advc_flags_2(:,nys-1,:) = advc_flags_2(:,nys,:)
1028	ENDIF
1029
1030	ENDIF
1031
1032
1033	END SUBROUTINE ws_init_flags
1034
1035
1036	!------------------------------------------------------------------------------!
1037	! Description:
1038	! ------------
1039	!> Initialize variables used for storing statistic quantities (fluxes, variances)
1040	!------------------------------------------------------------------------------!
1041	SUBROUTINE ws_statistics
1042
1043	USE control_parameters, &
1044	ONLY: cloud_physics, humidity, passive_scalar, ocean, &
1045	microphysics_morrison, microphysics_seifert, ws_scheme_mom, &
1046	ws_scheme_sca
1047
1048	USE kinds
1049
1050	USE statistics, &
1051	ONLY: sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsncs_ws_l,&
1052	sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqcs_ws_l, &
1053	sums_wsqrs_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l, &
1054	sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l
1055
1056
1057	IMPLICIT NONE
1058
1059	!
1060	!-- The arrays needed for statistical evaluation are set to to 0 at the
1061	!-- beginning of prognostic_equations.
1062	IF ( ws_scheme_mom ) THEN
1063	sums_wsus_ws_l = 0.0_wp
1064	sums_wsvs_ws_l = 0.0_wp
1065	sums_us2_ws_l = 0.0_wp
1066	sums_vs2_ws_l = 0.0_wp
1067	sums_ws2_ws_l = 0.0_wp
1068	ENDIF
1069
1070	IF ( ws_scheme_sca ) THEN
1071	sums_wspts_ws_l = 0.0_wp
1072	IF ( humidity ) sums_wsqs_ws_l = 0.0_wp
1073	IF ( passive_scalar ) sums_wsss_ws_l = 0.0_wp
1074	IF ( cloud_physics .AND. microphysics_morrison ) THEN
1075	sums_wsqcs_ws_l = 0.0_wp
1076	sums_wsncs_ws_l = 0.0_wp
1077	ENDIF
1078	IF ( cloud_physics .AND. microphysics_seifert ) THEN
1079	sums_wsqrs_ws_l = 0.0_wp
1080	sums_wsnrs_ws_l = 0.0_wp
1081	ENDIF
1082	IF ( ocean ) sums_wssas_ws_l = 0.0_wp
1083
1084	ENDIF
1085
1086	END SUBROUTINE ws_statistics
1087
1088
1089	!------------------------------------------------------------------------------!
1090	! Description:
1091	! ------------
1092	!> Scalar advection - Call for grid point i,j
1093	!------------------------------------------------------------------------------!
1094	SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local, &
1095	swap_diss_y_local, swap_flux_x_local, &
1096	swap_diss_x_local, i_omp, tn )
1097
1098	USE arrays_3d, &
1099	ONLY: ddzw, drho_air, tend, u, v, w, rho_air_zw
1100
1101	USE constants, &
1102	ONLY: adv_sca_1, adv_sca_3, adv_sca_5
1103
1104	USE control_parameters, &
1105	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
1106
1107	USE grid_variables, &
1108	ONLY: ddx, ddy
1109
1110	USE indices, &
1111	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, &
1112	nzt, advc_flags_1
1113
1114	USE kinds
1115
1116	USE pegrid
1117
1118	USE statistics, &
1119	ONLY: hom, sums_wsncs_ws_l, sums_wsnrs_ws_l, sums_wspts_ws_l, &
1120	sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l, &
1121	sums_wssas_ws_l, sums_wsss_ws_l, weight_substep
1122
1123
1124	IMPLICIT NONE
1125
1126	CHARACTER (LEN = *), INTENT(IN) :: sk_char !<
1127
1128	INTEGER(iwp) :: i !<
1129	INTEGER(iwp) :: ibit0 !<
1130	INTEGER(iwp) :: ibit1 !<
1131	INTEGER(iwp) :: ibit2 !<
1132	INTEGER(iwp) :: ibit3 !<
1133	INTEGER(iwp) :: ibit4 !<
1134	INTEGER(iwp) :: ibit5 !<
1135	INTEGER(iwp) :: ibit6 !<
1136	INTEGER(iwp) :: ibit7 !<
1137	INTEGER(iwp) :: ibit8 !<
1138	INTEGER(iwp) :: i_omp !<
1139	INTEGER(iwp) :: j !<
1140	INTEGER(iwp) :: k !<
1141	INTEGER(iwp) :: k_mm !<
1142	INTEGER(iwp) :: k_mmm !<
1143	INTEGER(iwp) :: k_pp !<
1144	INTEGER(iwp) :: k_ppp !<
1145	INTEGER(iwp) :: tn !<
1146
1147	REAL(wp) :: diss_d !<
1148	REAL(wp) :: div !<
1149	REAL(wp) :: flux_d !<
1150	REAL(wp) :: u_comp !<
1151	REAL(wp) :: v_comp !<
1152
1153	#if defined( __nopointer )
1154	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !<
1155	#else
1156	REAL(wp), DIMENSION(:,:,:), POINTER :: sk !<
1157	#endif
1158	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !<
1159	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !<
1160	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !<
1161	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !<
1162	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !<
1163	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !<
1164
1165	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local !<
1166	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local !<
1167
1168	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_diss_x_local !<
1169	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_flux_x_local !<
1170
1171
1172	!
1173	!-- Compute southside fluxes of the respective PE bounds.
1174	IF ( j == nys ) THEN
1175	!
1176	!-- Up to the top of the highest topography.
1177	DO k = nzb+1, nzb_max
1178
1179	ibit5 = IBITS(advc_flags_1(k,j-1,i),5,1)
1180	ibit4 = IBITS(advc_flags_1(k,j-1,i),4,1)
1181	ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
1182
1183	v_comp = v(k,j,i) - v_gtrans
1184	swap_flux_y_local(k,tn) = v_comp * ( &
1185	( 37.0_wp * ibit5 * adv_sca_5 &
1186	+ 7.0_wp * ibit4 * adv_sca_3 &
1187	+ ibit3 * adv_sca_1 &
1188	) * &
1189	( sk(k,j,i) + sk(k,j-1,i) ) &
1190	- ( 8.0_wp * ibit5 * adv_sca_5 &
1191	+ ibit4 * adv_sca_3 &
1192	) * &
1193	( sk(k,j+1,i) + sk(k,j-2,i) ) &
1194	+ ( ibit5 * adv_sca_5 &
1195	) * &
1196	( sk(k,j+2,i) + sk(k,j-3,i) ) &
1197	)
1198
1199	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
1200	( 10.0_wp * ibit5 * adv_sca_5 &
1201	+ 3.0_wp * ibit4 * adv_sca_3 &
1202	+ ibit3 * adv_sca_1 &
1203	) * &
1204	( sk(k,j,i) - sk(k,j-1,i) ) &
1205	- ( 5.0_wp * ibit5 * adv_sca_5 &
1206	+ ibit4 * adv_sca_3 &
1207	) * &
1208	( sk(k,j+1,i) - sk(k,j-2,i) ) &
1209	+ ( ibit5 * adv_sca_5 &
1210	) * &
1211	( sk(k,j+2,i) - sk(k,j-3,i) ) &
1212	)
1213
1214	ENDDO
1215	!
1216	!-- Above to the top of the highest topography. No degradation necessary.
1217	DO k = nzb_max+1, nzt
1218
1219	v_comp = v(k,j,i) - v_gtrans
1220	swap_flux_y_local(k,tn) = v_comp * ( &
1221	37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) &
1222	- 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
1223	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) &
1224	) * adv_sca_5
1225	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
1226	10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) &
1227	- 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
1228	+ sk(k,j+2,i) - sk(k,j-3,i) &
1229	) * adv_sca_5
1230
1231	ENDDO
1232
1233	ENDIF
1234	!
1235	!-- Compute leftside fluxes of the respective PE bounds.
1236	IF ( i == i_omp ) THEN
1237
1238	DO k = nzb+1, nzb_max
1239
1240	ibit2 = IBITS(advc_flags_1(k,j,i-1),2,1)
1241	ibit1 = IBITS(advc_flags_1(k,j,i-1),1,1)
1242	ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
1243
1244	u_comp = u(k,j,i) - u_gtrans
1245	swap_flux_x_local(k,j,tn) = u_comp * ( &
1246	( 37.0_wp * ibit2 * adv_sca_5 &
1247	+ 7.0_wp * ibit1 * adv_sca_3 &
1248	+ ibit0 * adv_sca_1 &
1249	) * &
1250	( sk(k,j,i) + sk(k,j,i-1) ) &
1251	- ( 8.0_wp * ibit2 * adv_sca_5 &
1252	+ ibit1 * adv_sca_3 &
1253	) * &
1254	( sk(k,j,i+1) + sk(k,j,i-2) ) &
1255	+ ( ibit2 * adv_sca_5 &
1256	) * &
1257	( sk(k,j,i+2) + sk(k,j,i-3) ) &
1258	)
1259
1260	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
1261	( 10.0_wp * ibit2 * adv_sca_5 &
1262	+ 3.0_wp * ibit1 * adv_sca_3 &
1263	+ ibit0 * adv_sca_1 &
1264	) * &
1265	( sk(k,j,i) - sk(k,j,i-1) ) &
1266	- ( 5.0_wp * ibit2 * adv_sca_5 &
1267	+ ibit1 * adv_sca_3 &
1268	) * &
1269	( sk(k,j,i+1) - sk(k,j,i-2) ) &
1270	+ ( ibit2 * adv_sca_5 &
1271	) * &
1272	( sk(k,j,i+2) - sk(k,j,i-3) ) &
1273	)
1274
1275	ENDDO
1276
1277	DO k = nzb_max+1, nzt
1278
1279	u_comp = u(k,j,i) - u_gtrans
1280	swap_flux_x_local(k,j,tn) = u_comp * ( &
1281	37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) &
1282	- 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
1283	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) &
1284	) * adv_sca_5
1285
1286	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
1287	10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) &
1288	- 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
1289	+ ( sk(k,j,i+2) - sk(k,j,i-3) ) &
1290	) * adv_sca_5
1291
1292	ENDDO
1293
1294	ENDIF
1295
1296	flux_t(0) = 0.0_wp
1297	diss_t(0) = 0.0_wp
1298	flux_d = 0.0_wp
1299	diss_d = 0.0_wp
1300	!
1301	!-- Now compute the fluxes and tendency terms for the horizontal and
1302	!-- vertical parts up to the top of the highest topography.
1303	DO k = nzb+1, nzb_max
1304	!
1305	!-- Note: It is faster to conduct all multiplications explicitly, e.g.
1306	!-- * adv_sca_5 ... than to determine a factor and multiplicate the
1307	!-- flux at the end.
1308
1309	ibit2 = IBITS(advc_flags_1(k,j,i),2,1)
1310	ibit1 = IBITS(advc_flags_1(k,j,i),1,1)
1311	ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
1312
1313	u_comp = u(k,j,i+1) - u_gtrans
1314	flux_r(k) = u_comp * ( &
1315	( 37.0_wp * ibit2 * adv_sca_5 &
1316	+ 7.0_wp * ibit1 * adv_sca_3 &
1317	+ ibit0 * adv_sca_1 &
1318	) * &
1319	( sk(k,j,i+1) + sk(k,j,i) ) &
1320	- ( 8.0_wp * ibit2 * adv_sca_5 &
1321	+ ibit1 * adv_sca_3 &
1322	) * &
1323	( sk(k,j,i+2) + sk(k,j,i-1) ) &
1324	+ ( ibit2 * adv_sca_5 &
1325	) * &
1326	( sk(k,j,i+3) + sk(k,j,i-2) ) &
1327	)
1328
1329	diss_r(k) = -ABS( u_comp ) * ( &
1330	( 10.0_wp * ibit2 * adv_sca_5 &
1331	+ 3.0_wp * ibit1 * adv_sca_3 &
1332	+ ibit0 * adv_sca_1 &
1333	) * &
1334	( sk(k,j,i+1) - sk(k,j,i) ) &
1335	- ( 5.0_wp * ibit2 * adv_sca_5 &
1336	+ ibit1 * adv_sca_3 &
1337	) * &
1338	( sk(k,j,i+2) - sk(k,j,i-1) ) &
1339	+ ( ibit2 * adv_sca_5 &
1340	) * &
1341	( sk(k,j,i+3) - sk(k,j,i-2) ) &
1342	)
1343
1344	ibit5 = IBITS(advc_flags_1(k,j,i),5,1)
1345	ibit4 = IBITS(advc_flags_1(k,j,i),4,1)
1346	ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
1347
1348	v_comp = v(k,j+1,i) - v_gtrans
1349	flux_n(k) = v_comp * ( &
1350	( 37.0_wp * ibit5 * adv_sca_5 &
1351	+ 7.0_wp * ibit4 * adv_sca_3 &
1352	+ ibit3 * adv_sca_1 &
1353	) * &
1354	( sk(k,j+1,i) + sk(k,j,i) ) &
1355	- ( 8.0_wp * ibit5 * adv_sca_5 &
1356	+ ibit4 * adv_sca_3 &
1357	) * &
1358	( sk(k,j+2,i) + sk(k,j-1,i) ) &
1359	+ ( ibit5 * adv_sca_5 &
1360	) * &
1361	( sk(k,j+3,i) + sk(k,j-2,i) ) &
1362	)
1363
1364	diss_n(k) = -ABS( v_comp ) * ( &
1365	( 10.0_wp * ibit5 * adv_sca_5 &
1366	+ 3.0_wp * ibit4 * adv_sca_3 &
1367	+ ibit3 * adv_sca_1 &
1368	) * &
1369	( sk(k,j+1,i) - sk(k,j,i) ) &
1370	- ( 5.0_wp * ibit5 * adv_sca_5 &
1371	+ ibit4 * adv_sca_3 &
1372	) * &
1373	( sk(k,j+2,i) - sk(k,j-1,i) ) &
1374	+ ( ibit5 * adv_sca_5 &
1375	) * &
1376	( sk(k,j+3,i) - sk(k,j-2,i) ) &
1377	)
1378	!
1379	!-- k index has to be modified near bottom and top, else array
1380	!-- subscripts will be exceeded.
1381	ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
1382	ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
1383	ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
1384
1385	k_ppp = k + 3 * ibit8
1386	k_pp = k + 2 * ( 1 - ibit6 )
1387	k_mm = k - 2 * ibit8
1388
1389
1390	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
1391	( 37.0_wp * ibit8 * adv_sca_5 &
1392	+ 7.0_wp * ibit7 * adv_sca_3 &
1393	+ ibit6 * adv_sca_1 &
1394	) * &
1395	( sk(k+1,j,i) + sk(k,j,i) ) &
1396	- ( 8.0_wp * ibit8 * adv_sca_5 &
1397	+ ibit7 * adv_sca_3 &
1398	) * &
1399	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
1400	+ ( ibit8 * adv_sca_5 &
1401	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
1402	)
1403
1404	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
1405	( 10.0_wp * ibit8 * adv_sca_5 &
1406	+ 3.0_wp * ibit7 * adv_sca_3 &
1407	+ ibit6 * adv_sca_1 &
1408	) * &
1409	( sk(k+1,j,i) - sk(k,j,i) ) &
1410	- ( 5.0_wp * ibit8 * adv_sca_5 &
1411	+ ibit7 * adv_sca_3 &
1412	) * &
1413	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
1414	+ ( ibit8 * adv_sca_5 &
1415	) * &
1416	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
1417	)
1418	!
1419	!-- Calculate the divergence of the velocity field. A respective
1420	!-- correction is needed to overcome numerical instabilities caused
1421	!-- by a not sufficient reduction of divergences near topography.
1422	div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) &
1423	- u(k,j,i) * ( IBITS(advc_flags_1(k,j,i-1),0,1) &
1424	+ IBITS(advc_flags_1(k,j,i-1),1,1) &
1425	+ IBITS(advc_flags_1(k,j,i-1),2,1) &
1426	) &
1427	) * ddx &
1428	+ ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) &
1429	- v(k,j,i) * ( IBITS(advc_flags_1(k,j-1,i),3,1) &
1430	+ IBITS(advc_flags_1(k,j-1,i),4,1) &
1431	+ IBITS(advc_flags_1(k,j-1,i),5,1) &
1432	) &
1433	) * ddy &
1434	+ ( w(k,j,i) * rho_air_zw(k) * &
1435	( ibit6 + ibit7 + ibit8 ) &
1436	- w(k-1,j,i) * rho_air_zw(k-1) * &
1437	( IBITS(advc_flags_1(k-1,j,i),6,1) &
1438	+ IBITS(advc_flags_1(k-1,j,i),7,1) &
1439	+ IBITS(advc_flags_1(k-1,j,i),8,1) &
1440	) &
1441	) * drho_air(k) * ddzw(k)
1442
1443
1444	tend(k,j,i) = tend(k,j,i) - ( &
1445	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
1446	swap_diss_x_local(k,j,tn) ) * ddx &
1447	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
1448	swap_diss_y_local(k,tn) ) * ddy &
1449	+ ( ( flux_t(k) + diss_t(k) ) - &
1450	( flux_d + diss_d ) &
1451	) * drho_air(k) * ddzw(k) &
1452	) + sk(k,j,i) * div
1453
1454	swap_flux_y_local(k,tn) = flux_n(k)
1455	swap_diss_y_local(k,tn) = diss_n(k)
1456	swap_flux_x_local(k,j,tn) = flux_r(k)
1457	swap_diss_x_local(k,j,tn) = diss_r(k)
1458	flux_d = flux_t(k)
1459	diss_d = diss_t(k)
1460
1461	ENDDO
1462	!
1463	!-- Now compute the fluxes and tendency terms for the horizontal and
1464	!-- vertical parts above the top of the highest topography. No degradation
1465	!-- for the horizontal parts, but for the vertical it is stell needed.
1466	DO k = nzb_max+1, nzt
1467
1468	u_comp = u(k,j,i+1) - u_gtrans
1469	flux_r(k) = u_comp * ( &
1470	37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) &
1471	- 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
1472	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
1473	diss_r(k) = -ABS( u_comp ) * ( &
1474	10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) &
1475	- 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
1476	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
1477
1478	v_comp = v(k,j+1,i) - v_gtrans
1479	flux_n(k) = v_comp * ( &
1480	37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) &
1481	- 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
1482	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
1483	diss_n(k) = -ABS( v_comp ) * ( &
1484	10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) &
1485	- 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
1486	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
1487	!
1488	!-- k index has to be modified near bottom and top, else array
1489	!-- subscripts will be exceeded.
1490	ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
1491	ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
1492	ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
1493
1494	k_ppp = k + 3 * ibit8
1495	k_pp = k + 2 * ( 1 - ibit6 )
1496	k_mm = k - 2 * ibit8
1497
1498	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
1499	( 37.0_wp * ibit8 * adv_sca_5 &
1500	+ 7.0_wp * ibit7 * adv_sca_3 &
1501	+ ibit6 * adv_sca_1 &
1502	) * &
1503	( sk(k+1,j,i) + sk(k,j,i) ) &
1504	- ( 8.0_wp * ibit8 * adv_sca_5 &
1505	+ ibit7 * adv_sca_3 &
1506	) * &
1507	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
1508	+ ( ibit8 * adv_sca_5 &
1509	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
1510	)
1511
1512	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
1513	( 10.0_wp * ibit8 * adv_sca_5 &
1514	+ 3.0_wp * ibit7 * adv_sca_3 &
1515	+ ibit6 * adv_sca_1 &
1516	) * &
1517	( sk(k+1,j,i) - sk(k,j,i) ) &
1518	- ( 5.0_wp * ibit8 * adv_sca_5 &
1519	+ ibit7 * adv_sca_3 &
1520	) * &
1521	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
1522	+ ( ibit8 * adv_sca_5 &
1523	) * &
1524	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
1525	)
1526	!
1527	!-- Calculate the divergence of the velocity field. A respective
1528	!-- correction is needed to overcome numerical instabilities introduced
1529	!-- by a not sufficient reduction of divergences near topography.
1530	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
1531	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
1532	+ ( w(k,j,i) * rho_air_zw(k) - &
1533	w(k-1,j,i) * rho_air_zw(k-1) &
1534	) * drho_air(k) * ddzw(k)
1535
1536	tend(k,j,i) = tend(k,j,i) - ( &
1537	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
1538	swap_diss_x_local(k,j,tn) ) * ddx &
1539	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
1540	swap_diss_y_local(k,tn) ) * ddy &
1541	+ ( ( flux_t(k) + diss_t(k) ) - &
1542	( flux_d + diss_d ) &
1543	) * drho_air(k) * ddzw(k) &
1544	) + sk(k,j,i) * div
1545
1546
1547	swap_flux_y_local(k,tn) = flux_n(k)
1548	swap_diss_y_local(k,tn) = diss_n(k)
1549	swap_flux_x_local(k,j,tn) = flux_r(k)
1550	swap_diss_x_local(k,j,tn) = diss_r(k)
1551	flux_d = flux_t(k)
1552	diss_d = diss_t(k)
1553
1554	ENDDO
1555
1556	!
1557	!-- Evaluation of statistics.
1558	SELECT CASE ( sk_char )
1559
1560	CASE ( 'pt' )
1561
1562	DO k = nzb, nzt
1563	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + &
1564	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1565	* ( w(k,j,i) - hom(k,1,3,0) ) &
1566	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1567	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1568	) * weight_substep(intermediate_timestep_count)
1569	ENDDO
1570
1571	CASE ( 'sa' )
1572
1573	DO k = nzb, nzt
1574	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) + &
1575	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1576	* ( w(k,j,i) - hom(k,1,3,0) ) &
1577	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1578	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1579	) * weight_substep(intermediate_timestep_count)
1580	ENDDO
1581
1582	CASE ( 'q' )
1583
1584	DO k = nzb, nzt
1585	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) + &
1586	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1587	* ( w(k,j,i) - hom(k,1,3,0) ) &
1588	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1589	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1590	) * weight_substep(intermediate_timestep_count)
1591	ENDDO
1592
1593	CASE ( 'qc' )
1594
1595	DO k = nzb, nzt
1596	sums_wsqcs_ws_l(k,tn) = sums_wsqcs_ws_l(k,tn) + &
1597	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1598	* ( w(k,j,i) - hom(k,1,3,0) ) &
1599	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1600	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1601	) * weight_substep(intermediate_timestep_count)
1602	ENDDO
1603
1604
1605	CASE ( 'qr' )
1606
1607	DO k = nzb, nzt
1608	sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) + &
1609	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1610	* ( w(k,j,i) - hom(k,1,3,0) ) &
1611	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1612	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1613	) * weight_substep(intermediate_timestep_count)
1614	ENDDO
1615
1616	CASE ( 'nc' )
1617
1618	DO k = nzb, nzt
1619	sums_wsncs_ws_l(k,tn) = sums_wsncs_ws_l(k,tn) + &
1620	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1621	* ( w(k,j,i) - hom(k,1,3,0) ) &
1622	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1623	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1624	) * weight_substep(intermediate_timestep_count)
1625	ENDDO
1626
1627	CASE ( 'nr' )
1628
1629	DO k = nzb, nzt
1630	sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) + &
1631	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1632	* ( w(k,j,i) - hom(k,1,3,0) ) &
1633	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1634	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1635	) * weight_substep(intermediate_timestep_count)
1636	ENDDO
1637
1638	CASE ( 's' )
1639
1640	DO k = nzb, nzt
1641	sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn) + &
1642	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1643	* ( w(k,j,i) - hom(k,1,3,0) ) &
1644	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1645	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1646	) * weight_substep(intermediate_timestep_count)
1647	ENDDO
1648
1649	END SELECT
1650
1651	END SUBROUTINE advec_s_ws_ij
1652
1653
1654
1655
1656	!------------------------------------------------------------------------------!
1657	! Description:
1658	! ------------
1659	!> Advection of u-component - Call for grid point i,j
1660	!------------------------------------------------------------------------------!
1661	SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
1662
1663	USE arrays_3d, &
1664	ONLY: ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w,&
1665	drho_air, rho_air_zw
1666
1667	USE constants, &
1668	ONLY: adv_mom_1, adv_mom_3, adv_mom_5
1669
1670	USE control_parameters, &
1671	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
1672
1673	USE grid_variables, &
1674	ONLY: ddx, ddy
1675
1676	USE indices, &
1677	ONLY: nxl, nxr, nyn, nys, nzb, nzb_max, nzt, advc_flags_1
1678
1679	USE kinds
1680
1681	USE statistics, &
1682	ONLY: hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep
1683
1684	IMPLICIT NONE
1685
1686	INTEGER(iwp) :: i !<
1687	INTEGER(iwp) :: ibit9 !<
1688	INTEGER(iwp) :: ibit10 !<
1689	INTEGER(iwp) :: ibit11 !<
1690	INTEGER(iwp) :: ibit12 !<
1691	INTEGER(iwp) :: ibit13 !<
1692	INTEGER(iwp) :: ibit14 !<
1693	INTEGER(iwp) :: ibit15 !<
1694	INTEGER(iwp) :: ibit16 !<
1695	INTEGER(iwp) :: ibit17 !<
1696	INTEGER(iwp) :: i_omp !<
1697	INTEGER(iwp) :: j !<
1698	INTEGER(iwp) :: k !<
1699	INTEGER(iwp) :: k_mm !<
1700	INTEGER(iwp) :: k_pp !<
1701	INTEGER(iwp) :: k_ppp !<
1702	INTEGER(iwp) :: tn !<
1703
1704	REAL(wp) :: diss_d !<
1705	REAL(wp) :: div !<
1706	REAL(wp) :: flux_d !<
1707	REAL(wp) :: gu !<
1708	REAL(wp) :: gv !<
1709	REAL(wp) :: u_comp_l !<
1710	REAL(wp) :: v_comp !<
1711	REAL(wp) :: w_comp !<
1712
1713	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !<
1714	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !<
1715	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !<
1716	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !<
1717	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !<
1718	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !<
1719	REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !<
1720
1721	gu = 2.0_wp * u_gtrans
1722	gv = 2.0_wp * v_gtrans
1723	!
1724	!-- Compute southside fluxes for the respective boundary of PE
1725	IF ( j == nys ) THEN
1726
1727	DO k = nzb+1, nzb_max
1728
1729	ibit14 = IBITS(advc_flags_1(k,j-1,i),14,1)
1730	ibit13 = IBITS(advc_flags_1(k,j-1,i),13,1)
1731	ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)
1732
1733	v_comp = v(k,j,i) + v(k,j,i-1) - gv
1734	flux_s_u(k,tn) = v_comp * ( &
1735	( 37.0_wp * ibit14 * adv_mom_5 &
1736	+ 7.0_wp * ibit13 * adv_mom_3 &
1737	+ ibit12 * adv_mom_1 &
1738	) * &
1739	( u(k,j,i) + u(k,j-1,i) ) &
1740	- ( 8.0_wp * ibit14 * adv_mom_5 &
1741	+ ibit13 * adv_mom_3 &
1742	) * &
1743	( u(k,j+1,i) + u(k,j-2,i) ) &
1744	+ ( ibit14 * adv_mom_5 &
1745	) * &
1746	( u(k,j+2,i) + u(k,j-3,i) ) &
1747	)
1748
1749	diss_s_u(k,tn) = - ABS ( v_comp ) * ( &
1750	( 10.0_wp * ibit14 * adv_mom_5 &
1751	+ 3.0_wp * ibit13 * adv_mom_3 &
1752	+ ibit12 * adv_mom_1 &
1753	) * &
1754	( u(k,j,i) - u(k,j-1,i) ) &
1755	- ( 5.0_wp * ibit14 * adv_mom_5 &
1756	+ ibit13 * adv_mom_3 &
1757	) * &
1758	( u(k,j+1,i) - u(k,j-2,i) ) &
1759	+ ( ibit14 * adv_mom_5 &
1760	) * &
1761	( u(k,j+2,i) - u(k,j-3,i) ) &
1762	)
1763
1764	ENDDO
1765
1766	DO k = nzb_max+1, nzt
1767
1768	v_comp = v(k,j,i) + v(k,j,i-1) - gv
1769	flux_s_u(k,tn) = v_comp * ( &
1770	37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) &
1771	- 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) &
1772	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
1773	diss_s_u(k,tn) = - ABS(v_comp) * ( &
1774	10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) &
1775	- 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) &
1776	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
1777
1778	ENDDO
1779
1780	ENDIF
1781	!
1782	!-- Compute leftside fluxes for the respective boundary of PE
1783	IF ( i == i_omp ) THEN
1784
1785	DO k = nzb+1, nzb_max
1786
1787	ibit11 = IBITS(advc_flags_1(k,j,i-1),11,1)
1788	ibit10 = IBITS(advc_flags_1(k,j,i-1),10,1)
1789	ibit9 = IBITS(advc_flags_1(k,j,i-1),9,1)
1790
1791	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
1792	flux_l_u(k,j,tn) = u_comp_l * ( &
1793	( 37.0_wp * ibit11 * adv_mom_5 &
1794	+ 7.0_wp * ibit10 * adv_mom_3 &
1795	+ ibit9 * adv_mom_1 &
1796	) * &
1797	( u(k,j,i) + u(k,j,i-1) ) &
1798	- ( 8.0_wp * ibit11 * adv_mom_5 &
1799	+ ibit10 * adv_mom_3 &
1800	) * &
1801	( u(k,j,i+1) + u(k,j,i-2) ) &
1802	+ ( ibit11 * adv_mom_5 &
1803	) * &
1804	( u(k,j,i+2) + u(k,j,i-3) ) &
1805	)
1806
1807	diss_l_u(k,j,tn) = - ABS( u_comp_l ) * ( &
1808	( 10.0_wp * ibit11 * adv_mom_5 &
1809	+ 3.0_wp * ibit10 * adv_mom_3 &
1810	+ ibit9 * adv_mom_1 &
1811	) * &
1812	( u(k,j,i) - u(k,j,i-1) ) &
1813	- ( 5.0_wp * ibit11 * adv_mom_5 &
1814	+ ibit10 * adv_mom_3 &
1815	) * &
1816	( u(k,j,i+1) - u(k,j,i-2) ) &
1817	+ ( ibit11 * adv_mom_5 &
1818	) * &
1819	( u(k,j,i+2) - u(k,j,i-3) ) &
1820	)
1821
1822	ENDDO
1823
1824	DO k = nzb_max+1, nzt
1825
1826	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
1827	flux_l_u(k,j,tn) = u_comp_l * ( &
1828	37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) &
1829	- 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) &
1830	+ ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
1831	diss_l_u(k,j,tn) = - ABS(u_comp_l) * ( &
1832	10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) &
1833	- 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) &
1834	+ ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
1835
1836	ENDDO
1837
1838	ENDIF
1839
1840	flux_t(0) = 0.0_wp
1841	diss_t(0) = 0.0_wp
1842	flux_d = 0.0_wp
1843	diss_d = 0.0_wp
1844	!
1845	!-- Now compute the fluxes tendency terms for the horizontal and
1846	!-- vertical parts.
1847	DO k = nzb+1, nzb_max
1848
1849	ibit11 = IBITS(advc_flags_1(k,j,i),11,1)
1850	ibit10 = IBITS(advc_flags_1(k,j,i),10,1)
1851	ibit9 = IBITS(advc_flags_1(k,j,i),9,1)
1852
1853	u_comp(k) = u(k,j,i+1) + u(k,j,i)
1854	flux_r(k) = ( u_comp(k) - gu ) * ( &
1855	( 37.0_wp * ibit11 * adv_mom_5 &
1856	+ 7.0_wp * ibit10 * adv_mom_3 &
1857	+ ibit9 * adv_mom_1 &
1858	) * &
1859	( u(k,j,i+1) + u(k,j,i) ) &
1860	- ( 8.0_wp * ibit11 * adv_mom_5 &
1861	+ ibit10 * adv_mom_3 &
1862	) * &
1863	( u(k,j,i+2) + u(k,j,i-1) ) &
1864	+ ( ibit11 * adv_mom_5 &
1865	) * &
1866	( u(k,j,i+3) + u(k,j,i-2) ) &
1867	)
1868
1869	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
1870	( 10.0_wp * ibit11 * adv_mom_5 &
1871	+ 3.0_wp * ibit10 * adv_mom_3 &
1872	+ ibit9 * adv_mom_1 &
1873	) * &
1874	( u(k,j,i+1) - u(k,j,i) ) &
1875	- ( 5.0_wp * ibit11 * adv_mom_5 &
1876	+ ibit10 * adv_mom_3 &
1877	) * &
1878	( u(k,j,i+2) - u(k,j,i-1) ) &
1879	+ ( ibit11 * adv_mom_5 &
1880	) * &
1881	( u(k,j,i+3) - u(k,j,i-2) ) &
1882	)
1883
1884	ibit14 = IBITS(advc_flags_1(k,j,i),14,1)
1885	ibit13 = IBITS(advc_flags_1(k,j,i),13,1)
1886	ibit12 = IBITS(advc_flags_1(k,j,i),12,1)
1887
1888	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
1889	flux_n(k) = v_comp * ( &
1890	( 37.0_wp * ibit14 * adv_mom_5 &
1891	+ 7.0_wp * ibit13 * adv_mom_3 &
1892	+ ibit12 * adv_mom_1 &
1893	) * &
1894	( u(k,j+1,i) + u(k,j,i) ) &
1895	- ( 8.0_wp * ibit14 * adv_mom_5 &
1896	+ ibit13 * adv_mom_3 &
1897	) * &
1898	( u(k,j+2,i) + u(k,j-1,i) ) &
1899	+ ( ibit14 * adv_mom_5 &
1900	) * &
1901	( u(k,j+3,i) + u(k,j-2,i) ) &
1902	)
1903
1904	diss_n(k) = - ABS ( v_comp ) * ( &
1905	( 10.0_wp * ibit14 * adv_mom_5 &
1906	+ 3.0_wp * ibit13 * adv_mom_3 &
1907	+ ibit12 * adv_mom_1 &
1908	) * &
1909	( u(k,j+1,i) - u(k,j,i) ) &
1910	- ( 5.0_wp * ibit14 * adv_mom_5 &
1911	+ ibit13 * adv_mom_3 &
1912	) * &
1913	( u(k,j+2,i) - u(k,j-1,i) ) &
1914	+ ( ibit14 * adv_mom_5 &
1915	) * &
1916	( u(k,j+3,i) - u(k,j-2,i) ) &
1917	)
1918	!
1919	!-- k index has to be modified near bottom and top, else array
1920	!-- subscripts will be exceeded.
1921	ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
1922	ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
1923	ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
1924
1925	k_ppp = k + 3 * ibit17
1926	k_pp = k + 2 * ( 1 - ibit15 )
1927	k_mm = k - 2 * ibit17
1928
1929	w_comp = w(k,j,i) + w(k,j,i-1)
1930	flux_t(k) = w_comp * rho_air_zw(k) * ( &
1931	( 37.0_wp * ibit17 * adv_mom_5 &
1932	+ 7.0_wp * ibit16 * adv_mom_3 &
1933	+ ibit15 * adv_mom_1 &
1934	) * &
1935	( u(k+1,j,i) + u(k,j,i) ) &
1936	- ( 8.0_wp * ibit17 * adv_mom_5 &
1937	+ ibit16 * adv_mom_3 &
1938	) * &
1939	( u(k_pp,j,i) + u(k-1,j,i) ) &
1940	+ ( ibit17 * adv_mom_5 &
1941	) * &
1942	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
1943	)
1944
1945	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
1946	( 10.0_wp * ibit17 * adv_mom_5 &
1947	+ 3.0_wp * ibit16 * adv_mom_3 &
1948	+ ibit15 * adv_mom_1 &
1949	) * &
1950	( u(k+1,j,i) - u(k,j,i) ) &
1951	- ( 5.0_wp * ibit17 * adv_mom_5 &
1952	+ ibit16 * adv_mom_3 &
1953	) * &
1954	( u(k_pp,j,i) - u(k-1,j,i) ) &
1955	+ ( ibit17 * adv_mom_5 &
1956	) * &
1957	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
1958	)
1959	!
1960	!-- Calculate the divergence of the velocity field. A respective
1961	!-- correction is needed to overcome numerical instabilities introduced
1962	!-- by a not sufficient reduction of divergences near topography.
1963	div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 ) &
1964	- ( u(k,j,i) + u(k,j,i-1) ) &
1965	* ( IBITS(advc_flags_1(k,j,i-1),9,1) &
1966	+ IBITS(advc_flags_1(k,j,i-1),10,1) &
1967	+ IBITS(advc_flags_1(k,j,i-1),11,1) &
1968	) &
1969	) * ddx &
1970	+ ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 ) &
1971	- ( v(k,j,i) + v(k,j,i-1 ) ) &
1972	* ( IBITS(advc_flags_1(k,j-1,i),12,1) &
1973	+ IBITS(advc_flags_1(k,j-1,i),13,1) &
1974	+ IBITS(advc_flags_1(k,j-1,i),14,1) &
1975	) &
1976	) * ddy &
1977	+ ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 ) &
1978	- ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) &
1979	* ( IBITS(advc_flags_1(k-1,j,i),15,1) &
1980	+ IBITS(advc_flags_1(k-1,j,i),16,1) &
1981	+ IBITS(advc_flags_1(k-1,j,i),17,1) &
1982	) &
1983	) * drho_air(k) * ddzw(k) &
1984	) * 0.5_wp
1985
1986
1987	tend(k,j,i) = tend(k,j,i) - ( &
1988	( flux_r(k) + diss_r(k) &
1989	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
1990	+ ( flux_n(k) + diss_n(k) &
1991	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &
1992	+ ( ( flux_t(k) + diss_t(k) ) &
1993	- ( flux_d + diss_d ) &
1994	) * drho_air(k) * ddzw(k) &
1995	) + div * u(k,j,i)
1996
1997	flux_l_u(k,j,tn) = flux_r(k)
1998	diss_l_u(k,j,tn) = diss_r(k)
1999	flux_s_u(k,tn) = flux_n(k)
2000	diss_s_u(k,tn) = diss_n(k)
2001	flux_d = flux_t(k)
2002	diss_d = diss_t(k)
2003	!
2004	!-- Statistical Evaluation of u'u'. The factor has to be applied for
2005	!-- right evaluation when gallilei_trans = .T. .
2006	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
2007	+ ( flux_r(k) &
2008	* ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2009	/ ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &
2010	+ diss_r(k) &
2011	* ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2012	/ ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &
2013	) * weight_substep(intermediate_timestep_count)
2014	!
2015	!-- Statistical Evaluation of w'u'.
2016	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
2017	+ ( flux_t(k) &
2018	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2019	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
2020	+ diss_t(k) &
2021	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2022	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
2023	) * weight_substep(intermediate_timestep_count)
2024	ENDDO
2025
2026	DO k = nzb_max+1, nzt
2027
2028	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2029	flux_r(k) = ( u_comp(k) - gu ) * ( &
2030	37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) &
2031	- 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) &
2032	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
2033	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
2034	10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) &
2035	- 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) &
2036	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
2037
2038	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2039	flux_n(k) = v_comp * ( &
2040	37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) &
2041	- 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) &
2042	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
2043	diss_n(k) = - ABS( v_comp ) * ( &
2044	10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) &
2045	- 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) &
2046	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
2047	!
2048	!-- k index has to be modified near bottom and top, else array
2049	!-- subscripts will be exceeded.
2050	ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
2051	ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
2052	ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
2053
2054	k_ppp = k + 3 * ibit17
2055	k_pp = k + 2 * ( 1 - ibit15 )
2056	k_mm = k - 2 * ibit17
2057
2058	w_comp = w(k,j,i) + w(k,j,i-1)
2059	flux_t(k) = w_comp * rho_air_zw(k) * ( &
2060	( 37.0_wp * ibit17 * adv_mom_5 &
2061	+ 7.0_wp * ibit16 * adv_mom_3 &
2062	+ ibit15 * adv_mom_1 &
2063	) * &
2064	( u(k+1,j,i) + u(k,j,i) ) &
2065	- ( 8.0_wp * ibit17 * adv_mom_5 &
2066	+ ibit16 * adv_mom_3 &
2067	) * &
2068	( u(k_pp,j,i) + u(k-1,j,i) ) &
2069	+ ( ibit17 * adv_mom_5 &
2070	) * &
2071	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
2072	)
2073
2074	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
2075	( 10.0_wp * ibit17 * adv_mom_5 &
2076	+ 3.0_wp * ibit16 * adv_mom_3 &
2077	+ ibit15 * adv_mom_1 &
2078	) * &
2079	( u(k+1,j,i) - u(k,j,i) ) &
2080	- ( 5.0_wp * ibit17 * adv_mom_5 &
2081	+ ibit16 * adv_mom_3 &
2082	) * &
2083	( u(k_pp,j,i) - u(k-1,j,i) ) &
2084	+ ( ibit17 * adv_mom_5 &
2085	) * &
2086	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
2087	)
2088	!
2089	!-- Calculate the divergence of the velocity field. A respective
2090	!-- correction is needed to overcome numerical instabilities introduced
2091	!-- by a not sufficient reduction of divergences near topography.
2092	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
2093	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
2094	+ ( w_comp * rho_air_zw(k) - &
2095	( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) &
2096	) * drho_air(k) * ddzw(k) &
2097	) * 0.5_wp
2098
2099	tend(k,j,i) = tend(k,j,i) - ( &
2100	( flux_r(k) + diss_r(k) &
2101	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
2102	+ ( flux_n(k) + diss_n(k) &
2103	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &
2104	+ ( ( flux_t(k) + diss_t(k) ) &
2105	- ( flux_d + diss_d ) &
2106	) * drho_air(k) * ddzw(k) &
2107	) + div * u(k,j,i)
2108
2109	flux_l_u(k,j,tn) = flux_r(k)
2110	diss_l_u(k,j,tn) = diss_r(k)
2111	flux_s_u(k,tn) = flux_n(k)
2112	diss_s_u(k,tn) = diss_n(k)
2113	flux_d = flux_t(k)
2114	diss_d = diss_t(k)
2115	!
2116	!-- Statistical Evaluation of u'u'. The factor has to be applied for
2117	!-- right evaluation when gallilei_trans = .T. .
2118	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
2119	+ ( flux_r(k) &
2120	* ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2121	/ ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &
2122	+ diss_r(k) &
2123	* ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2124	/ ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &
2125	) * weight_substep(intermediate_timestep_count)
2126	!
2127	!-- Statistical Evaluation of w'u'.
2128	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
2129	+ ( flux_t(k) &
2130	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2131	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
2132	+ diss_t(k) &
2133	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2134	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
2135	) * weight_substep(intermediate_timestep_count)
2136	ENDDO
2137
2138	sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
2139
2140
2141
2142	END SUBROUTINE advec_u_ws_ij
2143
2144
2145
2146	!-----------------------------------------------------------------------------!
2147	! Description:
2148	! ------------
2149	!> Advection of v-component - Call for grid point i,j
2150	!-----------------------------------------------------------------------------!
2151	SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn )
2152
2153	USE arrays_3d, &
2154	ONLY: ddzw, diss_l_v, diss_s_v, flux_l_v, flux_s_v, tend, u, v, w, &
2155	drho_air, rho_air_zw
2156
2157	USE constants, &
2158	ONLY: adv_mom_1, adv_mom_3, adv_mom_5
2159
2160	USE control_parameters, &
2161	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
2162
2163	USE grid_variables, &
2164	ONLY: ddx, ddy
2165
2166	USE indices, &
2167	ONLY: nxl, nxr, nyn, nys, nysv, nzb, nzb_max, nzt, advc_flags_1
2168
2169	USE kinds
2170
2171	USE statistics, &
2172	ONLY: hom, sums_vs2_ws_l, sums_wsvs_ws_l, weight_substep
2173
2174	IMPLICIT NONE
2175
2176	INTEGER(iwp) :: i !<
2177	INTEGER(iwp) :: ibit18 !<
2178	INTEGER(iwp) :: ibit19 !<
2179	INTEGER(iwp) :: ibit20 !<
2180	INTEGER(iwp) :: ibit21 !<
2181	INTEGER(iwp) :: ibit22 !<
2182	INTEGER(iwp) :: ibit23 !<
2183	INTEGER(iwp) :: ibit24 !<
2184	INTEGER(iwp) :: ibit25 !<
2185	INTEGER(iwp) :: ibit26 !<
2186	INTEGER(iwp) :: i_omp !<
2187	INTEGER(iwp) :: j !<
2188	INTEGER(iwp) :: k !<
2189	INTEGER(iwp) :: k_mm !<
2190	INTEGER(iwp) :: k_pp !<
2191	INTEGER(iwp) :: k_ppp !<
2192	INTEGER(iwp) :: tn !<
2193
2194	REAL(wp) :: diss_d !<
2195	REAL(wp) :: div !<
2196	REAL(wp) :: flux_d !<
2197	REAL(wp) :: gu !<
2198	REAL(wp) :: gv !<
2199	REAL(wp) :: u_comp !<
2200	REAL(wp) :: v_comp_l !<
2201	REAL(wp) :: w_comp !<
2202
2203	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !<
2204	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !<
2205	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !<
2206	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !<
2207	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !<
2208	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !<
2209	REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !<
2210
2211	gu = 2.0_wp * u_gtrans
2212	gv = 2.0_wp * v_gtrans
2213
2214	!
2215	!-- Compute leftside fluxes for the respective boundary.
2216	IF ( i == i_omp ) THEN
2217
2218	DO k = nzb+1, nzb_max
2219
2220	ibit20 = IBITS(advc_flags_1(k,j,i-1),20,1)
2221	ibit19 = IBITS(advc_flags_1(k,j,i-1),19,1)
2222	ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)
2223
2224	u_comp = u(k,j-1,i) + u(k,j,i) - gu
2225	flux_l_v(k,j,tn) = u_comp * ( &
2226	( 37.0_wp * ibit20 * adv_mom_5 &
2227	+ 7.0_wp * ibit19 * adv_mom_3 &
2228	+ ibit18 * adv_mom_1 &
2229	) * &
2230	( v(k,j,i) + v(k,j,i-1) ) &
2231	- ( 8.0_wp * ibit20 * adv_mom_5 &
2232	+ ibit19 * adv_mom_3 &
2233	) * &
2234	( v(k,j,i+1) + v(k,j,i-2) ) &
2235	+ ( ibit20 * adv_mom_5 &
2236	) * &
2237	( v(k,j,i+2) + v(k,j,i-3) ) &
2238	)
2239
2240	diss_l_v(k,j,tn) = - ABS( u_comp ) * ( &
2241	( 10.0_wp * ibit20 * adv_mom_5 &
2242	+ 3.0_wp * ibit19 * adv_mom_3 &
2243	+ ibit18 * adv_mom_1 &
2244	) * &
2245	( v(k,j,i) - v(k,j,i-1) ) &
2246	- ( 5.0_wp * ibit20 * adv_mom_5 &
2247	+ ibit19 * adv_mom_3 &
2248	) * &
2249	( v(k,j,i+1) - v(k,j,i-2) ) &
2250	+ ( ibit20 * adv_mom_5 &
2251	) * &
2252	( v(k,j,i+2) - v(k,j,i-3) ) &
2253	)
2254
2255	ENDDO
2256
2257	DO k = nzb_max+1, nzt
2258
2259	u_comp = u(k,j-1,i) + u(k,j,i) - gu
2260	flux_l_v(k,j,tn) = u_comp * ( &
2261	37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) &
2262	- 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) &
2263	+ ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
2264	diss_l_v(k,j,tn) = - ABS( u_comp ) * ( &
2265	10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) &
2266	- 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) &
2267	+ ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
2268
2269	ENDDO
2270
2271	ENDIF
2272	!
2273	!-- Compute southside fluxes for the respective boundary.
2274	IF ( j == nysv ) THEN
2275
2276	DO k = nzb+1, nzb_max
2277
2278	ibit23 = IBITS(advc_flags_1(k,j-1,i),23,1)
2279	ibit22 = IBITS(advc_flags_1(k,j-1,i),22,1)
2280	ibit21 = IBITS(advc_flags_1(k,j-1,i),21,1)
2281
2282	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
2283	flux_s_v(k,tn) = v_comp_l * ( &
2284	( 37.0_wp * ibit23 * adv_mom_5 &
2285	+ 7.0_wp * ibit22 * adv_mom_3 &
2286	+ ibit21 * adv_mom_1 &
2287	) * &
2288	( v(k,j,i) + v(k,j-1,i) ) &
2289	- ( 8.0_wp * ibit23 * adv_mom_5 &
2290	+ ibit22 * adv_mom_3 &
2291	) * &
2292	( v(k,j+1,i) + v(k,j-2,i) ) &
2293	+ ( ibit23 * adv_mom_5 &
2294	) * &
2295	( v(k,j+2,i) + v(k,j-3,i) ) &
2296	)
2297
2298	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
2299	( 10.0_wp * ibit23 * adv_mom_5 &
2300	+ 3.0_wp * ibit22 * adv_mom_3 &
2301	+ ibit21 * adv_mom_1 &
2302	) * &
2303	( v(k,j,i) - v(k,j-1,i) ) &
2304	- ( 5.0_wp * ibit23 * adv_mom_5 &
2305	+ ibit22 * adv_mom_3 &
2306	) * &
2307	( v(k,j+1,i) - v(k,j-2,i) ) &
2308	+ ( ibit23 * adv_mom_5 &
2309	) * &
2310	( v(k,j+2,i) - v(k,j-3,i) ) &
2311	)
2312
2313	ENDDO
2314
2315	DO k = nzb_max+1, nzt
2316
2317	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
2318	flux_s_v(k,tn) = v_comp_l * ( &
2319	37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) &
2320	- 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) &
2321	+ ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
2322	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
2323	10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) &
2324	- 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) &
2325	+ ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
2326
2327	ENDDO
2328
2329	ENDIF
2330
2331	flux_t(0) = 0.0_wp
2332	diss_t(0) = 0.0_wp
2333	flux_d = 0.0_wp
2334	diss_d = 0.0_wp
2335	!
2336	!-- Now compute the fluxes and tendency terms for the horizontal and
2337	!-- verical parts.
2338	DO k = nzb+1, nzb_max
2339
2340	ibit20 = IBITS(advc_flags_1(k,j,i),20,1)
2341	ibit19 = IBITS(advc_flags_1(k,j,i),19,1)
2342	ibit18 = IBITS(advc_flags_1(k,j,i),18,1)
2343
2344	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
2345	flux_r(k) = u_comp * ( &
2346	( 37.0_wp * ibit20 * adv_mom_5 &
2347	+ 7.0_wp * ibit19 * adv_mom_3 &
2348	+ ibit18 * adv_mom_1 &
2349	) * &
2350	( v(k,j,i+1) + v(k,j,i) ) &
2351	- ( 8.0_wp * ibit20 * adv_mom_5 &
2352	+ ibit19 * adv_mom_3 &
2353	) * &
2354	( v(k,j,i+2) + v(k,j,i-1) ) &
2355	+ ( ibit20 * adv_mom_5 &
2356	) * &
2357	( v(k,j,i+3) + v(k,j,i-2) ) &
2358	)
2359
2360	diss_r(k) = - ABS( u_comp ) * ( &
2361	( 10.0_wp * ibit20 * adv_mom_5 &
2362	+ 3.0_wp * ibit19 * adv_mom_3 &
2363	+ ibit18 * adv_mom_1 &
2364	) * &
2365	( v(k,j,i+1) - v(k,j,i) ) &
2366	- ( 5.0_wp * ibit20 * adv_mom_5 &
2367	+ ibit19 * adv_mom_3 &
2368	) * &
2369	( v(k,j,i+2) - v(k,j,i-1) ) &
2370	+ ( ibit20 * adv_mom_5 &
2371	) * &
2372	( v(k,j,i+3) - v(k,j,i-2) ) &
2373	)
2374
2375	ibit23 = IBITS(advc_flags_1(k,j,i),23,1)
2376	ibit22 = IBITS(advc_flags_1(k,j,i),22,1)
2377	ibit21 = IBITS(advc_flags_1(k,j,i),21,1)
2378
2379
2380	v_comp(k) = v(k,j+1,i) + v(k,j,i)
2381	flux_n(k) = ( v_comp(k) - gv ) * ( &
2382	( 37.0_wp * ibit23 * adv_mom_5 &
2383	+ 7.0_wp * ibit22 * adv_mom_3 &
2384	+ ibit21 * adv_mom_1 &
2385	) * &
2386	( v(k,j+1,i) + v(k,j,i) ) &
2387	- ( 8.0_wp * ibit23 * adv_mom_5 &
2388	+ ibit22 * adv_mom_3 &
2389	) * &
2390	( v(k,j+2,i) + v(k,j-1,i) ) &
2391	+ ( ibit23 * adv_mom_5 &
2392	) * &
2393	( v(k,j+3,i) + v(k,j-2,i) ) &
2394	)
2395
2396	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
2397	( 10.0_wp * ibit23 * adv_mom_5 &
2398	+ 3.0_wp * ibit22 * adv_mom_3 &
2399	+ ibit21 * adv_mom_1 &
2400	) * &
2401	( v(k,j+1,i) - v(k,j,i) ) &
2402	- ( 5.0_wp * ibit23 * adv_mom_5 &
2403	+ ibit22 * adv_mom_3 &
2404	) * &
2405	( v(k,j+2,i) - v(k,j-1,i) ) &
2406	+ ( ibit23 * adv_mom_5 &
2407	) * &
2408	( v(k,j+3,i) - v(k,j-2,i) ) &
2409	)
2410	!
2411	!-- k index has to be modified near bottom and top, else array
2412	!-- subscripts will be exceeded.
2413	ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
2414	ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
2415	ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
2416
2417	k_ppp = k + 3 * ibit26
2418	k_pp = k + 2 * ( 1 - ibit24 )
2419	k_mm = k - 2 * ibit26
2420
2421	w_comp = w(k,j-1,i) + w(k,j,i)
2422	flux_t(k) = w_comp * rho_air_zw(k) * ( &
2423	( 37.0_wp * ibit26 * adv_mom_5 &
2424	+ 7.0_wp * ibit25 * adv_mom_3 &
2425	+ ibit24 * adv_mom_1 &
2426	) * &
2427	( v(k+1,j,i) + v(k,j,i) ) &
2428	- ( 8.0_wp * ibit26 * adv_mom_5 &
2429	+ ibit25 * adv_mom_3 &
2430	) * &
2431	( v(k_pp,j,i) + v(k-1,j,i) ) &
2432	+ ( ibit26 * adv_mom_5 &
2433	) * &
2434	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
2435	)
2436
2437	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
2438	( 10.0_wp * ibit26 * adv_mom_5 &
2439	+ 3.0_wp * ibit25 * adv_mom_3 &
2440	+ ibit24 * adv_mom_1 &
2441	) * &
2442	( v(k+1,j,i) - v(k,j,i) ) &
2443	- ( 5.0_wp * ibit26 * adv_mom_5 &
2444	+ ibit25 * adv_mom_3 &
2445	) * &
2446	( v(k_pp,j,i) - v(k-1,j,i) ) &
2447	+ ( ibit26 * adv_mom_5 &
2448	) * &
2449	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
2450	)
2451	!
2452	!-- Calculate the divergence of the velocity field. A respective
2453	!-- correction is needed to overcome numerical instabilities introduced
2454	!-- by a not sufficient reduction of divergences near topography.
2455	div = ( ( ( u_comp + gu ) &
2456	* ( ibit18 + ibit19 + ibit20 ) &
2457	- ( u(k,j-1,i) + u(k,j,i) ) &
2458	* ( IBITS(advc_flags_1(k,j,i-1),18,1) &
2459	+ IBITS(advc_flags_1(k,j,i-1),19,1) &
2460	+ IBITS(advc_flags_1(k,j,i-1),20,1) &
2461	) &
2462	) * ddx &
2463	+ ( v_comp(k) &
2464	* ( ibit21 + ibit22 + ibit23 ) &
2465	- ( v(k,j,i) + v(k,j-1,i) ) &
2466	* ( IBITS(advc_flags_1(k,j-1,i),21,1) &
2467	+ IBITS(advc_flags_1(k,j-1,i),22,1) &
2468	+ IBITS(advc_flags_1(k,j-1,i),23,1) &
2469	) &
2470	) * ddy &
2471	+ ( w_comp * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 ) &
2472	- ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) &
2473	* ( IBITS(advc_flags_1(k-1,j,i),24,1) &
2474	+ IBITS(advc_flags_1(k-1,j,i),25,1) &
2475	+ IBITS(advc_flags_1(k-1,j,i),26,1) &
2476	) &
2477	) * drho_air(k) * ddzw(k) &
2478	) * 0.5_wp
2479
2480
2481	tend(k,j,i) = tend(k,j,i) - ( &
2482	( flux_r(k) + diss_r(k) &
2483	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
2484	+ ( flux_n(k) + diss_n(k) &
2485	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &
2486	+ ( ( flux_t(k) + diss_t(k) ) &
2487	- ( flux_d + diss_d ) &
2488	) * drho_air(k) * ddzw(k) &
2489	) + v(k,j,i) * div
2490
2491	flux_l_v(k,j,tn) = flux_r(k)
2492	diss_l_v(k,j,tn) = diss_r(k)
2493	flux_s_v(k,tn) = flux_n(k)
2494	diss_s_v(k,tn) = diss_n(k)
2495	flux_d = flux_t(k)
2496	diss_d = diss_t(k)
2497
2498	!
2499	!-- Statistical Evaluation of v'v'. The factor has to be applied for
2500	!-- right evaluation when gallilei_trans = .T. .
2501	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
2502	+ ( flux_n(k) &
2503	* ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2504	/ ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &
2505	+ diss_n(k) &
2506	* ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2507	/ ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &
2508	) * weight_substep(intermediate_timestep_count)
2509	!
2510	!-- Statistical Evaluation of w'u'.
2511	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
2512	+ ( flux_t(k) &
2513	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2514	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
2515	+ diss_t(k) &
2516	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2517	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
2518	) * weight_substep(intermediate_timestep_count)
2519
2520	ENDDO
2521
2522	DO k = nzb_max+1, nzt
2523
2524	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
2525	flux_r(k) = u_comp * ( &
2526	37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) &
2527	- 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) &
2528	+ ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
2529
2530	diss_r(k) = - ABS( u_comp ) * ( &
2531	10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) &
2532	- 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) &
2533	+ ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
2534
2535
2536	v_comp(k) = v(k,j+1,i) + v(k,j,i)
2537	flux_n(k) = ( v_comp(k) - gv ) * ( &
2538	37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) &
2539	- 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) &
2540	+ ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
2541
2542	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
2543	10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) &
2544	- 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) &
2545	+ ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
2546	!
2547	!-- k index has to be modified near bottom and top, else array
2548	!-- subscripts will be exceeded.
2549	ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
2550	ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
2551	ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
2552
2553	k_ppp = k + 3 * ibit26
2554	k_pp = k + 2 * ( 1 - ibit24 )
2555	k_mm = k - 2 * ibit26
2556
2557	w_comp = w(k,j-1,i) + w(k,j,i)
2558	flux_t(k) = w_comp * rho_air_zw(k) * ( &
2559	( 37.0_wp * ibit26 * adv_mom_5 &
2560	+ 7.0_wp * ibit25 * adv_mom_3 &
2561	+ ibit24 * adv_mom_1 &
2562	) * &
2563	( v(k+1,j,i) + v(k,j,i) ) &
2564	- ( 8.0_wp * ibit26 * adv_mom_5 &
2565	+ ibit25 * adv_mom_3 &
2566	) * &
2567	( v(k_pp,j,i) + v(k-1,j,i) ) &
2568	+ ( ibit26 * adv_mom_5 &
2569	) * &
2570	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
2571	)
2572
2573	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
2574	( 10.0_wp * ibit26 * adv_mom_5 &
2575	+ 3.0_wp * ibit25 * adv_mom_3 &
2576	+ ibit24 * adv_mom_1 &
2577	) * &
2578	( v(k+1,j,i) - v(k,j,i) ) &
2579	- ( 5.0_wp * ibit26 * adv_mom_5 &
2580	+ ibit25 * adv_mom_3 &
2581	) * &
2582	( v(k_pp,j,i) - v(k-1,j,i) ) &
2583	+ ( ibit26 * adv_mom_5 &
2584	) * &
2585	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
2586	)
2587	!
2588	!-- Calculate the divergence of the velocity field. A respective
2589	!-- correction is needed to overcome numerical instabilities introduced
2590	!-- by a not sufficient reduction of divergences near topography.
2591	div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
2592	+ ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
2593	+ ( w_comp * rho_air_zw(k) - &
2594	( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) &
2595	) * drho_air(k) * ddzw(k) &
2596	) * 0.5_wp
2597
2598	tend(k,j,i) = tend(k,j,i) - ( &
2599	( flux_r(k) + diss_r(k) &
2600	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
2601	+ ( flux_n(k) + diss_n(k) &
2602	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &
2603	+ ( ( flux_t(k) + diss_t(k) ) &
2604	- ( flux_d + diss_d ) &
2605	) * drho_air(k) * ddzw(k) &
2606	) + v(k,j,i) * div
2607
2608	flux_l_v(k,j,tn) = flux_r(k)
2609	diss_l_v(k,j,tn) = diss_r(k)
2610	flux_s_v(k,tn) = flux_n(k)
2611	diss_s_v(k,tn) = diss_n(k)
2612	flux_d = flux_t(k)
2613	diss_d = diss_t(k)
2614
2615	!
2616	!-- Statistical Evaluation of v'v'. The factor has to be applied for
2617	!-- right evaluation when gallilei_trans = .T. .
2618	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
2619	+ ( flux_n(k) &
2620	* ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2621	/ ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &
2622	+ diss_n(k) &
2623	* ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2624	/ ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &
2625	) * weight_substep(intermediate_timestep_count)
2626	!
2627	!-- Statistical Evaluation of w'u'.
2628	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
2629	+ ( flux_t(k) &
2630	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2631	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
2632	+ diss_t(k) &
2633	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2634	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
2635	) * weight_substep(intermediate_timestep_count)
2636
2637	ENDDO
2638	sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
2639
2640
2641	END SUBROUTINE advec_v_ws_ij
2642
2643
2644
2645	!------------------------------------------------------------------------------!
2646	! Description:
2647	! ------------
2648	!> Advection of w-component - Call for grid point i,j
2649	!------------------------------------------------------------------------------!
2650	SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn )
2651
2652	USE arrays_3d, &
2653	ONLY: ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w,&
2654	drho_air_zw, rho_air
2655
2656	USE constants, &
2657	ONLY: adv_mom_1, adv_mom_3, adv_mom_5
2658
2659	USE control_parameters, &
2660	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
2661
2662	USE grid_variables, &
2663	ONLY: ddx, ddy
2664
2665	USE indices, &
2666	ONLY: nxl, nxr, nyn, nys, nzb, nzb_max, nzt, advc_flags_1, &
2667	advc_flags_2
2668
2669	USE kinds
2670
2671	USE statistics, &
2672	ONLY: hom, sums_ws2_ws_l, weight_substep
2673
2674	IMPLICIT NONE
2675
2676	INTEGER(iwp) :: i !<
2677	INTEGER(iwp) :: ibit27 !<
2678	INTEGER(iwp) :: ibit28 !<
2679	INTEGER(iwp) :: ibit29 !<
2680	INTEGER(iwp) :: ibit30 !<
2681	INTEGER(iwp) :: ibit31 !<
2682	INTEGER(iwp) :: ibit32 !<
2683	INTEGER(iwp) :: ibit33 !<
2684	INTEGER(iwp) :: ibit34 !<
2685	INTEGER(iwp) :: ibit35 !<
2686	INTEGER(iwp) :: i_omp !<
2687	INTEGER(iwp) :: j !<
2688	INTEGER(iwp) :: k !<
2689	INTEGER(iwp) :: k_mm !<
2690	INTEGER(iwp) :: k_pp !<
2691	INTEGER(iwp) :: k_ppp !<
2692	INTEGER(iwp) :: tn !<
2693
2694	REAL(wp) :: diss_d !<
2695	REAL(wp) :: div !<
2696	REAL(wp) :: flux_d !<
2697	REAL(wp) :: gu !<
2698	REAL(wp) :: gv !<
2699	REAL(wp) :: u_comp !<
2700	REAL(wp) :: v_comp !<
2701	REAL(wp) :: w_comp !<
2702
2703	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !<
2704	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !<
2705	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !<
2706	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !<
2707	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !<
2708	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !<
2709
2710	gu = 2.0_wp * u_gtrans
2711	gv = 2.0_wp * v_gtrans
2712
2713	!
2714	!-- Compute southside fluxes for the respective boundary.
2715	IF ( j == nys ) THEN
2716
2717	DO k = nzb+1, nzb_max
2718	ibit32 = IBITS(advc_flags_2(k,j-1,i),0,1)
2719	ibit31 = IBITS(advc_flags_1(k,j-1,i),31,1)
2720	ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)
2721
2722	v_comp = v(k+1,j,i) + v(k,j,i) - gv
2723	flux_s_w(k,tn) = v_comp * ( &
2724	( 37.0_wp * ibit32 * adv_mom_5 &
2725	+ 7.0_wp * ibit31 * adv_mom_3 &
2726	+ ibit30 * adv_mom_1 &
2727	) * &
2728	( w(k,j,i) + w(k,j-1,i) ) &
2729	- ( 8.0_wp * ibit32 * adv_mom_5 &
2730	+ ibit31 * adv_mom_3 &
2731	) * &
2732	( w(k,j+1,i) + w(k,j-2,i) ) &
2733	+ ( ibit32 * adv_mom_5 &
2734	) * &
2735	( w(k,j+2,i) + w(k,j-3,i) ) &
2736	)
2737
2738	diss_s_w(k,tn) = - ABS( v_comp ) * ( &
2739	( 10.0_wp * ibit32 * adv_mom_5 &
2740	+ 3.0_wp * ibit31 * adv_mom_3 &
2741	+ ibit30 * adv_mom_1 &
2742	) * &
2743	( w(k,j,i) - w(k,j-1,i) ) &
2744	- ( 5.0_wp * ibit32 * adv_mom_5 &
2745	+ ibit31 * adv_mom_3 &
2746	) * &
2747	( w(k,j+1,i) - w(k,j-2,i) ) &
2748	+ ( ibit32 * adv_mom_5 &
2749	) * &
2750	( w(k,j+2,i) - w(k,j-3,i) ) &
2751	)
2752
2753	ENDDO
2754
2755	DO k = nzb_max+1, nzt
2756
2757	v_comp = v(k+1,j,i) + v(k,j,i) - gv
2758	flux_s_w(k,tn) = v_comp * ( &
2759	37.0_wp * ( w(k,j,i) + w(k,j-1,i) ) &
2760	- 8.0_wp * ( w(k,j+1,i) +w(k,j-2,i) ) &
2761	+ ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
2762	diss_s_w(k,tn) = - ABS( v_comp ) * ( &
2763	10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) &
2764	- 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) &
2765	+ ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
2766
2767	ENDDO
2768
2769	ENDIF
2770	!
2771	!-- Compute leftside fluxes for the respective boundary.
2772	IF ( i == i_omp ) THEN
2773
2774	DO k = nzb+1, nzb_max
2775
2776	ibit29 = IBITS(advc_flags_1(k,j,i-1),29,1)
2777	ibit28 = IBITS(advc_flags_1(k,j,i-1),28,1)
2778	ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)
2779
2780	u_comp = u(k+1,j,i) + u(k,j,i) - gu
2781	flux_l_w(k,j,tn) = u_comp * ( &
2782	( 37.0_wp * ibit29 * adv_mom_5 &
2783	+ 7.0_wp * ibit28 * adv_mom_3 &
2784	+ ibit27 * adv_mom_1 &
2785	) * &
2786	( w(k,j,i) + w(k,j,i-1) ) &
2787	- ( 8.0_wp * ibit29 * adv_mom_5 &
2788	+ ibit28 * adv_mom_3 &
2789	) * &
2790	( w(k,j,i+1) + w(k,j,i-2) ) &
2791	+ ( ibit29 * adv_mom_5 &
2792	) * &
2793	( w(k,j,i+2) + w(k,j,i-3) ) &
2794	)
2795
2796	diss_l_w(k,j,tn) = - ABS( u_comp ) * ( &
2797	( 10.0_wp * ibit29 * adv_mom_5 &
2798	+ 3.0_wp * ibit28 * adv_mom_3 &
2799	+ ibit27 * adv_mom_1 &
2800	) * &
2801	( w(k,j,i) - w(k,j,i-1) ) &
2802	- ( 5.0_wp * ibit29 * adv_mom_5 &
2803	+ ibit28 * adv_mom_3 &
2804	) * &
2805	( w(k,j,i+1) - w(k,j,i-2) ) &
2806	+ ( ibit29 * adv_mom_5 &
2807	) * &
2808	( w(k,j,i+2) - w(k,j,i-3) ) &
2809	)
2810
2811	ENDDO
2812
2813	DO k = nzb_max+1, nzt
2814
2815	u_comp = u(k+1,j,i) + u(k,j,i) - gu
2816	flux_l_w(k,j,tn) = u_comp * ( &
2817	37.0_wp * ( w(k,j,i) + w(k,j,i-1) ) &
2818	- 8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) ) &
2819	+ ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
2820	diss_l_w(k,j,tn) = - ABS( u_comp ) * ( &
2821	10.0_wp * ( w(k,j,i) - w(k,j,i-1) ) &
2822	- 5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) ) &
2823	+ ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
2824
2825	ENDDO
2826
2827	ENDIF
2828	!
2829	!-- The lower flux has to be calculated explicetely for the tendency at
2830	!-- the first w-level. For topography wall this is done implicitely by
2831	!-- advc_flags_1.
2832	k = nzb + 1
2833	w_comp = w(k,j,i) + w(k-1,j,i)
2834	flux_t(0) = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
2835	diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
2836	flux_d = flux_t(0)
2837	diss_d = diss_t(0)
2838	!
2839	!-- Now compute the fluxes and tendency terms for the horizontal
2840	!-- and vertical parts.
2841	DO k = nzb+1, nzb_max
2842
2843	ibit29 = IBITS(advc_flags_1(k,j,i),29,1)
2844	ibit28 = IBITS(advc_flags_1(k,j,i),28,1)
2845	ibit27 = IBITS(advc_flags_1(k,j,i),27,1)
2846
2847	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
2848	flux_r(k) = u_comp * ( &
2849	( 37.0_wp * ibit29 * adv_mom_5 &
2850	+ 7.0_wp * ibit28 * adv_mom_3 &
2851	+ ibit27 * adv_mom_1 &
2852	) * &
2853	( w(k,j,i+1) + w(k,j,i) ) &
2854	- ( 8.0_wp * ibit29 * adv_mom_5 &
2855	+ ibit28 * adv_mom_3 &
2856	) * &
2857	( w(k,j,i+2) + w(k,j,i-1) ) &
2858	+ ( ibit29 * adv_mom_5 &
2859	) * &
2860	( w(k,j,i+3) + w(k,j,i-2) ) &
2861	)
2862
2863	diss_r(k) = - ABS( u_comp ) * ( &
2864	( 10.0_wp * ibit29 * adv_mom_5 &
2865	+ 3.0_wp * ibit28 * adv_mom_3 &
2866	+ ibit27 * adv_mom_1 &
2867	) * &
2868	( w(k,j,i+1) - w(k,j,i) ) &
2869	- ( 5.0_wp * ibit29 * adv_mom_5 &
2870	+ ibit28 * adv_mom_3 &
2871	) * &
2872	( w(k,j,i+2) - w(k,j,i-1) ) &
2873	+ ( ibit29 * adv_mom_5 &
2874	) * &
2875	( w(k,j,i+3) - w(k,j,i-2) ) &
2876	)
2877
2878	ibit32 = IBITS(advc_flags_2(k,j,i),0,1)
2879	ibit31 = IBITS(advc_flags_1(k,j,i),31,1)
2880	ibit30 = IBITS(advc_flags_1(k,j,i),30,1)
2881
2882	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
2883	flux_n(k) = v_comp * ( &
2884	( 37.0_wp * ibit32 * adv_mom_5 &
2885	+ 7.0_wp * ibit31 * adv_mom_3 &
2886	+ ibit30 * adv_mom_1 &
2887	) * &
2888	( w(k,j+1,i) + w(k,j,i) ) &
2889	- ( 8.0_wp * ibit32 * adv_mom_5 &
2890	+ ibit31 * adv_mom_3 &
2891	) * &
2892	( w(k,j+2,i) + w(k,j-1,i) ) &
2893	+ ( ibit32 * adv_mom_5 &
2894	) * &
2895	( w(k,j+3,i) + w(k,j-2,i) ) &
2896	)
2897
2898	diss_n(k) = - ABS( v_comp ) * ( &
2899	( 10.0_wp * ibit32 * adv_mom_5 &
2900	+ 3.0_wp * ibit31 * adv_mom_3 &
2901	+ ibit30 * adv_mom_1 &
2902	) * &
2903	( w(k,j+1,i) - w(k,j,i) ) &
2904	- ( 5.0_wp * ibit32 * adv_mom_5 &
2905	+ ibit31 * adv_mom_3 &
2906	) * &
2907	( w(k,j+2,i) - w(k,j-1,i) ) &
2908	+ ( ibit32 * adv_mom_5 &
2909	) * &
2910	( w(k,j+3,i) - w(k,j-2,i) ) &
2911	)
2912	!
2913	!-- k index has to be modified near bottom and top, else array
2914	!-- subscripts will be exceeded.
2915	ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
2916	ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
2917	ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
2918
2919	k_ppp = k + 3 * ibit35
2920	k_pp = k + 2 * ( 1 - ibit33 )
2921	k_mm = k - 2 * ibit35
2922
2923	w_comp = w(k+1,j,i) + w(k,j,i)
2924	flux_t(k) = w_comp * rho_air(k+1) * ( &
2925	( 37.0_wp * ibit35 * adv_mom_5 &
2926	+ 7.0_wp * ibit34 * adv_mom_3 &
2927	+ ibit33 * adv_mom_1 &
2928	) * &
2929	( w(k+1,j,i) + w(k,j,i) ) &
2930	- ( 8.0_wp * ibit35 * adv_mom_5 &
2931	+ ibit34 * adv_mom_3 &
2932	) * &
2933	( w(k_pp,j,i) + w(k-1,j,i) ) &
2934	+ ( ibit35 * adv_mom_5 &
2935	) * &
2936	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
2937	)
2938
2939	diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * ( &
2940	( 10.0_wp * ibit35 * adv_mom_5 &
2941	+ 3.0_wp * ibit34 * adv_mom_3 &
2942	+ ibit33 * adv_mom_1 &
2943	) * &
2944	( w(k+1,j,i) - w(k,j,i) ) &
2945	- ( 5.0_wp * ibit35 * adv_mom_5 &
2946	+ ibit34 * adv_mom_3 &
2947	) * &
2948	( w(k_pp,j,i) - w(k-1,j,i) ) &
2949	+ ( ibit35 * adv_mom_5 &
2950	) * &
2951	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
2952	)
2953
2954	!
2955	!-- Calculate the divergence of the velocity field. A respective
2956	!-- correction is needed to overcome numerical instabilities introduced
2957	!-- by a not sufficient reduction of divergences near topography.
2958	div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 ) &
2959	- ( u(k+1,j,i) + u(k,j,i) ) &
2960	* ( IBITS(advc_flags_1(k,j,i-1),27,1) &
2961	+ IBITS(advc_flags_1(k,j,i-1),28,1) &
2962	+ IBITS(advc_flags_1(k,j,i-1),29,1) &
2963	) &
2964	) * ddx &
2965	+ ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) &
2966	- ( v(k+1,j,i) + v(k,j,i) ) &
2967	* ( IBITS(advc_flags_1(k,j-1,i),30,1) &
2968	+ IBITS(advc_flags_1(k,j-1,i),31,1) &
2969	+ IBITS(advc_flags_2(k,j-1,i),0,1) &
2970	) &
2971	) * ddy &
2972	+ ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 ) &
2973	- ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) &
2974	* ( IBITS(advc_flags_2(k-1,j,i),1,1) &
2975	+ IBITS(advc_flags_2(k-1,j,i),2,1) &
2976	+ IBITS(advc_flags_2(k-1,j,i),3,1) &
2977	) &
2978	) * drho_air_zw(k) * ddzu(k+1) &
2979	) * 0.5_wp
2980
2981
2982	tend(k,j,i) = tend(k,j,i) - ( &
2983	( flux_r(k) + diss_r(k) &
2984	- flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &
2985	+ ( flux_n(k) + diss_n(k) &
2986	- flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy &
2987	+ ( ( flux_t(k) + diss_t(k) ) &
2988	- ( flux_d + diss_d ) &
2989	) * drho_air_zw(k) * ddzu(k+1) &
2990	) + div * w(k,j,i)
2991
2992	flux_l_w(k,j,tn) = flux_r(k)
2993	diss_l_w(k,j,tn) = diss_r(k)
2994	flux_s_w(k,tn) = flux_n(k)
2995	diss_s_w(k,tn) = diss_n(k)
2996	flux_d = flux_t(k)
2997	diss_d = diss_t(k)
2998	!
2999	!-- Statistical Evaluation of w'w'.
3000	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
3001	+ ( flux_t(k) &
3002	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
3003	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
3004	+ diss_t(k) &
3005	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
3006	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
3007	) * weight_substep(intermediate_timestep_count)
3008
3009	ENDDO
3010
3011	DO k = nzb_max+1, nzt
3012
3013	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
3014	flux_r(k) = u_comp * ( &
3015	37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) &
3016	- 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) &
3017	+ ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
3018
3019	diss_r(k) = - ABS( u_comp ) * ( &
3020	10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) &
3021	- 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) &
3022	+ ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
3023
3024	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
3025	flux_n(k) = v_comp * ( &
3026	37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) &
3027	- 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) &
3028	+ ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
3029
3030	diss_n(k) = - ABS( v_comp ) * ( &
3031	10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) &
3032	- 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) &
3033	+ ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
3034	!
3035	!-- k index has to be modified near bottom and top, else array
3036	!-- subscripts will be exceeded.
3037	ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
3038	ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
3039	ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
3040
3041	k_ppp = k + 3 * ibit35
3042	k_pp = k + 2 * ( 1 - ibit33 )
3043	k_mm = k - 2 * ibit35
3044
3045	w_comp = w(k+1,j,i) + w(k,j,i)
3046	flux_t(k) = w_comp * rho_air(k+1) * ( &
3047	( 37.0_wp * ibit35 * adv_mom_5 &
3048	+ 7.0_wp * ibit34 * adv_mom_3 &
3049	+ ibit33 * adv_mom_1 &
3050	) * &
3051	( w(k+1,j,i) + w(k,j,i) ) &
3052	- ( 8.0_wp * ibit35 * adv_mom_5 &
3053	+ ibit34 * adv_mom_3 &
3054	) * &
3055	( w(k_pp,j,i) + w(k-1,j,i) ) &
3056	+ ( ibit35 * adv_mom_5 &
3057	) * &
3058	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
3059	)
3060
3061	diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * ( &
3062	( 10.0_wp * ibit35 * adv_mom_5 &
3063	+ 3.0_wp * ibit34 * adv_mom_3 &
3064	+ ibit33 * adv_mom_1 &
3065	) * &
3066	( w(k+1,j,i) - w(k,j,i) ) &
3067	- ( 5.0_wp * ibit35 * adv_mom_5 &
3068	+ ibit34 * adv_mom_3 &
3069	) * &
3070	( w(k_pp,j,i) - w(k-1,j,i) ) &
3071	+ ( ibit35 * adv_mom_5 &
3072	) * &
3073	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
3074	)
3075	!
3076	!-- Calculate the divergence of the velocity field. A respective
3077	!-- correction is needed to overcome numerical instabilities introduced
3078	!-- by a not sufficient reduction of divergences near topography.
3079	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
3080	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
3081	+ ( w_comp * rho_air(k+1) - &
3082	( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) &
3083	) * drho_air_zw(k) * ddzu(k+1) &
3084	) * 0.5_wp
3085
3086	tend(k,j,i) = tend(k,j,i) - ( &
3087	( flux_r(k) + diss_r(k) &
3088	- flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &
3089	+ ( flux_n(k) + diss_n(k) &
3090	- flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy &
3091	+ ( ( flux_t(k) + diss_t(k) ) &
3092	- ( flux_d + diss_d ) &
3093	) * drho_air_zw(k) * ddzu(k+1) &
3094	) + div * w(k,j,i)
3095
3096	flux_l_w(k,j,tn) = flux_r(k)
3097	diss_l_w(k,j,tn) = diss_r(k)
3098	flux_s_w(k,tn) = flux_n(k)
3099	diss_s_w(k,tn) = diss_n(k)
3100	flux_d = flux_t(k)
3101	diss_d = diss_t(k)
3102	!
3103	!-- Statistical Evaluation of w'w'.
3104	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
3105	+ ( flux_t(k) &
3106	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
3107	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
3108	+ diss_t(k) &
3109	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
3110	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
3111	) * weight_substep(intermediate_timestep_count)
3112
3113	ENDDO
3114
3115
3116	END SUBROUTINE advec_w_ws_ij
3117
3118
3119	!------------------------------------------------------------------------------!
3120	! Description:
3121	! ------------
3122	!> Scalar advection - Call for all grid points
3123	!------------------------------------------------------------------------------!
3124	SUBROUTINE advec_s_ws( sk, sk_char )
3125
3126	USE arrays_3d, &
3127	ONLY: ddzw, drho_air, tend, u, v, w, rho_air_zw
3128
3129	USE constants, &
3130	ONLY: adv_sca_1, adv_sca_3, adv_sca_5
3131
3132	USE control_parameters, &
3133	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
3134
3135	USE grid_variables, &
3136	ONLY: ddx, ddy
3137
3138	USE indices, &
3139	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, &
3140	nzt, advc_flags_1
3141
3142	USE kinds
3143
3144	USE statistics, &
3145	ONLY: hom, sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l, &
3146	sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsncs_ws_l, &
3147	sums_wsnrs_ws_l, sums_wsss_ws_l, weight_substep
3148
3149
3150
3151	IMPLICIT NONE
3152
3153	CHARACTER (LEN = *), INTENT(IN) :: sk_char !<
3154
3155	INTEGER(iwp) :: i !<
3156	INTEGER(iwp) :: ibit0 !<
3157	INTEGER(iwp) :: ibit1 !<
3158	INTEGER(iwp) :: ibit2 !<
3159	INTEGER(iwp) :: ibit3 !<
3160	INTEGER(iwp) :: ibit4 !<
3161	INTEGER(iwp) :: ibit5 !<
3162	INTEGER(iwp) :: ibit6 !<
3163	INTEGER(iwp) :: ibit7 !<
3164	INTEGER(iwp) :: ibit8 !<
3165	INTEGER(iwp) :: j !<
3166	INTEGER(iwp) :: k !<
3167	INTEGER(iwp) :: k_mm !<
3168	INTEGER(iwp) :: k_mmm !<
3169	INTEGER(iwp) :: k_pp !<
3170	INTEGER(iwp) :: k_ppp !<
3171	INTEGER(iwp) :: tn = 0 !<
3172
3173	#if defined( __nopointer )
3174	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !<
3175	#else
3176	REAL(wp), DIMENSION(:,:,:), POINTER :: sk !<
3177	#endif
3178
3179	REAL(wp) :: diss_d !<
3180	REAL(wp) :: div !<
3181	REAL(wp) :: flux_d !<
3182	REAL(wp) :: u_comp !<
3183	REAL(wp) :: v_comp !<
3184
3185	REAL(wp), DIMENSION(nzb:nzt) :: diss_n !<
3186	REAL(wp), DIMENSION(nzb:nzt) :: diss_r !<
3187	REAL(wp), DIMENSION(nzb:nzt) :: diss_t !<
3188	REAL(wp), DIMENSION(nzb:nzt) :: flux_n !<
3189	REAL(wp), DIMENSION(nzb:nzt) :: flux_r !<
3190	REAL(wp), DIMENSION(nzb:nzt) :: flux_t !<
3191
3192	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local !<
3193	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local !<
3194
3195	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local !<
3196	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local !<
3197
3198
3199	!
3200	!-- Compute the fluxes for the whole left boundary of the processor domain.
3201	i = nxl
3202	DO j = nys, nyn
3203
3204	DO k = nzb+1, nzb_max
3205
3206	ibit2 = IBITS(advc_flags_1(k,j,i-1),2,1)
3207	ibit1 = IBITS(advc_flags_1(k,j,i-1),1,1)
3208	ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
3209
3210	u_comp = u(k,j,i) - u_gtrans
3211	swap_flux_x_local(k,j) = u_comp * ( &
3212	( 37.0_wp * ibit2 * adv_sca_5 &
3213	+ 7.0_wp * ibit1 * adv_sca_3 &
3214	+ ibit0 * adv_sca_1 &
3215	) * &
3216	( sk(k,j,i) + sk(k,j,i-1) ) &
3217	- ( 8.0_wp * ibit2 * adv_sca_5 &
3218	+ ibit1 * adv_sca_3 &
3219	) * &
3220	( sk(k,j,i+1) + sk(k,j,i-2) ) &
3221	+ ( ibit2 * adv_sca_5 &
3222	) * &
3223	( sk(k,j,i+2) + sk(k,j,i-3) ) &
3224	)
3225
3226	swap_diss_x_local(k,j) = -ABS( u_comp ) * ( &
3227	( 10.0_wp * ibit2 * adv_sca_5 &
3228	+ 3.0_wp * ibit1 * adv_sca_3 &
3229	+ ibit0 * adv_sca_1 &
3230	) * &
3231	( sk(k,j,i) - sk(k,j,i-1) ) &
3232	- ( 5.0_wp * ibit2 * adv_sca_5 &
3233	+ ibit1 * adv_sca_3 &
3234	) * &
3235	( sk(k,j,i+1) - sk(k,j,i-2) ) &
3236	+ ( ibit2 * adv_sca_5 &
3237	) * &
3238	( sk(k,j,i+2) - sk(k,j,i-3) ) &
3239	)
3240
3241	ENDDO
3242
3243	DO k = nzb_max+1, nzt
3244
3245	u_comp = u(k,j,i) - u_gtrans
3246	swap_flux_x_local(k,j) = u_comp * ( &
3247	37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) &
3248	- 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
3249	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) &
3250	) * adv_sca_5
3251
3252	swap_diss_x_local(k,j) = -ABS( u_comp ) * ( &
3253	10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) &
3254	- 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
3255	+ ( sk(k,j,i+2) - sk(k,j,i-3) ) &
3256	) * adv_sca_5
3257
3258	ENDDO
3259
3260	ENDDO
3261
3262	DO i = nxl, nxr
3263
3264	j = nys
3265	DO k = nzb+1, nzb_max
3266
3267	ibit5 = IBITS(advc_flags_1(k,j-1,i),5,1)
3268	ibit4 = IBITS(advc_flags_1(k,j-1,i),4,1)
3269	ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
3270
3271	v_comp = v(k,j,i) - v_gtrans
3272	swap_flux_y_local(k) = v_comp * ( &
3273	( 37.0_wp * ibit5 * adv_sca_5 &
3274	+ 7.0_wp * ibit4 * adv_sca_3 &
3275	+ ibit3 * adv_sca_1 &
3276	) * &
3277	( sk(k,j,i) + sk(k,j-1,i) ) &
3278	- ( 8.0_wp * ibit5 * adv_sca_5 &
3279	+ ibit4 * adv_sca_3 &
3280	) * &
3281	( sk(k,j+1,i) + sk(k,j-2,i) ) &
3282	+ ( ibit5 * adv_sca_5 &
3283	) * &
3284	( sk(k,j+2,i) + sk(k,j-3,i) ) &
3285	)
3286
3287	swap_diss_y_local(k) = -ABS( v_comp ) * ( &
3288	( 10.0_wp * ibit5 * adv_sca_5 &
3289	+ 3.0_wp * ibit4 * adv_sca_3 &
3290	+ ibit3 * adv_sca_1 &
3291	) * &
3292	( sk(k,j,i) - sk(k,j-1,i) ) &
3293	- ( 5.0_wp * ibit5 * adv_sca_5 &
3294	+ ibit4 * adv_sca_3 &
3295	) * &
3296	( sk(k,j+1,i) - sk(k,j-2,i) ) &
3297	+ ( ibit5 * adv_sca_5 &
3298	) * &
3299	( sk(k,j+2,i) - sk(k,j-3,i) ) &
3300	)
3301
3302	ENDDO
3303	!
3304	!-- Above to the top of the highest topography. No degradation necessary.
3305	DO k = nzb_max+1, nzt
3306
3307	v_comp = v(k,j,i) - v_gtrans
3308	swap_flux_y_local(k) = v_comp * ( &
3309	37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) &
3310	- 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
3311	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) &
3312	) * adv_sca_5
3313	swap_diss_y_local(k) = -ABS( v_comp ) * ( &
3314	10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) &
3315	- 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
3316	+ sk(k,j+2,i) - sk(k,j-3,i) &
3317	) * adv_sca_5
3318
3319	ENDDO
3320
3321	DO j = nys, nyn
3322
3323	flux_t(0) = 0.0_wp
3324	diss_t(0) = 0.0_wp
3325	flux_d = 0.0_wp
3326	diss_d = 0.0_wp
3327
3328	DO k = nzb+1, nzb_max
3329
3330	ibit2 = IBITS(advc_flags_1(k,j,i),2,1)
3331	ibit1 = IBITS(advc_flags_1(k,j,i),1,1)
3332	ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
3333
3334	u_comp = u(k,j,i+1) - u_gtrans
3335	flux_r(k) = u_comp * ( &
3336	( 37.0_wp * ibit2 * adv_sca_5 &
3337	+ 7.0_wp * ibit1 * adv_sca_3 &
3338	+ ibit0 * adv_sca_1 &
3339	) * &
3340	( sk(k,j,i+1) + sk(k,j,i) ) &
3341	- ( 8.0_wp * ibit2 * adv_sca_5 &
3342	+ ibit1 * adv_sca_3 &
3343	) * &
3344	( sk(k,j,i+2) + sk(k,j,i-1) ) &
3345	+ ( ibit2 * adv_sca_5 &
3346	) * &
3347	( sk(k,j,i+3) + sk(k,j,i-2) ) &
3348	)
3349
3350	diss_r(k) = -ABS( u_comp ) * ( &
3351	( 10.0_wp * ibit2 * adv_sca_5 &
3352	+ 3.0_wp * ibit1 * adv_sca_3 &
3353	+ ibit0 * adv_sca_1 &
3354	) * &
3355	( sk(k,j,i+1) - sk(k,j,i) ) &
3356	- ( 5.0_wp * ibit2 * adv_sca_5 &
3357	+ ibit1 * adv_sca_3 &
3358	) * &
3359	( sk(k,j,i+2) - sk(k,j,i-1) ) &
3360	+ ( ibit2 * adv_sca_5 &
3361	) * &
3362	( sk(k,j,i+3) - sk(k,j,i-2) ) &
3363	)
3364
3365	ibit5 = IBITS(advc_flags_1(k,j,i),5,1)
3366	ibit4 = IBITS(advc_flags_1(k,j,i),4,1)
3367	ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
3368
3369	v_comp = v(k,j+1,i) - v_gtrans
3370	flux_n(k) = v_comp * ( &
3371	( 37.0_wp * ibit5 * adv_sca_5 &
3372	+ 7.0_wp * ibit4 * adv_sca_3 &
3373	+ ibit3 * adv_sca_1 &
3374	) * &
3375	( sk(k,j+1,i) + sk(k,j,i) ) &
3376	- ( 8.0_wp * ibit5 * adv_sca_5 &
3377	+ ibit4 * adv_sca_3 &
3378	) * &
3379	( sk(k,j+2,i) + sk(k,j-1,i) ) &
3380	+ ( ibit5 * adv_sca_5 &
3381	) * &
3382	( sk(k,j+3,i) + sk(k,j-2,i) ) &
3383	)
3384
3385	diss_n(k) = -ABS( v_comp ) * ( &
3386	( 10.0_wp * ibit5 * adv_sca_5 &
3387	+ 3.0_wp * ibit4 * adv_sca_3 &
3388	+ ibit3 * adv_sca_1 &
3389	) * &
3390	( sk(k,j+1,i) - sk(k,j,i) ) &
3391	- ( 5.0_wp * ibit5 * adv_sca_5 &
3392	+ ibit4 * adv_sca_3 &
3393	) * &
3394	( sk(k,j+2,i) - sk(k,j-1,i) ) &
3395	+ ( ibit5 * adv_sca_5 &
3396	) * &
3397	( sk(k,j+3,i) - sk(k,j-2,i) ) &
3398	)
3399	!
3400	!-- k index has to be modified near bottom and top, else array
3401	!-- subscripts will be exceeded.
3402	ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
3403	ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
3404	ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
3405
3406	k_ppp = k + 3 * ibit8
3407	k_pp = k + 2 * ( 1 - ibit6 )
3408	k_mm = k - 2 * ibit8
3409
3410
3411	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
3412	( 37.0_wp * ibit8 * adv_sca_5 &
3413	+ 7.0_wp * ibit7 * adv_sca_3 &
3414	+ ibit6 * adv_sca_1 &
3415	) * &
3416	( sk(k+1,j,i) + sk(k,j,i) ) &
3417	- ( 8.0_wp * ibit8 * adv_sca_5 &
3418	+ ibit7 * adv_sca_3 &
3419	) * &
3420	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
3421	+ ( ibit8 * adv_sca_5 &
3422	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
3423	)
3424
3425	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
3426	( 10.0_wp * ibit8 * adv_sca_5 &
3427	+ 3.0_wp * ibit7 * adv_sca_3 &
3428	+ ibit6 * adv_sca_1 &
3429	) * &
3430	( sk(k+1,j,i) - sk(k,j,i) ) &
3431	- ( 5.0_wp * ibit8 * adv_sca_5 &
3432	+ ibit7 * adv_sca_3 &
3433	) * &
3434	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
3435	+ ( ibit8 * adv_sca_5 &
3436	) * &
3437	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
3438	)
3439	!
3440	!-- Calculate the divergence of the velocity field. A respective
3441	!-- correction is needed to overcome numerical instabilities caused
3442	!-- by a not sufficient reduction of divergences near topography.
3443	div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) &
3444	- u(k,j,i) * ( IBITS(advc_flags_1(k,j,i-1),0,1) &
3445	+ IBITS(advc_flags_1(k,j,i-1),1,1) &
3446	+ IBITS(advc_flags_1(k,j,i-1),2,1) &
3447	) &
3448	) * ddx &
3449	+ ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) &
3450	- v(k,j,i) * ( IBITS(advc_flags_1(k,j-1,i),3,1) &
3451	+ IBITS(advc_flags_1(k,j-1,i),4,1) &
3452	+ IBITS(advc_flags_1(k,j-1,i),5,1) &
3453	) &
3454	) * ddy &
3455	+ ( w(k,j,i) * rho_air_zw(k) * &
3456	( ibit6 + ibit7 + ibit8 ) &
3457	- w(k-1,j,i) * rho_air_zw(k-1) * &
3458	( IBITS(advc_flags_1(k-1,j,i),6,1) &
3459	+ IBITS(advc_flags_1(k-1,j,i),7,1) &
3460	+ IBITS(advc_flags_1(k-1,j,i),8,1) &
3461	) &
3462	) * drho_air(k) * ddzw(k)
3463
3464
3465	tend(k,j,i) = tend(k,j,i) - ( &
3466	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
3467	swap_diss_x_local(k,j) ) * ddx &
3468	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - &
3469	swap_diss_y_local(k) ) * ddy &
3470	+ ( ( flux_t(k) + diss_t(k) ) - &
3471	( flux_d + diss_d ) &
3472	) * drho_air(k) * ddzw(k) &
3473	) + sk(k,j,i) * div
3474
3475	swap_flux_y_local(k) = flux_n(k)
3476	swap_diss_y_local(k) = diss_n(k)
3477	swap_flux_x_local(k,j) = flux_r(k)
3478	swap_diss_x_local(k,j) = diss_r(k)
3479	flux_d = flux_t(k)
3480	diss_d = diss_t(k)
3481
3482	ENDDO
3483
3484	DO k = nzb_max+1, nzt
3485
3486	u_comp = u(k,j,i+1) - u_gtrans
3487	flux_r(k) = u_comp * ( &
3488	37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) &
3489	- 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
3490	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
3491	diss_r(k) = -ABS( u_comp ) * ( &
3492	10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) &
3493	- 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
3494	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
3495
3496	v_comp = v(k,j+1,i) - v_gtrans
3497	flux_n(k) = v_comp * ( &
3498	37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) &
3499	- 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
3500	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
3501	diss_n(k) = -ABS( v_comp ) * ( &
3502	10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) &
3503	- 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
3504	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
3505	!
3506	!-- k index has to be modified near bottom and top, else array
3507	!-- subscripts will be exceeded.
3508	ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
3509	ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
3510	ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
3511
3512	k_ppp = k + 3 * ibit8
3513	k_pp = k + 2 * ( 1 - ibit6 )
3514	k_mm = k - 2 * ibit8
3515
3516
3517	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
3518	( 37.0_wp * ibit8 * adv_sca_5 &
3519	+ 7.0_wp * ibit7 * adv_sca_3 &
3520	+ ibit6 * adv_sca_1 &
3521	) * &
3522	( sk(k+1,j,i) + sk(k,j,i) ) &
3523	- ( 8.0_wp * ibit8 * adv_sca_5 &
3524	+ ibit7 * adv_sca_3 &
3525	) * &
3526	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
3527	+ ( ibit8 * adv_sca_5 &
3528	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
3529	)
3530
3531	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
3532	( 10.0_wp * ibit8 * adv_sca_5 &
3533	+ 3.0_wp * ibit7 * adv_sca_3 &
3534	+ ibit6 * adv_sca_1 &
3535	) * &
3536	( sk(k+1,j,i) - sk(k,j,i) ) &
3537	- ( 5.0_wp * ibit8 * adv_sca_5 &
3538	+ ibit7 * adv_sca_3 &
3539	) * &
3540	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
3541	+ ( ibit8 * adv_sca_5 &
3542	) * &
3543	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
3544	)
3545	!
3546	!-- Calculate the divergence of the velocity field. A respective
3547	!-- correction is needed to overcome numerical instabilities introduced
3548	!-- by a not sufficient reduction of divergences near topography.
3549	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
3550	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
3551	+ ( w(k,j,i) * rho_air_zw(k) - &
3552	w(k-1,j,i) * rho_air_zw(k-1) &
3553	) * drho_air(k) * ddzw(k)
3554
3555	tend(k,j,i) = tend(k,j,i) - ( &
3556	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
3557	swap_diss_x_local(k,j) ) * ddx &
3558	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - &
3559	swap_diss_y_local(k) ) * ddy &
3560	+ ( ( flux_t(k) + diss_t(k) ) - &
3561	( flux_d + diss_d ) &
3562	) * drho_air(k) * ddzw(k) &
3563	) + sk(k,j,i) * div
3564
3565	swap_flux_y_local(k) = flux_n(k)
3566	swap_diss_y_local(k) = diss_n(k)
3567	swap_flux_x_local(k,j) = flux_r(k)
3568	swap_diss_x_local(k,j) = diss_r(k)
3569	flux_d = flux_t(k)
3570	diss_d = diss_t(k)
3571
3572	ENDDO
3573	!
3574	!-- Evaluation of statistics.
3575	SELECT CASE ( sk_char )
3576
3577	CASE ( 'pt' )
3578	DO k = nzb, nzt
3579	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) &
3580	+ ( flux_t(k) &
3581	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3582	* ( w(k,j,i) - hom(k,1,3,0) ) &
3583	+ diss_t(k) &
3584	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3585	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3586	) * weight_substep(intermediate_timestep_count)
3587	ENDDO
3588	CASE ( 'sa' )
3589	DO k = nzb, nzt
3590	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) &
3591	+ ( flux_t(k) &
3592	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3593	* ( w(k,j,i) - hom(k,1,3,0) ) &
3594	+ diss_t(k) &
3595	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3596	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3597	) * weight_substep(intermediate_timestep_count)
3598	ENDDO
3599	CASE ( 'q' )
3600	DO k = nzb, nzt
3601	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) &
3602	+ ( flux_t(k) &
3603	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3604	* ( w(k,j,i) - hom(k,1,3,0) ) &
3605	+ diss_t(k) &
3606	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3607	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3608	) * weight_substep(intermediate_timestep_count)
3609	ENDDO
3610	CASE ( 'qc' )
3611	DO k = nzb, nzt
3612	sums_wsqcs_ws_l(k,tn) = sums_wsqcs_ws_l(k,tn) &
3613	+ ( flux_t(k) &
3614	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3615	* ( w(k,j,i) - hom(k,1,3,0) ) &
3616	+ diss_t(k) &
3617	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3618	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3619	) * weight_substep(intermediate_timestep_count)
3620	ENDDO
3621	CASE ( 'qr' )
3622	DO k = nzb, nzt
3623	sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) &
3624	+ ( flux_t(k) &
3625	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3626	* ( w(k,j,i) - hom(k,1,3,0) ) &
3627	+ diss_t(k) &
3628	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3629	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3630	) * weight_substep(intermediate_timestep_count)
3631	ENDDO
3632	CASE ( 'nc' )
3633	DO k = nzb, nzt
3634	sums_wsncs_ws_l(k,tn) = sums_wsncs_ws_l(k,tn) &
3635	+ ( flux_t(k) &
3636	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3637	* ( w(k,j,i) - hom(k,1,3,0) ) &
3638	+ diss_t(k) &
3639	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3640	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3641	) * weight_substep(intermediate_timestep_count)
3642	ENDDO
3643	CASE ( 'nr' )
3644	DO k = nzb, nzt
3645	sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) &
3646	+ ( flux_t(k) &
3647	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3648	* ( w(k,j,i) - hom(k,1,3,0) ) &
3649	+ diss_t(k) &
3650	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3651	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3652	) * weight_substep(intermediate_timestep_count)
3653	ENDDO
3654	CASE ( 's' )
3655	DO k = nzb, nzt
3656	sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn) &
3657	+ ( flux_t(k) &
3658	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3659	* ( w(k,j,i) - hom(k,1,3,0) ) &
3660	+ diss_t(k) &
3661	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3662	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3663	) * weight_substep(intermediate_timestep_count)
3664	ENDDO
3665
3666
3667	END SELECT
3668
3669	ENDDO
3670	ENDDO
3671
3672	END SUBROUTINE advec_s_ws
3673
3674
3675	!------------------------------------------------------------------------------!
3676	! Description:
3677	! ------------
3678	!> Advection of u - Call for all grid points
3679	!------------------------------------------------------------------------------!
3680	SUBROUTINE advec_u_ws
3681
3682	USE arrays_3d, &
3683	ONLY: ddzw, drho_air, tend, u, v, w, rho_air_zw
3684
3685	USE constants, &
3686	ONLY: adv_mom_1, adv_mom_3, adv_mom_5
3687
3688	USE control_parameters, &
3689	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
3690
3691	USE grid_variables, &
3692	ONLY: ddx, ddy
3693
3694	USE indices, &
3695	ONLY: nxl, nxlu, nxr, nyn, nys, nzb, nzb_max, nzt, advc_flags_1
3696
3697	USE kinds
3698
3699	USE statistics, &
3700	ONLY: hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep
3701
3702	IMPLICIT NONE
3703
3704	INTEGER(iwp) :: i !<
3705	INTEGER(iwp) :: ibit9 !<
3706	INTEGER(iwp) :: ibit10 !<
3707	INTEGER(iwp) :: ibit11 !<
3708	INTEGER(iwp) :: ibit12 !<
3709	INTEGER(iwp) :: ibit13 !<
3710	INTEGER(iwp) :: ibit14 !<
3711	INTEGER(iwp) :: ibit15 !<
3712	INTEGER(iwp) :: ibit16 !<
3713	INTEGER(iwp) :: ibit17 !<
3714	INTEGER(iwp) :: j !<
3715	INTEGER(iwp) :: k !<
3716	INTEGER(iwp) :: k_mm !<
3717	INTEGER(iwp) :: k_pp !<
3718	INTEGER(iwp) :: k_ppp !<
3719	INTEGER(iwp) :: tn = 0 !<
3720
3721	REAL(wp) :: diss_d !<
3722	REAL(wp) :: div !<
3723	REAL(wp) :: flux_d !<
3724	REAL(wp) :: gu !<
3725	REAL(wp) :: gv !<
3726	REAL(wp) :: v_comp !<
3727	REAL(wp) :: w_comp !<
3728
3729	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local_u !<
3730	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local_u !<
3731
3732	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_u !<
3733	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_u !<
3734
3735	REAL(wp), DIMENSION(nzb:nzt) :: diss_n !<
3736	REAL(wp), DIMENSION(nzb:nzt) :: diss_r !<
3737	REAL(wp), DIMENSION(nzb:nzt) :: diss_t !<
3738	REAL(wp), DIMENSION(nzb:nzt) :: flux_n !<
3739	REAL(wp), DIMENSION(nzb:nzt) :: flux_r !<
3740	REAL(wp), DIMENSION(nzb:nzt) :: flux_t !<
3741	REAL(wp), DIMENSION(nzb:nzt) :: u_comp !<
3742
3743	gu = 2.0_wp * u_gtrans
3744	gv = 2.0_wp * v_gtrans
3745
3746	!
3747	!-- Compute the fluxes for the whole left boundary of the processor domain.
3748	i = nxlu
3749	DO j = nys, nyn
3750	DO k = nzb+1, nzb_max
3751
3752	ibit11 = IBITS(advc_flags_1(k,j,i-1),11,1)
3753	ibit10 = IBITS(advc_flags_1(k,j,i-1),10,1)
3754	ibit9 = IBITS(advc_flags_1(k,j,i-1),9,1)
3755
3756	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
3757	swap_flux_x_local_u(k,j) = u_comp(k) * ( &
3758	( 37.0_wp * ibit11 * adv_mom_5 &
3759	+ 7.0_wp * ibit10 * adv_mom_3 &
3760	+ ibit9 * adv_mom_1 &
3761	) * &
3762	( u(k,j,i) + u(k,j,i-1) ) &
3763	- ( 8.0_wp * ibit11 * adv_mom_5 &
3764	+ ibit10 * adv_mom_3 &
3765	) * &
3766	( u(k,j,i+1) + u(k,j,i-2) ) &
3767	+ ( ibit11 * adv_mom_5 &
3768	) * &
3769	( u(k,j,i+2) + u(k,j,i-3) ) &
3770	)
3771
3772	swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * ( &
3773	( 10.0_wp * ibit11 * adv_mom_5 &
3774	+ 3.0_wp * ibit10 * adv_mom_3 &
3775	+ ibit9 * adv_mom_1 &
3776	) * &
3777	( u(k,j,i) - u(k,j,i-1) ) &
3778	- ( 5.0_wp * ibit11 * adv_mom_5 &
3779	+ ibit10 * adv_mom_3 &
3780	) * &
3781	( u(k,j,i+1) - u(k,j,i-2) ) &
3782	+ ( ibit11 * adv_mom_5 &
3783	) * &
3784	( u(k,j,i+2) - u(k,j,i-3) ) &
3785	)
3786
3787	ENDDO
3788
3789	DO k = nzb_max+1, nzt
3790
3791	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
3792	swap_flux_x_local_u(k,j) = u_comp(k) * ( &
3793	37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) &
3794	- 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) &
3795	+ ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
3796	swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * ( &
3797	10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) &
3798	- 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) &
3799	+ ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
3800
3801	ENDDO
3802	ENDDO
3803
3804	DO i = nxlu, nxr
3805	!
3806	!-- The following loop computes the fluxes for the south boundary points
3807	j = nys
3808	DO k = nzb+1, nzb_max
3809
3810	ibit14 = IBITS(advc_flags_1(k,j-1,i),14,1)
3811	ibit13 = IBITS(advc_flags_1(k,j-1,i),13,1)
3812	ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)
3813
3814	v_comp = v(k,j,i) + v(k,j,i-1) - gv
3815	swap_flux_y_local_u(k) = v_comp * ( &
3816	( 37.0_wp * ibit14 * adv_mom_5 &
3817	+ 7.0_wp * ibit13 * adv_mom_3 &
3818	+ ibit12 * adv_mom_1 &
3819	) * &
3820	( u(k,j,i) + u(k,j-1,i) ) &
3821	- ( 8.0_wp * ibit14 * adv_mom_5 &
3822	+ ibit13 * adv_mom_3 &
3823	) * &
3824	( u(k,j+1,i) + u(k,j-2,i) ) &
3825	+ ( ibit14 * adv_mom_5 &
3826	) * &
3827	( u(k,j+2,i) + u(k,j-3,i) ) &
3828	)
3829
3830	swap_diss_y_local_u(k) = - ABS ( v_comp ) * ( &
3831	( 10.0_wp * ibit14 * adv_mom_5 &
3832	+ 3.0_wp * ibit13 * adv_mom_3 &
3833	+ ibit12 * adv_mom_1 &
3834	) * &
3835	( u(k,j,i) - u(k,j-1,i) ) &
3836	- ( 5.0_wp * ibit14 * adv_mom_5 &
3837	+ ibit13 * adv_mom_3 &
3838	) * &
3839	( u(k,j+1,i) - u(k,j-2,i) ) &
3840	+ ( ibit14 * adv_mom_5 &
3841	) * &
3842	( u(k,j+2,i) - u(k,j-3,i) ) &
3843	)
3844
3845	ENDDO
3846
3847	DO k = nzb_max+1, nzt
3848
3849	v_comp = v(k,j,i) + v(k,j,i-1) - gv
3850	swap_flux_y_local_u(k) = v_comp * ( &
3851	37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) &
3852	- 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) &
3853	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
3854	swap_diss_y_local_u(k) = - ABS(v_comp) * ( &
3855	10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) &
3856	- 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) &
3857	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
3858
3859	ENDDO
3860	!
3861	!-- Computation of interior fluxes and tendency terms
3862	DO j = nys, nyn
3863
3864	flux_t(0) = 0.0_wp
3865	diss_t(0) = 0.0_wp
3866	flux_d = 0.0_wp
3867	diss_d = 0.0_wp
3868
3869	DO k = nzb+1, nzb_max
3870
3871	ibit11 = IBITS(advc_flags_1(k,j,i),11,1)
3872	ibit10 = IBITS(advc_flags_1(k,j,i),10,1)
3873	ibit9 = IBITS(advc_flags_1(k,j,i),9,1)
3874
3875	u_comp(k) = u(k,j,i+1) + u(k,j,i)
3876	flux_r(k) = ( u_comp(k) - gu ) * ( &
3877	( 37.0_wp * ibit11 * adv_mom_5 &
3878	+ 7.0_wp * ibit10 * adv_mom_3 &
3879	+ ibit9 * adv_mom_1 &
3880	) * &
3881	( u(k,j,i+1) + u(k,j,i) ) &
3882	- ( 8.0_wp * ibit11 * adv_mom_5 &
3883	+ ibit10 * adv_mom_3 &
3884	) * &
3885	( u(k,j,i+2) + u(k,j,i-1) ) &
3886	+ ( ibit11 * adv_mom_5 &
3887	) * &
3888	( u(k,j,i+3) + u(k,j,i-2) ) &
3889	)
3890
3891	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
3892	( 10.0_wp * ibit11 * adv_mom_5 &
3893	+ 3.0_wp * ibit10 * adv_mom_3 &
3894	+ ibit9 * adv_mom_1 &
3895	) * &
3896	( u(k,j,i+1) - u(k,j,i) ) &
3897	- ( 5.0_wp * ibit11 * adv_mom_5 &
3898	+ ibit10 * adv_mom_3 &
3899	) * &
3900	( u(k,j,i+2) - u(k,j,i-1) ) &
3901	+ ( ibit11 * adv_mom_5 &
3902	) * &
3903	( u(k,j,i+3) - u(k,j,i-2) ) &
3904	)
3905
3906	ibit14 = IBITS(advc_flags_1(k,j,i),14,1)
3907	ibit13 = IBITS(advc_flags_1(k,j,i),13,1)
3908	ibit12 = IBITS(advc_flags_1(k,j,i),12,1)
3909
3910	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
3911	flux_n(k) = v_comp * ( &
3912	( 37.0_wp * ibit14 * adv_mom_5 &
3913	+ 7.0_wp * ibit13 * adv_mom_3 &
3914	+ ibit12 * adv_mom_1 &
3915	) * &
3916	( u(k,j+1,i) + u(k,j,i) ) &
3917	- ( 8.0_wp * ibit14 * adv_mom_5 &
3918	+ ibit13 * adv_mom_3 &
3919	) * &
3920	( u(k,j+2,i) + u(k,j-1,i) ) &
3921	+ ( ibit14 * adv_mom_5 &
3922	) * &
3923	( u(k,j+3,i) + u(k,j-2,i) ) &
3924	)
3925
3926	diss_n(k) = - ABS ( v_comp ) * ( &
3927	( 10.0_wp * ibit14 * adv_mom_5 &
3928	+ 3.0_wp * ibit13 * adv_mom_3 &
3929	+ ibit12 * adv_mom_1 &
3930	) * &
3931	( u(k,j+1,i) - u(k,j,i) ) &
3932	- ( 5.0_wp * ibit14 * adv_mom_5 &
3933	+ ibit13 * adv_mom_3 &
3934	) * &
3935	( u(k,j+2,i) - u(k,j-1,i) ) &
3936	+ ( ibit14 * adv_mom_5 &
3937	) * &
3938	( u(k,j+3,i) - u(k,j-2,i) ) &
3939	)
3940	!
3941	!-- k index has to be modified near bottom and top, else array
3942	!-- subscripts will be exceeded.
3943	ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
3944	ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
3945	ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
3946
3947	k_ppp = k + 3 * ibit17
3948	k_pp = k + 2 * ( 1 - ibit15 )
3949	k_mm = k - 2 * ibit17
3950
3951	w_comp = w(k,j,i) + w(k,j,i-1)
3952	flux_t(k) = w_comp * rho_air_zw(k) * ( &
3953	( 37.0_wp * ibit17 * adv_mom_5 &
3954	+ 7.0_wp * ibit16 * adv_mom_3 &
3955	+ ibit15 * adv_mom_1 &
3956	) * &
3957	( u(k+1,j,i) + u(k,j,i) ) &
3958	- ( 8.0_wp * ibit17 * adv_mom_5 &
3959	+ ibit16 * adv_mom_3 &
3960	) * &
3961	( u(k_pp,j,i) + u(k-1,j,i) ) &
3962	+ ( ibit17 * adv_mom_5 &
3963	) * &
3964	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
3965	)
3966
3967	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
3968	( 10.0_wp * ibit17 * adv_mom_5 &
3969	+ 3.0_wp * ibit16 * adv_mom_3 &
3970	+ ibit15 * adv_mom_1 &
3971	) * &
3972	( u(k+1,j,i) - u(k,j,i) ) &
3973	- ( 5.0_wp * ibit17 * adv_mom_5 &
3974	+ ibit16 * adv_mom_3 &
3975	) * &
3976	( u(k_pp,j,i) - u(k-1,j,i) ) &
3977	+ ( ibit17 * adv_mom_5 &
3978	) * &
3979	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
3980	)
3981	!
3982	!-- Calculate the divergence of the velocity field. A respective
3983	!-- correction is needed to overcome numerical instabilities caused
3984	!-- by a not sufficient reduction of divergences near topography.
3985	div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 ) &
3986	- ( u(k,j,i) + u(k,j,i-1) ) &
3987	* ( IBITS(advc_flags_1(k,j,i-1),9,1) &
3988	+ IBITS(advc_flags_1(k,j,i-1),10,1) &
3989	+ IBITS(advc_flags_1(k,j,i-1),11,1) &
3990	) &
3991	) * ddx &
3992	+ ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 ) &
3993	- ( v(k,j,i) + v(k,j,i-1 ) ) &
3994	* ( IBITS(advc_flags_1(k,j-1,i),12,1) &
3995	+ IBITS(advc_flags_1(k,j-1,i),13,1) &
3996	+ IBITS(advc_flags_1(k,j-1,i),14,1) &
3997	) &
3998	) * ddy &
3999	+ ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 ) &
4000	- ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) &
4001	* ( IBITS(advc_flags_1(k-1,j,i),15,1) &
4002	+ IBITS(advc_flags_1(k-1,j,i),16,1) &
4003	+ IBITS(advc_flags_1(k-1,j,i),17,1) &
4004	) &
4005	) * drho_air(k) * ddzw(k) &
4006	) * 0.5_wp
4007
4008
4009
4010	tend(k,j,i) = tend(k,j,i) - ( &
4011	( flux_r(k) + diss_r(k) &
4012	- swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
4013	+ ( flux_n(k) + diss_n(k) &
4014	- swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy &
4015	+ ( ( flux_t(k) + diss_t(k) ) &
4016	- ( flux_d + diss_d ) &
4017	) * drho_air(k) * ddzw(k) &
4018	) + div * u(k,j,i)
4019
4020	swap_flux_x_local_u(k,j) = flux_r(k)
4021	swap_diss_x_local_u(k,j) = diss_r(k)
4022	swap_flux_y_local_u(k) = flux_n(k)
4023	swap_diss_y_local_u(k) = diss_n(k)
4024	flux_d = flux_t(k)
4025	diss_d = diss_t(k)
4026	!
4027	!-- Statistical Evaluation of u'u'. The factor has to be applied
4028	!-- for right evaluation when gallilei_trans = .T. .
4029	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
4030	+ ( flux_r(k) &
4031	* ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
4032	/ ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &
4033	+ diss_r(k) &
4034	* ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
4035	/ ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &
4036	) * weight_substep(intermediate_timestep_count)
4037	!
4038	!-- Statistical Evaluation of w'u'.
4039	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
4040	+ ( flux_t(k) &
4041	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4042	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
4043	+ diss_t(k) &
4044	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4045	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
4046	) * weight_substep(intermediate_timestep_count)
4047
4048	ENDDO
4049
4050	DO k = nzb_max+1, nzt
4051
4052	u_comp(k) = u(k,j,i+1) + u(k,j,i)
4053	flux_r(k) = ( u_comp(k) - gu ) * ( &
4054	37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) &
4055	- 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) &
4056	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
4057	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
4058	10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) &
4059	- 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) &
4060	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
4061
4062	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
4063	flux_n(k) = v_comp * ( &
4064	37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) &
4065	- 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) &
4066	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
4067	diss_n(k) = - ABS( v_comp ) * ( &
4068	10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) &
4069	- 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) &
4070	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
4071	!
4072	!-- k index has to be modified near bottom and top, else array
4073	!-- subscripts will be exceeded.
4074	ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
4075	ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
4076	ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
4077
4078	k_ppp = k + 3 * ibit17
4079	k_pp = k + 2 * ( 1 - ibit15 )
4080	k_mm = k - 2 * ibit17
4081
4082	w_comp = w(k,j,i) + w(k,j,i-1)
4083	flux_t(k) = w_comp * rho_air_zw(k) * ( &
4084	( 37.0_wp * ibit17 * adv_mom_5 &
4085	+ 7.0_wp * ibit16 * adv_mom_3 &
4086	+ ibit15 * adv_mom_1 &
4087	) * &
4088	( u(k+1,j,i) + u(k,j,i) ) &
4089	- ( 8.0_wp * ibit17 * adv_mom_5 &
4090	+ ibit16 * adv_mom_3 &
4091	) * &
4092	( u(k_pp,j,i) + u(k-1,j,i) ) &
4093	+ ( ibit17 * adv_mom_5 &
4094	) * &
4095	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
4096	)
4097
4098	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
4099	( 10.0_wp * ibit17 * adv_mom_5 &
4100	+ 3.0_wp * ibit16 * adv_mom_3 &
4101	+ ibit15 * adv_mom_1 &
4102	) * &
4103	( u(k+1,j,i) - u(k,j,i) ) &
4104	- ( 5.0_wp * ibit17 * adv_mom_5 &
4105	+ ibit16 * adv_mom_3 &
4106	) * &
4107	( u(k_pp,j,i) - u(k-1,j,i) ) &
4108	+ ( ibit17 * adv_mom_5 &
4109	) * &
4110	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
4111	)
4112	!
4113	!-- Calculate the divergence of the velocity field. A respective
4114	!-- correction is needed to overcome numerical instabilities caused
4115	!-- by a not sufficient reduction of divergences near topography.
4116	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
4117	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
4118	+ ( w_comp * rho_air_zw(k) - &
4119	( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) &
4120	) * drho_air(k) * ddzw(k) &
4121	) * 0.5_wp
4122
4123	tend(k,j,i) = tend(k,j,i) - ( &
4124	( flux_r(k) + diss_r(k) &
4125	- swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
4126	+ ( flux_n(k) + diss_n(k) &
4127	- swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy &
4128	+ ( ( flux_t(k) + diss_t(k) ) &
4129	- ( flux_d + diss_d ) &
4130	) * drho_air(k) * ddzw(k) &
4131	) + div * u(k,j,i)
4132
4133	swap_flux_x_local_u(k,j) = flux_r(k)
4134	swap_diss_x_local_u(k,j) = diss_r(k)
4135	swap_flux_y_local_u(k) = flux_n(k)
4136	swap_diss_y_local_u(k) = diss_n(k)
4137	flux_d = flux_t(k)
4138	diss_d = diss_t(k)
4139	!
4140	!-- Statistical Evaluation of u'u'. The factor has to be applied
4141	!-- for right evaluation when gallilei_trans = .T. .
4142	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
4143	+ ( flux_r(k) &
4144	* ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
4145	/ ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &
4146	+ diss_r(k) &
4147	* ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
4148	/ ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &
4149	) * weight_substep(intermediate_timestep_count)
4150	!
4151	!-- Statistical Evaluation of w'u'.
4152	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
4153	+ ( flux_t(k) &
4154	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4155	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
4156	+ diss_t(k) &
4157	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4158	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
4159	) * weight_substep(intermediate_timestep_count)
4160	ENDDO
4161	ENDDO
4162	ENDDO
4163	sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
4164
4165
4166	END SUBROUTINE advec_u_ws
4167
4168
4169	!------------------------------------------------------------------------------!
4170	! Description:
4171	! ------------
4172	!> Advection of v - Call for all grid points
4173	!------------------------------------------------------------------------------!
4174	SUBROUTINE advec_v_ws
4175
4176	USE arrays_3d, &
4177	ONLY: ddzw, drho_air, tend, u, v, w, rho_air_zw
4178
4179	USE constants, &
4180	ONLY: adv_mom_1, adv_mom_3, adv_mom_5
4181
4182	USE control_parameters, &
4183	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
4184
4185	USE grid_variables, &
4186	ONLY: ddx, ddy
4187
4188	USE indices, &
4189	ONLY: nxl, nxr, nyn, nys, nysv, nzb, nzb_max, nzt, advc_flags_1
4190
4191	USE kinds
4192
4193	USE statistics, &
4194	ONLY: hom, sums_vs2_ws_l, sums_wsvs_ws_l, weight_substep
4195
4196	IMPLICIT NONE
4197
4198
4199	INTEGER(iwp) :: i !<
4200	INTEGER(iwp) :: ibit18 !<
4201	INTEGER(iwp) :: ibit19 !<
4202	INTEGER(iwp) :: ibit20 !<
4203	INTEGER(iwp) :: ibit21 !<
4204	INTEGER(iwp) :: ibit22 !<
4205	INTEGER(iwp) :: ibit23 !<
4206	INTEGER(iwp) :: ibit24 !<
4207	INTEGER(iwp) :: ibit25 !<
4208	INTEGER(iwp) :: ibit26 !<
4209	INTEGER(iwp) :: j !<
4210	INTEGER(iwp) :: k !<
4211	INTEGER(iwp) :: k_mm !<
4212	INTEGER(iwp) :: k_pp !<
4213	INTEGER(iwp) :: k_ppp !<
4214	INTEGER(iwp) :: tn = 0 !<
4215
4216	REAL(wp) :: diss_d !<
4217	REAL(wp) :: div !<
4218	REAL(wp) :: flux_d !<
4219	REAL(wp) :: gu !<
4220	REAL(wp) :: gv !<
4221	REAL(wp) :: u_comp !<
4222	REAL(wp) :: w_comp !<
4223
4224	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local_v !<
4225	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local_v !<
4226
4227	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_v !<
4228	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_v !<
4229
4230	REAL(wp), DIMENSION(nzb:nzt) :: diss_n !<
4231	REAL(wp), DIMENSION(nzb:nzt) :: diss_r !<
4232	REAL(wp), DIMENSION(nzb:nzt) :: diss_t !<
4233	REAL(wp), DIMENSION(nzb:nzt) :: flux_n !<
4234	REAL(wp), DIMENSION(nzb:nzt) :: flux_r !<
4235	REAL(wp), DIMENSION(nzb:nzt) :: flux_t !<
4236	REAL(wp), DIMENSION(nzb:nzt) :: v_comp !<
4237
4238	gu = 2.0_wp * u_gtrans
4239	gv = 2.0_wp * v_gtrans
4240	!
4241	!-- First compute the whole left boundary of the processor domain
4242	i = nxl
4243	DO j = nysv, nyn
4244	DO k = nzb+1, nzb_max
4245
4246	ibit20 = IBITS(advc_flags_1(k,j,i-1),20,1)
4247	ibit19 = IBITS(advc_flags_1(k,j,i-1),19,1)
4248	ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)
4249
4250	u_comp = u(k,j-1,i) + u(k,j,i) - gu
4251	swap_flux_x_local_v(k,j) = u_comp * ( &
4252	( 37.0_wp * ibit20 * adv_mom_5 &
4253	+ 7.0_wp * ibit19 * adv_mom_3 &
4254	+ ibit18 * adv_mom_1 &
4255	) * &
4256	( v(k,j,i) + v(k,j,i-1) ) &
4257	- ( 8.0_wp * ibit20 * adv_mom_5 &
4258	+ ibit19 * adv_mom_3 &
4259	) * &
4260	( v(k,j,i+1) + v(k,j,i-2) ) &
4261	+ ( ibit20 * adv_mom_5 &
4262	) * &
4263	( v(k,j,i+2) + v(k,j,i-3) ) &
4264	)
4265
4266	swap_diss_x_local_v(k,j) = - ABS( u_comp ) * ( &
4267	( 10.0_wp * ibit20 * adv_mom_5 &
4268	+ 3.0_wp * ibit19 * adv_mom_3 &
4269	+ ibit18 * adv_mom_1 &
4270	) * &
4271	( v(k,j,i) - v(k,j,i-1) ) &
4272	- ( 5.0_wp * ibit20 * adv_mom_5 &
4273	+ ibit19 * adv_mom_3 &
4274	) * &
4275	( v(k,j,i+1) - v(k,j,i-2) ) &
4276	+ ( ibit20 * adv_mom_5 &
4277	) * &
4278	( v(k,j,i+2) - v(k,j,i-3) ) &
4279	)
4280
4281	ENDDO
4282
4283	DO k = nzb_max+1, nzt
4284
4285	u_comp = u(k,j-1,i) + u(k,j,i) - gu
4286	swap_flux_x_local_v(k,j) = u_comp * ( &
4287	37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) &
4288	- 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) &
4289	+ ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
4290	swap_diss_x_local_v(k,j) = - ABS( u_comp ) * ( &
4291	10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) &
4292	- 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) &
4293	+ ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
4294
4295	ENDDO
4296
4297	ENDDO
4298
4299	DO i = nxl, nxr
4300
4301	j = nysv
4302	DO k = nzb+1, nzb_max
4303
4304	ibit23 = IBITS(advc_flags_1(k,j-1,i),23,1)
4305	ibit22 = IBITS(advc_flags_1(k,j-1,i),22,1)
4306	ibit21 = IBITS(advc_flags_1(k,j-1,i),21,1)
4307
4308	v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv
4309	swap_flux_y_local_v(k) = v_comp(k) * ( &
4310	( 37.0_wp * ibit23 * adv_mom_5 &
4311	+ 7.0_wp * ibit22 * adv_mom_3 &
4312	+ ibit21 * adv_mom_1 &
4313	) * &
4314	( v(k,j,i) + v(k,j-1,i) ) &
4315	- ( 8.0_wp * ibit23 * adv_mom_5 &
4316	+ ibit22 * adv_mom_3 &
4317	) * &
4318	( v(k,j+1,i) + v(k,j-2,i) ) &
4319	+ ( ibit23 * adv_mom_5 &
4320	) * &
4321	( v(k,j+2,i) + v(k,j-3,i) ) &
4322	)
4323
4324	swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * ( &
4325	( 10.0_wp * ibit23 * adv_mom_5 &
4326	+ 3.0_wp * ibit22 * adv_mom_3 &
4327	+ ibit21 * adv_mom_1 &
4328	) * &
4329	( v(k,j,i) - v(k,j-1,i) ) &
4330	- ( 5.0_wp * ibit23 * adv_mom_5 &
4331	+ ibit22 * adv_mom_3 &
4332	) * &
4333	( v(k,j+1,i) - v(k,j-2,i) ) &
4334	+ ( ibit23 * adv_mom_5 &
4335	) * &
4336	( v(k,j+2,i) - v(k,j-3,i) ) &
4337	)
4338
4339	ENDDO
4340
4341	DO k = nzb_max+1, nzt
4342
4343	v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv
4344	swap_flux_y_local_v(k) = v_comp(k) * ( &
4345	37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) &
4346	- 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) &
4347	+ ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
4348	swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * ( &
4349	10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) &
4350	- 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) &
4351	+ ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
4352
4353	ENDDO
4354
4355	DO j = nysv, nyn
4356
4357	flux_t(0) = 0.0_wp
4358	diss_t(0) = 0.0_wp
4359	flux_d = 0.0_wp
4360	diss_d = 0.0_wp
4361
4362	DO k = nzb+1, nzb_max
4363
4364	ibit20 = IBITS(advc_flags_1(k,j,i),20,1)
4365	ibit19 = IBITS(advc_flags_1(k,j,i),19,1)
4366	ibit18 = IBITS(advc_flags_1(k,j,i),18,1)
4367
4368	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
4369	flux_r(k) = u_comp * ( &
4370	( 37.0_wp * ibit20 * adv_mom_5 &
4371	+ 7.0_wp * ibit19 * adv_mom_3 &
4372	+ ibit18 * adv_mom_1 &
4373	) * &
4374	( v(k,j,i+1) + v(k,j,i) ) &
4375	- ( 8.0_wp * ibit20 * adv_mom_5 &
4376	+ ibit19 * adv_mom_3 &
4377	) * &
4378	( v(k,j,i+2) + v(k,j,i-1) ) &
4379	+ ( ibit20 * adv_mom_5 &
4380	) * &
4381	( v(k,j,i+3) + v(k,j,i-2) ) &
4382	)
4383
4384	diss_r(k) = - ABS( u_comp ) * ( &
4385	( 10.0_wp * ibit20 * adv_mom_5 &
4386	+ 3.0_wp * ibit19 * adv_mom_3 &
4387	+ ibit18 * adv_mom_1 &
4388	) * &
4389	( v(k,j,i+1) - v(k,j,i) ) &
4390	- ( 5.0_wp * ibit20 * adv_mom_5 &
4391	+ ibit19 * adv_mom_3 &
4392	) * &
4393	( v(k,j,i+2) - v(k,j,i-1) ) &
4394	+ ( ibit20 * adv_mom_5 &
4395	) * &
4396	( v(k,j,i+3) - v(k,j,i-2) ) &
4397	)
4398
4399	ibit23 = IBITS(advc_flags_1(k,j,i),23,1)
4400	ibit22 = IBITS(advc_flags_1(k,j,i),22,1)
4401	ibit21 = IBITS(advc_flags_1(k,j,i),21,1)
4402
4403	v_comp(k) = v(k,j+1,i) + v(k,j,i)
4404	flux_n(k) = ( v_comp(k) - gv ) * ( &
4405	( 37.0_wp * ibit23 * adv_mom_5 &
4406	+ 7.0_wp * ibit22 * adv_mom_3 &
4407	+ ibit21 * adv_mom_1 &
4408	) * &
4409	( v(k,j+1,i) + v(k,j,i) ) &
4410	- ( 8.0_wp * ibit23 * adv_mom_5 &
4411	+ ibit22 * adv_mom_3 &
4412	) * &
4413	( v(k,j+2,i) + v(k,j-1,i) ) &
4414	+ ( ibit23 * adv_mom_5 &
4415	) * &
4416	( v(k,j+3,i) + v(k,j-2,i) ) &
4417	)
4418
4419	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
4420	( 10.0_wp * ibit23 * adv_mom_5 &
4421	+ 3.0_wp * ibit22 * adv_mom_3 &
4422	+ ibit21 * adv_mom_1 &
4423	) * &
4424	( v(k,j+1,i) - v(k,j,i) ) &
4425	- ( 5.0_wp * ibit23 * adv_mom_5 &
4426	+ ibit22 * adv_mom_3 &
4427	) * &
4428	( v(k,j+2,i) - v(k,j-1,i) ) &
4429	+ ( ibit23 * adv_mom_5 &
4430	) * &
4431	( v(k,j+3,i) - v(k,j-2,i) ) &
4432	)
4433	!
4434	!-- k index has to be modified near bottom and top, else array
4435	!-- subscripts will be exceeded.
4436	ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
4437	ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
4438	ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
4439
4440	k_ppp = k + 3 * ibit26
4441	k_pp = k + 2 * ( 1 - ibit24 )
4442	k_mm = k - 2 * ibit26
4443
4444	w_comp = w(k,j-1,i) + w(k,j,i)
4445	flux_t(k) = w_comp * rho_air_zw(k) * ( &
4446	( 37.0_wp * ibit26 * adv_mom_5 &
4447	+ 7.0_wp * ibit25 * adv_mom_3 &
4448	+ ibit24 * adv_mom_1 &
4449	) * &
4450	( v(k+1,j,i) + v(k,j,i) ) &
4451	- ( 8.0_wp * ibit26 * adv_mom_5 &
4452	+ ibit25 * adv_mom_3 &
4453	) * &
4454	( v(k_pp,j,i) + v(k-1,j,i) ) &
4455	+ ( ibit26 * adv_mom_5 &
4456	) * &
4457	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
4458	)
4459
4460	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
4461	( 10.0_wp * ibit26 * adv_mom_5 &
4462	+ 3.0_wp * ibit25 * adv_mom_3 &
4463	+ ibit24 * adv_mom_1 &
4464	) * &
4465	( v(k+1,j,i) - v(k,j,i) ) &
4466	- ( 5.0_wp * ibit26 * adv_mom_5 &
4467	+ ibit25 * adv_mom_3 &
4468	) * &
4469	( v(k_pp,j,i) - v(k-1,j,i) ) &
4470	+ ( ibit26 * adv_mom_5 &
4471	) * &
4472	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
4473	)
4474	!
4475	!-- Calculate the divergence of the velocity field. A respective
4476	!-- correction is needed to overcome numerical instabilities caused
4477	!-- by a not sufficient reduction of divergences near topography.
4478	div = ( ( ( u_comp + gu ) &
4479	* ( ibit18 + ibit19 + ibit20 ) &
4480	- ( u(k,j-1,i) + u(k,j,i) ) &
4481	* ( IBITS(advc_flags_1(k,j,i-1),18,1) &
4482	+ IBITS(advc_flags_1(k,j,i-1),19,1) &
4483	+ IBITS(advc_flags_1(k,j,i-1),20,1) &
4484	) &
4485	) * ddx &
4486	+ ( v_comp(k) &
4487	* ( ibit21 + ibit22 + ibit23 ) &
4488	- ( v(k,j,i) + v(k,j-1,i) ) &
4489	* ( IBITS(advc_flags_1(k,j-1,i),21,1) &
4490	+ IBITS(advc_flags_1(k,j-1,i),22,1) &
4491	+ IBITS(advc_flags_1(k,j-1,i),23,1) &
4492	) &
4493	) * ddy &
4494	+ ( w_comp * rho_air_zw(k) &
4495	* ( ibit24 + ibit25 + ibit26 ) &
4496	- ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) &
4497	* ( IBITS(advc_flags_1(k-1,j,i),24,1) &
4498	+ IBITS(advc_flags_1(k-1,j,i),25,1) &
4499	+ IBITS(advc_flags_1(k-1,j,i),26,1) &
4500	) &
4501	) * drho_air(k) * ddzw(k) &
4502	) * 0.5_wp
4503
4504
4505	tend(k,j,i) = tend(k,j,i) - ( &
4506	( flux_r(k) + diss_r(k) &
4507	- swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j) &
4508	) * ddx &
4509	+ ( flux_n(k) + diss_n(k) &
4510	- swap_flux_y_local_v(k) - swap_diss_y_local_v(k) &
4511	) * ddy &
4512	+ ( ( flux_t(k) + diss_t(k) ) &
4513	- ( flux_d + diss_d ) &
4514	) * drho_air(k) * ddzw(k) &
4515	) + v(k,j,i) * div
4516
4517	swap_flux_x_local_v(k,j) = flux_r(k)
4518	swap_diss_x_local_v(k,j) = diss_r(k)
4519	swap_flux_y_local_v(k) = flux_n(k)
4520	swap_diss_y_local_v(k) = diss_n(k)
4521	flux_d = flux_t(k)
4522	diss_d = diss_t(k)
4523
4524	!
4525	!-- Statistical Evaluation of v'v'. The factor has to be applied
4526	!-- for right evaluation when gallilei_trans = .T. .
4527	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
4528	+ ( flux_n(k) &
4529	* ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
4530	/ ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &
4531	+ diss_n(k) &
4532	* ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
4533	/ ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &
4534	) * weight_substep(intermediate_timestep_count)
4535	!
4536	!-- Statistical Evaluation of w'u'.
4537	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
4538	+ ( flux_t(k) &
4539	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4540	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
4541	+ diss_t(k) &
4542	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4543	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
4544	) * weight_substep(intermediate_timestep_count)
4545
4546	ENDDO
4547
4548	DO k = nzb_max+1, nzt
4549
4550	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
4551	flux_r(k) = u_comp * ( &
4552	37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) &
4553	- 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) &
4554	+ ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
4555
4556	diss_r(k) = - ABS( u_comp ) * ( &
4557	10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) &
4558	- 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) &
4559	+ ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
4560
4561
4562	v_comp(k) = v(k,j+1,i) + v(k,j,i)
4563	flux_n(k) = ( v_comp(k) - gv ) * ( &
4564	37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) &
4565	- 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) &
4566	+ ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
4567
4568	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
4569	10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) &
4570	- 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) &
4571	+ ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
4572	!
4573	!-- k index has to be modified near bottom and top, else array
4574	!-- subscripts will be exceeded.
4575	ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
4576	ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
4577	ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
4578
4579	k_ppp = k + 3 * ibit26
4580	k_pp = k + 2 * ( 1 - ibit24 )
4581	k_mm = k - 2 * ibit26
4582
4583	w_comp = w(k,j-1,i) + w(k,j,i)
4584	flux_t(k) = w_comp * rho_air_zw(k) * ( &
4585	( 37.0_wp * ibit26 * adv_mom_5 &
4586	+ 7.0_wp * ibit25 * adv_mom_3 &
4587	+ ibit24 * adv_mom_1 &
4588	) * &
4589	( v(k+1,j,i) + v(k,j,i) ) &
4590	- ( 8.0_wp * ibit26 * adv_mom_5 &
4591	+ ibit25 * adv_mom_3 &
4592	) * &
4593	( v(k_pp,j,i) + v(k-1,j,i) ) &
4594	+ ( ibit26 * adv_mom_5 &
4595	) * &
4596	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
4597	)
4598
4599	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
4600	( 10.0_wp * ibit26 * adv_mom_5 &
4601	+ 3.0_wp * ibit25 * adv_mom_3 &
4602	+ ibit24 * adv_mom_1 &
4603	) * &
4604	( v(k+1,j,i) - v(k,j,i) ) &
4605	- ( 5.0_wp * ibit26 * adv_mom_5 &
4606	+ ibit25 * adv_mom_3 &
4607	) * &
4608	( v(k_pp,j,i) - v(k-1,j,i) ) &
4609	+ ( ibit26 * adv_mom_5 &
4610	) * &
4611	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
4612	)
4613	!
4614	!-- Calculate the divergence of the velocity field. A respective
4615	!-- correction is needed to overcome numerical instabilities caused
4616	!-- by a not sufficient reduction of divergences near topography.
4617	div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
4618	+ ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
4619	+ ( w_comp * rho_air_zw(k) - &
4620	( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) &
4621	) * drho_air(k) * ddzw(k) &
4622	) * 0.5_wp
4623
4624	tend(k,j,i) = tend(k,j,i) - ( &
4625	( flux_r(k) + diss_r(k) &
4626	- swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j) &
4627	) * ddx &
4628	+ ( flux_n(k) + diss_n(k) &
4629	- swap_flux_y_local_v(k) - swap_diss_y_local_v(k) &
4630	) * ddy &
4631	+ ( ( flux_t(k) + diss_t(k) ) &
4632	- ( flux_d + diss_d ) &
4633	) * drho_air(k) * ddzw(k) &
4634	) + v(k,j,i) * div
4635
4636	swap_flux_x_local_v(k,j) = flux_r(k)
4637	swap_diss_x_local_v(k,j) = diss_r(k)
4638	swap_flux_y_local_v(k) = flux_n(k)
4639	swap_diss_y_local_v(k) = diss_n(k)
4640	flux_d = flux_t(k)
4641	diss_d = diss_t(k)
4642
4643	!
4644	!-- Statistical Evaluation of v'v'. The factor has to be applied
4645	!-- for right evaluation when gallilei_trans = .T. .
4646	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
4647	+ ( flux_n(k) &
4648	* ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
4649	/ ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &
4650	+ diss_n(k) &
4651	* ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
4652	/ ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &
4653	) * weight_substep(intermediate_timestep_count)
4654	!
4655	!-- Statistical Evaluation of w'u'.
4656	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
4657	+ ( flux_t(k) &
4658	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4659	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
4660	+ diss_t(k) &
4661	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4662	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
4663	) * weight_substep(intermediate_timestep_count)
4664
4665	ENDDO
4666	ENDDO
4667	ENDDO
4668	sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
4669
4670
4671	END SUBROUTINE advec_v_ws
4672
4673
4674	!------------------------------------------------------------------------------!
4675	! Description:
4676	! ------------
4677	!> Advection of w - Call for all grid points
4678	!------------------------------------------------------------------------------!
4679	SUBROUTINE advec_w_ws
4680
4681	USE arrays_3d, &
4682	ONLY: ddzu, drho_air_zw, tend, u, v, w, rho_air
4683
4684	USE constants, &
4685	ONLY: adv_mom_1, adv_mom_3, adv_mom_5
4686
4687	USE control_parameters, &
4688	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
4689
4690	USE grid_variables, &
4691	ONLY: ddx, ddy
4692
4693	USE indices, &
4694	ONLY: nxl, nxr, nyn, nys, nzb, nzb_max, nzt, advc_flags_1, &
4695	advc_flags_2
4696
4697	USE kinds
4698
4699	USE statistics, &
4700	ONLY: hom, sums_ws2_ws_l, weight_substep
4701
4702	IMPLICIT NONE
4703
4704	INTEGER(iwp) :: i !<
4705	INTEGER(iwp) :: ibit27 !<
4706	INTEGER(iwp) :: ibit28 !<
4707	INTEGER(iwp) :: ibit29 !<
4708	INTEGER(iwp) :: ibit30 !<
4709	INTEGER(iwp) :: ibit31 !<
4710	INTEGER(iwp) :: ibit32 !<
4711	INTEGER(iwp) :: ibit33 !<
4712	INTEGER(iwp) :: ibit34 !<
4713	INTEGER(iwp) :: ibit35 !<
4714	INTEGER(iwp) :: j !<
4715	INTEGER(iwp) :: k !<
4716	INTEGER(iwp) :: k_mm !<
4717	INTEGER(iwp) :: k_pp !<
4718	INTEGER(iwp) :: k_ppp !<
4719	INTEGER(iwp) :: tn = 0 !<
4720
4721	REAL(wp) :: diss_d !<
4722	REAL(wp) :: div !<
4723	REAL(wp) :: flux_d !<
4724	REAL(wp) :: gu !<
4725	REAL(wp) :: gv !<
4726	REAL(wp) :: u_comp !<
4727	REAL(wp) :: v_comp !<
4728	REAL(wp) :: w_comp !<
4729
4730	REAL(wp), DIMENSION(nzb:nzt) :: diss_t !<
4731	REAL(wp), DIMENSION(nzb:nzt) :: flux_t !<
4732
4733	REAL(wp), DIMENSION(nzb+1:nzt) :: diss_n !<
4734	REAL(wp), DIMENSION(nzb+1:nzt) :: diss_r !<
4735	REAL(wp), DIMENSION(nzb+1:nzt) :: flux_n !<
4736	REAL(wp), DIMENSION(nzb+1:nzt) :: flux_r !<
4737	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local_w !<
4738	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local_w !<
4739
4740	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_w !<
4741	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_w !<
4742
4743	gu = 2.0_wp * u_gtrans
4744	gv = 2.0_wp * v_gtrans
4745	!
4746	!-- compute the whole left boundary of the processor domain
4747	i = nxl
4748	DO j = nys, nyn
4749	DO k = nzb+1, nzb_max
4750
4751	ibit29 = IBITS(advc_flags_1(k,j,i-1),29,1)
4752	ibit28 = IBITS(advc_flags_1(k,j,i-1),28,1)
4753	ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)
4754
4755	u_comp = u(k+1,j,i) + u(k,j,i) - gu
4756	swap_flux_x_local_w(k,j) = u_comp * ( &
4757	( 37.0_wp * ibit29 * adv_mom_5 &
4758	+ 7.0_wp * ibit28 * adv_mom_3 &
4759	+ ibit27 * adv_mom_1 &
4760	) * &
4761	( w(k,j,i) + w(k,j,i-1) ) &
4762	- ( 8.0_wp * ibit29 * adv_mom_5 &
4763	+ ibit28 * adv_mom_3 &
4764	) * &
4765	( w(k,j,i+1) + w(k,j,i-2) ) &
4766	+ ( ibit29 * adv_mom_5 &
4767	) * &
4768	( w(k,j,i+2) + w(k,j,i-3) ) &
4769	)
4770
4771	swap_diss_x_local_w(k,j) = - ABS( u_comp ) * ( &
4772	( 10.0_wp * ibit29 * adv_mom_5 &
4773	+ 3.0_wp * ibit28 * adv_mom_3 &
4774	+ ibit27 * adv_mom_1 &
4775	) * &
4776	( w(k,j,i) - w(k,j,i-1) ) &
4777	- ( 5.0_wp * ibit29 * adv_mom_5 &
4778	+ ibit28 * adv_mom_3 &
4779	) * &
4780	( w(k,j,i+1) - w(k,j,i-2) ) &
4781	+ ( ibit29 * adv_mom_5 &
4782	) * &
4783	( w(k,j,i+2) - w(k,j,i-3) ) &
4784	)
4785
4786	ENDDO
4787
4788	DO k = nzb_max+1, nzt
4789
4790	u_comp = u(k+1,j,i) + u(k,j,i) - gu
4791	swap_flux_x_local_w(k,j) = u_comp * ( &
4792	37.0_wp * ( w(k,j,i) + w(k,j,i-1) ) &
4793	- 8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) ) &
4794	+ ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
4795	swap_diss_x_local_w(k,j) = - ABS( u_comp ) * ( &
4796	10.0_wp * ( w(k,j,i) - w(k,j,i-1) ) &
4797	- 5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) ) &
4798	+ ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
4799
4800	ENDDO
4801
4802	ENDDO
4803
4804	DO i = nxl, nxr
4805
4806	j = nys
4807	DO k = nzb+1, nzb_max
4808
4809	ibit32 = IBITS(advc_flags_2(k,j-1,i),0,1)
4810	ibit31 = IBITS(advc_flags_1(k,j-1,i),31,1)
4811	ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)
4812
4813	v_comp = v(k+1,j,i) + v(k,j,i) - gv
4814	swap_flux_y_local_w(k) = v_comp * ( &
4815	( 37.0_wp * ibit32 * adv_mom_5 &
4816	+ 7.0_wp * ibit31 * adv_mom_3 &
4817	+ ibit30 * adv_mom_1 &
4818	) * &
4819	( w(k,j,i) + w(k,j-1,i) ) &
4820	- ( 8.0_wp * ibit32 * adv_mom_5 &
4821	+ ibit31 * adv_mom_3 &
4822	) * &
4823	( w(k,j+1,i) + w(k,j-2,i) ) &
4824	+ ( ibit32 * adv_mom_5 &
4825	) * &
4826	( w(k,j+2,i) + w(k,j-3,i) ) &
4827	)
4828
4829	swap_diss_y_local_w(k) = - ABS( v_comp ) * ( &
4830	( 10.0_wp * ibit32 * adv_mom_5 &
4831	+ 3.0_wp * ibit31 * adv_mom_3 &
4832	+ ibit30 * adv_mom_1 &
4833	) * &
4834	( w(k,j,i) - w(k,j-1,i) ) &
4835	- ( 5.0_wp * ibit32 * adv_mom_5 &
4836	+ ibit31 * adv_mom_3 &
4837	) * &
4838	( w(k,j+1,i) - w(k,j-2,i) ) &
4839	+ ( ibit32 * adv_mom_5 &
4840	) * &
4841	( w(k,j+2,i) - w(k,j-3,i) ) &
4842	)
4843
4844	ENDDO
4845
4846	DO k = nzb_max+1, nzt
4847
4848	v_comp = v(k+1,j,i) + v(k,j,i) - gv
4849	swap_flux_y_local_w(k) = v_comp * ( &
4850	37.0_wp * ( w(k,j,i) + w(k,j-1,i) ) &
4851	- 8.0_wp * ( w(k,j+1,i) +w(k,j-2,i) ) &
4852	+ ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
4853	swap_diss_y_local_w(k) = - ABS( v_comp ) * ( &
4854	10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) &
4855	- 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) &
4856	+ ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
4857
4858	ENDDO
4859
4860	DO j = nys, nyn
4861
4862	!
4863	!-- The lower flux has to be calculated explicetely for the tendency
4864	!-- at the first w-level. For topography wall this is done implicitely
4865	!-- by advc_flags_1.
4866	k = nzb + 1
4867	w_comp = w(k,j,i) + w(k-1,j,i)
4868	flux_t(0) = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
4869	diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
4870	flux_d = flux_t(0)
4871	diss_d = diss_t(0)
4872
4873	DO k = nzb+1, nzb_max
4874
4875	ibit29 = IBITS(advc_flags_1(k,j,i),29,1)
4876	ibit28 = IBITS(advc_flags_1(k,j,i),28,1)
4877	ibit27 = IBITS(advc_flags_1(k,j,i),27,1)
4878
4879	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
4880	flux_r(k) = u_comp * ( &
4881	( 37.0_wp * ibit29 * adv_mom_5 &
4882	+ 7.0_wp * ibit28 * adv_mom_3 &
4883	+ ibit27 * adv_mom_1 &
4884	) * &
4885	( w(k,j,i+1) + w(k,j,i) ) &
4886	- ( 8.0_wp * ibit29 * adv_mom_5 &
4887	+ ibit28 * adv_mom_3 &
4888	) * &
4889	( w(k,j,i+2) + w(k,j,i-1) ) &
4890	+ ( ibit29 * adv_mom_5 &
4891	) * &
4892	( w(k,j,i+3) + w(k,j,i-2) ) &
4893	)
4894
4895	diss_r(k) = - ABS( u_comp ) * ( &
4896	( 10.0_wp * ibit29 * adv_mom_5 &
4897	+ 3.0_wp * ibit28 * adv_mom_3 &
4898	+ ibit27 * adv_mom_1 &
4899	) * &
4900	( w(k,j,i+1) - w(k,j,i) ) &
4901	- ( 5.0_wp * ibit29 * adv_mom_5 &
4902	+ ibit28 * adv_mom_3 &
4903	) * &
4904	( w(k,j,i+2) - w(k,j,i-1) ) &
4905	+ ( ibit29 * adv_mom_5 &
4906	) * &
4907	( w(k,j,i+3) - w(k,j,i-2) ) &
4908	)
4909
4910	ibit32 = IBITS(advc_flags_2(k,j,i),0,1)
4911	ibit31 = IBITS(advc_flags_1(k,j,i),31,1)
4912	ibit30 = IBITS(advc_flags_1(k,j,i),30,1)
4913
4914	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
4915	flux_n(k) = v_comp * ( &
4916	( 37.0_wp * ibit32 * adv_mom_5 &
4917	+ 7.0_wp * ibit31 * adv_mom_3 &
4918	+ ibit30 * adv_mom_1 &
4919	) * &
4920	( w(k,j+1,i) + w(k,j,i) ) &
4921	- ( 8.0_wp * ibit32 * adv_mom_5 &
4922	+ ibit31 * adv_mom_3 &
4923	) * &
4924	( w(k,j+2,i) + w(k,j-1,i) ) &
4925	+ ( ibit32 * adv_mom_5 &
4926	) * &
4927	( w(k,j+3,i) + w(k,j-2,i) ) &
4928	)
4929
4930	diss_n(k) = - ABS( v_comp ) * ( &
4931	( 10.0_wp * ibit32 * adv_mom_5 &
4932	+ 3.0_wp * ibit31 * adv_mom_3 &
4933	+ ibit30 * adv_mom_1 &
4934	) * &
4935	( w(k,j+1,i) - w(k,j,i) ) &
4936	- ( 5.0_wp * ibit32 * adv_mom_5 &
4937	+ ibit31 * adv_mom_3 &
4938	) * &
4939	( w(k,j+2,i) - w(k,j-1,i) ) &
4940	+ ( ibit32 * adv_mom_5 &
4941	) * &
4942	( w(k,j+3,i) - w(k,j-2,i) ) &
4943	)
4944	!
4945	!-- k index has to be modified near bottom and top, else array
4946	!-- subscripts will be exceeded.
4947	ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
4948	ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
4949	ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
4950
4951	k_ppp = k + 3 * ibit35
4952	k_pp = k + 2 * ( 1 - ibit33 )
4953	k_mm = k - 2 * ibit35
4954
4955	w_comp = w(k+1,j,i) + w(k,j,i)
4956	flux_t(k) = w_comp * rho_air(k+1) * ( &
4957	( 37.0_wp * ibit35 * adv_mom_5 &
4958	+ 7.0_wp * ibit34 * adv_mom_3 &
4959	+ ibit33 * adv_mom_1 &
4960	) * &
4961	( w(k+1,j,i) + w(k,j,i) ) &
4962	- ( 8.0_wp * ibit35 * adv_mom_5 &
4963	+ ibit34 * adv_mom_3 &
4964	) * &
4965	( w(k_pp,j,i) + w(k-1,j,i) ) &
4966	+ ( ibit35 * adv_mom_5 &
4967	) * &
4968	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
4969	)
4970
4971	diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * ( &
4972	( 10.0_wp * ibit35 * adv_mom_5 &
4973	+ 3.0_wp * ibit34 * adv_mom_3 &
4974	+ ibit33 * adv_mom_1 &
4975	) * &
4976	( w(k+1,j,i) - w(k,j,i) ) &
4977	- ( 5.0_wp * ibit35 * adv_mom_5 &
4978	+ ibit34 * adv_mom_3 &
4979	) * &
4980	( w(k_pp,j,i) - w(k-1,j,i) ) &
4981	+ ( ibit35 * adv_mom_5 &
4982	) * &
4983	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
4984	)
4985	!
4986	!-- Calculate the divergence of the velocity field. A respective
4987	!-- correction is needed to overcome numerical instabilities caused
4988	!-- by a not sufficient reduction of divergences near topography.
4989	div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 ) &
4990	- ( u(k+1,j,i) + u(k,j,i) ) &
4991	* ( IBITS(advc_flags_1(k,j,i-1),27,1) &
4992	+ IBITS(advc_flags_1(k,j,i-1),28,1) &
4993	+ IBITS(advc_flags_1(k,j,i-1),29,1) &
4994	) &
4995	) * ddx &
4996	+ ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) &
4997	- ( v(k+1,j,i) + v(k,j,i) ) &
4998	* ( IBITS(advc_flags_1(k,j-1,i),30,1) &
4999	+ IBITS(advc_flags_1(k,j-1,i),31,1) &
5000	+ IBITS(advc_flags_2(k,j-1,i),0,1) &
5001	) &
5002	) * ddy &
5003	+ ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 ) &
5004	- ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) &
5005	* ( IBITS(advc_flags_2(k-1,j,i),1,1) &
5006	+ IBITS(advc_flags_2(k-1,j,i),2,1) &
5007	+ IBITS(advc_flags_2(k-1,j,i),3,1) &
5008	) &
5009	) * drho_air_zw(k) * ddzu(k+1) &
5010	) * 0.5_wp
5011
5012
5013
5014	tend(k,j,i) = tend(k,j,i) - ( &
5015	( flux_r(k) + diss_r(k) &
5016	- swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j) &
5017	) * ddx &
5018	+ ( flux_n(k) + diss_n(k) &
5019	- swap_flux_y_local_w(k) - swap_diss_y_local_w(k) &
5020	) * ddy &
5021	+ ( ( flux_t(k) + diss_t(k) ) &
5022	- ( flux_d + diss_d ) &
5023	) * drho_air_zw(k) * ddzu(k+1) &
5024	) + div * w(k,j,i)
5025
5026	swap_flux_x_local_w(k,j) = flux_r(k)
5027	swap_diss_x_local_w(k,j) = diss_r(k)
5028	swap_flux_y_local_w(k) = flux_n(k)
5029	swap_diss_y_local_w(k) = diss_n(k)
5030	flux_d = flux_t(k)
5031	diss_d = diss_t(k)
5032
5033	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
5034	+ ( flux_t(k) &
5035	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
5036	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
5037	+ diss_t(k) &
5038	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
5039	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
5040	) * weight_substep(intermediate_timestep_count)
5041
5042	ENDDO
5043
5044	DO k = nzb_max+1, nzt
5045
5046	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
5047	flux_r(k) = u_comp * ( &
5048	37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) &
5049	- 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) &
5050	+ ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
5051
5052	diss_r(k) = - ABS( u_comp ) * ( &
5053	10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) &
5054	- 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) &
5055	+ ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
5056
5057	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
5058	flux_n(k) = v_comp * ( &
5059	37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) &
5060	- 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) &
5061	+ ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
5062
5063	diss_n(k) = - ABS( v_comp ) * ( &
5064	10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) &
5065	- 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) &
5066	+ ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
5067	!
5068	!-- k index has to be modified near bottom and top, else array
5069	!-- subscripts will be exceeded.
5070	ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
5071	ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
5072	ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
5073
5074	k_ppp = k + 3 * ibit35
5075	k_pp = k + 2 * ( 1 - ibit33 )
5076	k_mm = k - 2 * ibit35
5077
5078	w_comp = w(k+1,j,i) + w(k,j,i)
5079	flux_t(k) = w_comp * rho_air(k+1) * ( &
5080	( 37.0_wp * ibit35 * adv_mom_5 &
5081	+ 7.0_wp * ibit34 * adv_mom_3 &
5082	+ ibit33 * adv_mom_1 &
5083	) * &
5084	( w(k+1,j,i) + w(k,j,i) ) &
5085	- ( 8.0_wp * ibit35 * adv_mom_5 &
5086	+ ibit34 * adv_mom_3 &
5087	) * &
5088	( w(k_pp,j,i) + w(k-1,j,i) ) &
5089	+ ( ibit35 * adv_mom_5 &
5090	) * &
5091	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
5092	)
5093
5094	diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * ( &
5095	( 10.0_wp * ibit35 * adv_mom_5 &
5096	+ 3.0_wp * ibit34 * adv_mom_3 &
5097	+ ibit33 * adv_mom_1 &
5098	) * &
5099	( w(k+1,j,i) - w(k,j,i) ) &
5100	- ( 5.0_wp * ibit35 * adv_mom_5 &
5101	+ ibit34 * adv_mom_3 &
5102	) * &
5103	( w(k_pp,j,i) - w(k-1,j,i) ) &
5104	+ ( ibit35 * adv_mom_5 &
5105	) * &
5106	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
5107	)
5108	!
5109	!-- Calculate the divergence of the velocity field. A respective
5110	!-- correction is needed to overcome numerical instabilities caused
5111	!-- by a not sufficient reduction of divergences near topography.
5112	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
5113	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
5114	+ ( w_comp * rho_air(k+1) - &
5115	( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) &
5116	) * drho_air_zw(k) * ddzu(k+1) &
5117	) * 0.5_wp
5118
5119	tend(k,j,i) = tend(k,j,i) - ( &
5120	( flux_r(k) + diss_r(k) &
5121	- swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j) &
5122	) * ddx &
5123	+ ( flux_n(k) + diss_n(k) &
5124	- swap_flux_y_local_w(k) - swap_diss_y_local_w(k) &
5125	) * ddy &
5126	+ ( ( flux_t(k) + diss_t(k) ) &
5127	- ( flux_d + diss_d ) &
5128	) * drho_air_zw(k) * ddzu(k+1) &
5129	) + div * w(k,j,i)
5130
5131	swap_flux_x_local_w(k,j) = flux_r(k)
5132	swap_diss_x_local_w(k,j) = diss_r(k)
5133	swap_flux_y_local_w(k) = flux_n(k)
5134	swap_diss_y_local_w(k) = diss_n(k)
5135	flux_d = flux_t(k)
5136	diss_d = diss_t(k)
5137
5138	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
5139	+ ( flux_t(k) &
5140	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
5141	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
5142	+ diss_t(k) &
5143	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
5144	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
5145	) * weight_substep(intermediate_timestep_count)
5146
5147	ENDDO
5148	ENDDO
5149	ENDDO
5150
5151	END SUBROUTINE advec_w_ws
5152
5153	END MODULE advec_ws

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |