Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |