Home

Context Navigation

advec_ws.f90 @ 4204

Last change on this file since 4204 was 4204, checked in by knoop, 5 years ago
Bugfix in advec_s_ws: Changed sk_num initialization default to avoid implicit SAVE-attribut.
Property svn:keywords set to `Id`
File size: 374.6 KB

Line
1	!> @file advec_ws.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
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-2019 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: advec_ws.f90 4204 2019-08-30 12:30:17Z knoop $
27	! Bugfix: Changed sk_num initialization default to avoid implicit SAVE-Attribut
28	!
29	! 4182 2019-08-22 15:20:23Z scharf
30	! Corrected "Former revisions" section
31	!
32	! 4110 2019-07-22 17:05:21Z suehring
33	! - Separate initialization of advection flags for momentum and scalars. In this
34	! context, resort the bits and do some minor formatting.
35	! - Make flag initialization for scalars more flexible, introduce an
36	! arguemnt list to indicate non-cyclic boundaries (required for decycled
37	! scalars such as chemical species or aerosols)
38	! - Introduce extended 'degradation zones', where horizontal advection of
39	! passive scalars is discretized by first-order scheme at all grid points
40	! that in the vicinity of buildings (<= 3 grid points). Even though no
41	! building is within the numerical stencil, first-order scheme is used.
42	! At fourth and fifth grid point the order of the horizontal advection scheme
43	! is successively upgraded.
44	! These extended degradation zones are used to avoid stationary numerical
45	! oscillations, which are responsible for high concentration maxima that may
46	! appear under shear-free stable conditions.
47	! - Change interface for scalar advection routine.
48	! - Bugfix, avoid uninitialized value sk_num in vector version of scalar
49	! advection
50	!
51	! 4109 2019-07-22 17:00:34Z suehring
52	! Implementation of a flux limiter according to Skamarock (2006) for the
53	! vertical scalar advection. Please note, this is only implemented for the
54	! cache-optimized version at the moment. Implementation for the vector-
55	! optimized version will follow after critical issues concerning
56	! vectorization are fixed.
57	!
58	! 3873 2019-04-08 15:44:30Z knoop
59	! Moved ocean_mode specific code to ocean_mod
60	!
61	! 3872 2019-04-08 15:03:06Z knoop
62	! Moved all USE statements to module level + removed salsa dependency
63	!
64	! 3871 2019-04-08 14:38:39Z knoop
65	! Moving initialization of bcm specific flux arrays into bulk_cloud_model_mod
66	!
67	! 3864 2019-04-05 09:01:56Z monakurppa
68	! Remove tailing white spaces
69	!
70	! 3696 2019-01-24 16:37:35Z suehring
71	! Bugfix in degradation height
72	!
73	! 3661 2019-01-08 18:22:50Z suehring
74	! - Minor bugfix in divergence correction (only has implications at
75	! downward-facing wall surfaces)
76	! - Remove setting of Neumann condition for horizontal velocity variances
77	! - Split loops for tendency calculation and divergence correction in order to
78	! reduce bit queries
79	! - Introduce new parameter nzb_max_l to better control order degradation at
80	! non-cyclic boundaries
81	!
82	! 3655 2019-01-07 16:51:22Z knoop
83	! OpenACC port for SPEC
84	!
85	! 411 2009-12-11 12:31:43 Z suehring
86	! Initial revision
87	!
88	! Authors:
89	! --------
90	! @author Matthias Suehring
91	!
92	!
93	! Description:
94	! ------------
95	!> Advection scheme for scalars and momentum using the flux formulation of
96	!> Wicker and Skamarock 5th order. Additionally the module contains of a
97	!> routine using for initialisation and steering of the statical evaluation.
98	!> The computation of turbulent fluxes takes place inside the advection
99	!> routines.
100	!> Near non-cyclic boundaries the order of the applied advection scheme is
101	!> degraded.
102	!> A divergence correction is applied. It is necessary for topography, since
103	!> the divergence is not sufficiently reduced, resulting in erroneous fluxes
104	!> and could lead to numerical instabilities.
105	!>
106	!> @todo Implement monotonic flux limiter also for vector version.
107	!> @todo Move 3d arrays advc_flag, advc_flags_m from modules to advec_ws
108	!> @todo Move arrays flux_l_x from modules to advec_ws
109	!------------------------------------------------------------------------------!
110	MODULE advec_ws
111
112	USE arrays_3d, &
113	ONLY: ddzu, ddzw, tend, u, v, w, &
114	drho_air, drho_air_zw, rho_air, rho_air_zw, &
115	u_stokes_zu, v_stokes_zu, &
116	diss_l_diss, diss_l_e, diss_l_pt, diss_l_q, &
117	diss_l_s, diss_l_sa, diss_l_u, diss_l_v, diss_l_w, &
118	flux_l_diss, flux_l_e, flux_l_pt, flux_l_q, flux_l_s, &
119	flux_l_sa, flux_l_u, flux_l_v, flux_l_w, &
120	diss_s_diss, diss_s_e, diss_s_pt, diss_s_q, diss_s_s, &
121	diss_s_sa, diss_s_u, diss_s_v, diss_s_w, &
122	flux_s_diss, flux_s_e, flux_s_pt, flux_s_q, flux_s_s, &
123	flux_s_sa, flux_s_u, flux_s_v, flux_s_w
124
125	USE control_parameters, &
126	ONLY: air_chemistry, &
127	bc_dirichlet_l, &
128	bc_dirichlet_n, &
129	bc_dirichlet_r, &
130	bc_dirichlet_s, &
131	bc_radiation_l, &
132	bc_radiation_n, &
133	bc_radiation_r, &
134	bc_radiation_s, &
135	humidity, &
136	loop_optimization, &
137	passive_scalar, &
138	rans_tke_e, &
139	momentum_advec, &
140	salsa, &
141	scalar_advec, &
142	intermediate_timestep_count, &
143	u_gtrans, &
144	v_gtrans, &
145	ws_scheme_mom, &
146	ws_scheme_sca, &
147	dt_3d
148
149	USE indices, &
150	ONLY: advc_flags_m, &
151	advc_flags_s, &
152	nbgp, &
153	nx, &
154	nxl, &
155	nxlg, &
156	nxlu, &
157	nxr, &
158	nxrg, &
159	ny, &
160	nyn, &
161	nyng, &
162	nys, &
163	nysg, &
164	nysv, &
165	nzb, &
166	nzb_max, &
167	nzt, &
168	wall_flags_0
169
170	USE grid_variables, &
171	ONLY: ddx, ddy
172
173	USE pegrid, &
174	ONLY: threads_per_task
175
176	USE kinds
177
178	USE statistics, &
179	ONLY: sums_salsa_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, &
180	sums_wspts_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l, &
181	sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l, &
182	sums_wsqcs_ws_l, sums_wsqrs_ws_l, &
183	sums_wsncs_ws_l, sums_wsnrs_ws_l, &
184	hom, weight_substep
185
186	IMPLICIT NONE
187
188	REAL(wp) :: adv_mom_1 !< 1/4 - constant used in 5th-order advection scheme for momentum advection (1st-order part)
189	REAL(wp) :: adv_mom_3 !< 1/24 - constant used in 5th-order advection scheme for momentum advection (3rd-order part)
190	REAL(wp) :: adv_mom_5 !< 1/120 - constant used in 5th-order advection scheme for momentum advection (5th-order part)
191	REAL(wp) :: adv_sca_1 !< 1/2 - constant used in 5th-order advection scheme for scalar advection (1st-order part)
192	REAL(wp) :: adv_sca_3 !< 1/12 - constant used in 5th-order advection scheme for scalar advection (3rd-order part)
193	REAL(wp) :: adv_sca_5 !< 1/60 - constant used in 5th-order advection scheme for scalar advection (5th-order part)
194
195	PRIVATE
196	PUBLIC advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, ws_init, &
197	ws_init_flags_momentum, ws_init_flags_scalar, ws_statistics
198
199	INTERFACE ws_init
200	MODULE PROCEDURE ws_init
201	END INTERFACE ws_init
202
203	INTERFACE ws_init_flags_momentum
204	MODULE PROCEDURE ws_init_flags_momentum
205	END INTERFACE ws_init_flags_momentum
206
207	INTERFACE ws_init_flags_scalar
208	MODULE PROCEDURE ws_init_flags_scalar
209	END INTERFACE ws_init_flags_scalar
210
211	INTERFACE ws_statistics
212	MODULE PROCEDURE ws_statistics
213	END INTERFACE ws_statistics
214
215	INTERFACE advec_s_ws
216	MODULE PROCEDURE advec_s_ws
217	MODULE PROCEDURE advec_s_ws_ij
218	END INTERFACE advec_s_ws
219
220	INTERFACE advec_u_ws
221	MODULE PROCEDURE advec_u_ws
222	MODULE PROCEDURE advec_u_ws_ij
223	END INTERFACE advec_u_ws
224
225	INTERFACE advec_v_ws
226	MODULE PROCEDURE advec_v_ws
227	MODULE PROCEDURE advec_v_ws_ij
228	END INTERFACE advec_v_ws
229
230	INTERFACE advec_w_ws
231	MODULE PROCEDURE advec_w_ws
232	MODULE PROCEDURE advec_w_ws_ij
233	END INTERFACE advec_w_ws
234
235	CONTAINS
236
237
238	!------------------------------------------------------------------------------!
239	! Description:
240	! ------------
241	!> Initialization of WS-scheme
242	!------------------------------------------------------------------------------!
243	SUBROUTINE ws_init
244
245	!
246	!-- Set the appropriate factors for scalar and momentum advection.
247	adv_sca_5 = 1.0_wp / 60.0_wp
248	adv_sca_3 = 1.0_wp / 12.0_wp
249	adv_sca_1 = 1.0_wp / 2.0_wp
250	adv_mom_5 = 1.0_wp / 120.0_wp
251	adv_mom_3 = 1.0_wp / 24.0_wp
252	adv_mom_1 = 1.0_wp / 4.0_wp
253	!
254	!-- Arrays needed for statical evaluation of fluxes.
255	IF ( ws_scheme_mom ) THEN
256
257	ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1), &
258	sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1), &
259	sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
260	sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
261	sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
262
263	sums_wsus_ws_l = 0.0_wp
264	sums_wsvs_ws_l = 0.0_wp
265	sums_us2_ws_l = 0.0_wp
266	sums_vs2_ws_l = 0.0_wp
267	sums_ws2_ws_l = 0.0_wp
268
269	ENDIF
270
271	IF ( ws_scheme_sca ) THEN
272
273	ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
274	sums_wspts_ws_l = 0.0_wp
275
276	IF ( humidity ) THEN
277	ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
278	sums_wsqs_ws_l = 0.0_wp
279	ENDIF
280
281	IF ( passive_scalar ) THEN
282	ALLOCATE( sums_wsss_ws_l(nzb:nzt+1,0:threads_per_task-1) )
283	sums_wsss_ws_l = 0.0_wp
284	ENDIF
285
286	ENDIF
287
288	!
289	!-- Arrays needed for reasons of speed optimization for cache version.
290	!-- For the vector version the buffer arrays are not necessary,
291	!-- because the the fluxes can swapped directly inside the loops of the
292	!-- advection routines.
293	IF ( loop_optimization /= 'vector' ) THEN
294
295	IF ( ws_scheme_mom ) THEN
296
297	ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1), &
298	flux_s_v(nzb+1:nzt,0:threads_per_task-1), &
299	flux_s_w(nzb+1:nzt,0:threads_per_task-1), &
300	diss_s_u(nzb+1:nzt,0:threads_per_task-1), &
301	diss_s_v(nzb+1:nzt,0:threads_per_task-1), &
302	diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
303	ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
304	flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
305	flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
306	diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
307	diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
308	diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
309
310	ENDIF
311
312	IF ( ws_scheme_sca ) THEN
313
314	ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1), &
315	flux_s_e(nzb+1:nzt,0:threads_per_task-1), &
316	diss_s_pt(nzb+1:nzt,0:threads_per_task-1), &
317	diss_s_e(nzb+1:nzt,0:threads_per_task-1) )
318	ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
319	flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
320	diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
321	diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
322
323	IF ( rans_tke_e ) THEN
324	ALLOCATE( flux_s_diss(nzb+1:nzt,0:threads_per_task-1), &
325	diss_s_diss(nzb+1:nzt,0:threads_per_task-1) )
326	ALLOCATE( flux_l_diss(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
327	diss_l_diss(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
328	ENDIF
329
330	IF ( humidity ) THEN
331	ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1), &
332	diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
333	ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
334	diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
335	ENDIF
336
337	IF ( passive_scalar ) THEN
338	ALLOCATE( flux_s_s(nzb+1:nzt,0:threads_per_task-1), &
339	diss_s_s(nzb+1:nzt,0:threads_per_task-1) )
340	ALLOCATE( flux_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
341	diss_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
342	ENDIF
343
344	ENDIF
345
346	ENDIF
347
348	END SUBROUTINE ws_init
349
350	!------------------------------------------------------------------------------!
351	! Description:
352	! ------------
353	!> Initialization of flags to control the order of the advection scheme near
354	!> solid walls and non-cyclic inflow boundaries, where the order is sucessively
355	!> degraded.
356	!------------------------------------------------------------------------------!
357	SUBROUTINE ws_init_flags_momentum
358
359
360	INTEGER(iwp) :: i !< index variable along x
361	INTEGER(iwp) :: j !< index variable along y
362	INTEGER(iwp) :: k !< index variable along z
363	INTEGER(iwp) :: k_mm !< dummy index along z
364	INTEGER(iwp) :: k_pp !< dummy index along z
365	INTEGER(iwp) :: k_ppp !< dummy index along z
366
367	LOGICAL :: flag_set !< steering variable for advection flags
368
369	advc_flags_m = 0
370
371	!
372	!-- Set advc_flags_m to steer the degradation of the advection scheme in advec_ws
373	!-- near topography, inflow- and outflow boundaries as well as bottom and
374	!-- top of model domain. advc_flags_m remains zero for all non-prognostic
375	!-- grid points.
376	!-- u-component
377	DO i = nxl, nxr
378	DO j = nys, nyn
379	DO k = nzb+1, nzt
380	!
381	!-- At first, set flags to WS1.
382	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
383	!-- in order to handle the left/south flux.
384	!-- near vertical walls.
385	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 )
386	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 )
387	!
388	!-- u component - x-direction
389	!-- WS1 (0), WS3 (1), WS5 (2)
390	IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),1) .OR. &
391	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
392	.AND. i <= nxlu ) .OR. &
393	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
394	.AND. i == nxr ) ) &
395	THEN
396	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 )
397	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),1) .AND. &
398	BTEST(wall_flags_0(k,j,i+1),1) .OR. &
399	.NOT. BTEST(wall_flags_0(k,j,i-1),1) ) &
400	.OR. &
401	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
402	.AND. i == nxr-1 ) .OR. &
403	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
404	.AND. i == nxlu+1) ) &
405	THEN
406	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 1 )
407	!
408	!-- Clear flag for WS1
409	advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 )
410	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),1) .AND. &
411	BTEST(wall_flags_0(k,j,i+2),1) .AND. &
412	BTEST(wall_flags_0(k,j,i-1),1) ) &
413	THEN
414	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 2 )
415	!
416	!-- Clear flag for WS1
417	advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 )
418	ENDIF
419	!
420	!-- u component - y-direction
421	!-- WS1 (3), WS3 (4), WS5 (5)
422	IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),1) .OR. &
423	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
424	.AND. j == nys ) .OR. &
425	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
426	.AND. j == nyn ) ) &
427	THEN
428	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 )
429	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),1) .AND. &
430	BTEST(wall_flags_0(k,j+1,i),1) .OR. &
431	.NOT. BTEST(wall_flags_0(k,j-1,i),1) ) &
432	.OR. &
433	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
434	.AND. j == nysv ) .OR. &
435	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
436	.AND. j == nyn-1 ) ) &
437	THEN
438	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 4 )
439	!
440	!-- Clear flag for WS1
441	advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 )
442	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),1) .AND. &
443	BTEST(wall_flags_0(k,j+2,i),1) .AND. &
444	BTEST(wall_flags_0(k,j-1,i),1) ) &
445	THEN
446	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 5 )
447	!
448	!-- Clear flag for WS1
449	advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 )
450	ENDIF
451	!
452	!-- u component - z-direction
453	!-- WS1 (6), WS3 (7), WS5 (8)
454	IF ( k == nzb+1 ) THEN
455	k_mm = nzb
456	ELSE
457	k_mm = k - 2
458	ENDIF
459	IF ( k > nzt-1 ) THEN
460	k_pp = nzt+1
461	ELSE
462	k_pp = k + 2
463	ENDIF
464	IF ( k > nzt-2 ) THEN
465	k_ppp = nzt+1
466	ELSE
467	k_ppp = k + 3
468	ENDIF
469
470	flag_set = .FALSE.
471	IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),1) .AND. &
472	BTEST(wall_flags_0(k,j,i),1) .OR. &
473	.NOT. BTEST(wall_flags_0(k_pp,j,i),1) .AND. &
474	BTEST(wall_flags_0(k,j,i),1) .OR. &
475	k == nzt ) &
476	THEN
477	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 6 )
478	flag_set = .TRUE.
479	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),1) .OR. &
480	.NOT. BTEST(wall_flags_0(k_ppp,j,i),1) ) .AND. &
481	BTEST(wall_flags_0(k-1,j,i),1) .AND. &
482	BTEST(wall_flags_0(k,j,i),1) .AND. &
483	BTEST(wall_flags_0(k+1,j,i),1) .AND. &
484	BTEST(wall_flags_0(k_pp,j,i),1) .OR. &
485	k == nzt - 1 ) &
486	THEN
487	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 )
488	flag_set = .TRUE.
489	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),1) .AND. &
490	BTEST(wall_flags_0(k-1,j,i),1) .AND. &
491	BTEST(wall_flags_0(k,j,i),1) .AND. &
492	BTEST(wall_flags_0(k+1,j,i),1) .AND. &
493	BTEST(wall_flags_0(k_pp,j,i),1) .AND. &
494	BTEST(wall_flags_0(k_ppp,j,i),1) .AND. &
495	.NOT. flag_set ) &
496	THEN
497	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 8 )
498	ENDIF
499
500	ENDDO
501	ENDDO
502	ENDDO
503	!
504	!-- v-component
505	DO i = nxl, nxr
506	DO j = nys, nyn
507	DO k = nzb+1, nzt
508	!
509	!-- At first, set flags to WS1.
510	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
511	!-- in order to handle the left/south flux.
512	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 )
513	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 )
514	!
515	!-- v component - x-direction
516	!-- WS1 (9), WS3 (10), WS5 (11)
517	IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),2) .OR. &
518	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
519	.AND. i == nxl ) .OR. &
520	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
521	.AND. i == nxr ) ) &
522	THEN
523	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 )
524	!
525	!-- WS3
526	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),2) .AND. &
527	BTEST(wall_flags_0(k,j,i+1),2) ) .OR. &
528	.NOT. BTEST(wall_flags_0(k,j,i-1),2) &
529	.OR. &
530	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
531	.AND. i == nxr-1 ) .OR. &
532	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
533	.AND. i == nxlu ) ) &
534	THEN
535	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 10 )
536	!
537	!-- Clear flag for WS1
538	advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 )
539	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),2) .AND. &
540	BTEST(wall_flags_0(k,j,i+2),2) .AND. &
541	BTEST(wall_flags_0(k,j,i-1),2) ) &
542	THEN
543	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 11 )
544	!
545	!-- Clear flag for WS1
546	advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 )
547	ENDIF
548	!
549	!-- v component - y-direction
550	!-- WS1 (12), WS3 (13), WS5 (14)
551	IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),2) .OR. &
552	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
553	.AND. j <= nysv ) .OR. &
554	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
555	.AND. j == nyn ) ) &
556	THEN
557	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 )
558	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),2) .AND. &
559	BTEST(wall_flags_0(k,j+1,i),2) .OR. &
560	.NOT. BTEST(wall_flags_0(k,j-1,i),2) ) &
561	.OR. &
562	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
563	.AND. j == nysv+1) .OR. &
564	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
565	.AND. j == nyn-1 ) ) &
566	THEN
567	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 13 )
568	!
569	!-- Clear flag for WS1
570	advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 )
571	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),2) .AND. &
572	BTEST(wall_flags_0(k,j+2,i),2) .AND. &
573	BTEST(wall_flags_0(k,j-1,i),2) ) &
574	THEN
575	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 14 )
576	!
577	!-- Clear flag for WS1
578	advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 )
579	ENDIF
580	!
581	!-- v component - z-direction
582	!-- WS1 (15), WS3 (16), WS5 (17)
583	IF ( k == nzb+1 ) THEN
584	k_mm = nzb
585	ELSE
586	k_mm = k - 2
587	ENDIF
588	IF ( k > nzt-1 ) THEN
589	k_pp = nzt+1
590	ELSE
591	k_pp = k + 2
592	ENDIF
593	IF ( k > nzt-2 ) THEN
594	k_ppp = nzt+1
595	ELSE
596	k_ppp = k + 3
597	ENDIF
598
599	flag_set = .FALSE.
600	IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),2) .AND. &
601	BTEST(wall_flags_0(k,j,i),2) .OR. &
602	.NOT. BTEST(wall_flags_0(k_pp,j,i),2) .AND. &
603	BTEST(wall_flags_0(k,j,i),2) .OR. &
604	k == nzt ) &
605	THEN
606	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 15 )
607	flag_set = .TRUE.
608	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),2) .OR. &
609	.NOT. BTEST(wall_flags_0(k_ppp,j,i),2) ) .AND. &
610	BTEST(wall_flags_0(k-1,j,i),2) .AND. &
611	BTEST(wall_flags_0(k,j,i),2) .AND. &
612	BTEST(wall_flags_0(k+1,j,i),2) .AND. &
613	BTEST(wall_flags_0(k_pp,j,i),2) .OR. &
614	k == nzt - 1 ) &
615	THEN
616	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 )
617	flag_set = .TRUE.
618	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),2) .AND. &
619	BTEST(wall_flags_0(k-1,j,i),2) .AND. &
620	BTEST(wall_flags_0(k,j,i),2) .AND. &
621	BTEST(wall_flags_0(k+1,j,i),2) .AND. &
622	BTEST(wall_flags_0(k_pp,j,i),2) .AND. &
623	BTEST(wall_flags_0(k_ppp,j,i),2) .AND. &
624	.NOT. flag_set ) &
625	THEN
626	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 17 )
627	ENDIF
628
629	ENDDO
630	ENDDO
631	ENDDO
632	!
633	!-- w - component
634	DO i = nxl, nxr
635	DO j = nys, nyn
636	DO k = nzb+1, nzt
637	!
638	!-- At first, set flags to WS1.
639	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
640	!-- in order to handle the left/south flux.
641	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 )
642	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 )
643	!
644	!-- w component - x-direction
645	!-- WS1 (18), WS3 (19), WS5 (20)
646	IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),3) .OR. &
647	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
648	.AND. i == nxl ) .OR. &
649	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
650	.AND. i == nxr ) ) &
651	THEN
652	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 )
653	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),3) .AND. &
654	BTEST(wall_flags_0(k,j,i+1),3) .OR. &
655	.NOT. BTEST(wall_flags_0(k,j,i-1),3) ) &
656	.OR. &
657	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
658	.AND. i == nxr-1 ) .OR. &
659	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
660	.AND. i == nxlu ) ) &
661	THEN
662	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 19 )
663	!
664	!-- Clear flag for WS1
665	advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 )
666	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),3) .AND. &
667	BTEST(wall_flags_0(k,j,i+2),3) .AND. &
668	BTEST(wall_flags_0(k,j,i-1),3) ) &
669	THEN
670	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i),20 )
671	!
672	!-- Clear flag for WS1
673	advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 )
674	ENDIF
675	!
676	!-- w component - y-direction
677	!-- WS1 (21), WS3 (22), WS5 (23)
678	IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),3) .OR. &
679	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
680	.AND. j == nys ) .OR. &
681	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
682	.AND. j == nyn ) ) &
683	THEN
684	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 )
685	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),3) .AND. &
686	BTEST(wall_flags_0(k,j+1,i),3) .OR. &
687	.NOT. BTEST(wall_flags_0(k,j-1,i),3) ) &
688	.OR. &
689	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
690	.AND. j == nysv ) .OR. &
691	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
692	.AND. j == nyn-1 ) ) &
693	THEN
694	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 22 )
695	!
696	!-- Clear flag for WS1
697	advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 )
698	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),3) .AND. &
699	BTEST(wall_flags_0(k,j+2,i),3) .AND. &
700	BTEST(wall_flags_0(k,j-1,i),3) ) &
701	THEN
702	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 23 )
703	!
704	!-- Clear flag for WS1
705	advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 )
706	ENDIF
707	!
708	!-- w component - z-direction
709	!-- WS1 (24), WS3 (25), WS5 (26)
710	flag_set = .FALSE.
711	IF ( k == nzb+1 ) THEN
712	k_mm = nzb
713	ELSE
714	k_mm = k - 2
715	ENDIF
716	IF ( k > nzt-1 ) THEN
717	k_pp = nzt+1
718	ELSE
719	k_pp = k + 2
720	ENDIF
721	IF ( k > nzt-2 ) THEN
722	k_ppp = nzt+1
723	ELSE
724	k_ppp = k + 3
725	ENDIF
726
727	IF ( ( .NOT. BTEST(wall_flags_0(k-1,j,i),3) .AND. &
728	.NOT. BTEST(wall_flags_0(k,j,i),3) .AND. &
729	BTEST(wall_flags_0(k+1,j,i),3) ) .OR. &
730	( .NOT. BTEST(wall_flags_0(k-1,j,i),3) .AND. &
731	BTEST(wall_flags_0(k,j,i),3) ) .OR. &
732	( .NOT. BTEST(wall_flags_0(k+1,j,i),3) .AND. &
733	BTEST(wall_flags_0(k,j,i),3) ) .OR. &
734	k == nzt ) &
735	THEN
736	!
737	!-- Please note, at k == nzb_w_inner(j,i) a flag is explictely
738	!-- set, although this is not a prognostic level. However,
739	!-- contrary to the advection of u,v and s this is necessary
740	!-- because flux_t(nzb_w_inner(j,i)) is used for the tendency
741	!-- at k == nzb_w_inner(j,i)+1.
742	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 24 )
743	flag_set = .TRUE.
744	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),3) .OR. &
745	.NOT. BTEST(wall_flags_0(k_ppp,j,i),3) ) .AND. &
746	BTEST(wall_flags_0(k-1,j,i),3) .AND. &
747	BTEST(wall_flags_0(k,j,i),3) .AND. &
748	BTEST(wall_flags_0(k+1,j,i),3) .OR. &
749	k == nzt - 1 ) &
750	THEN
751	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 25 )
752	flag_set = .TRUE.
753	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),3) .AND. &
754	BTEST(wall_flags_0(k-1,j,i),3) .AND. &
755	BTEST(wall_flags_0(k,j,i),3) .AND. &
756	BTEST(wall_flags_0(k+1,j,i),3) .AND. &
757	BTEST(wall_flags_0(k_pp,j,i),3) .AND. &
758	BTEST(wall_flags_0(k_ppp,j,i),3) .AND. &
759	.NOT. flag_set ) &
760	THEN
761	advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 26 )
762	ENDIF
763
764	ENDDO
765	ENDDO
766	ENDDO
767	!
768	!-- Exchange ghost points for advection flags
769	CALL exchange_horiz_int( advc_flags_m, nys, nyn, nxl, nxr, nzt, nbgp )
770	!
771	!-- Set boundary flags at inflow and outflow boundary in case of
772	!-- non-cyclic boundary conditions.
773	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
774	advc_flags_m(:,:,nxl-1) = advc_flags_m(:,:,nxl)
775	ENDIF
776
777	IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN
778	advc_flags_m(:,:,nxr+1) = advc_flags_m(:,:,nxr)
779	ENDIF
780
781	IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN
782	advc_flags_m(:,nyn+1,:) = advc_flags_m(:,nyn,:)
783	ENDIF
784
785	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
786	advc_flags_m(:,nys-1,:) = advc_flags_m(:,nys,:)
787	ENDIF
788
789	END SUBROUTINE ws_init_flags_momentum
790
791
792	!------------------------------------------------------------------------------!
793	! Description:
794	! ------------
795	!> Initialization of flags to control the order of the advection scheme near
796	!> solid walls and non-cyclic inflow boundaries, where the order is sucessively
797	!> degraded.
798	!------------------------------------------------------------------------------!
799	SUBROUTINE ws_init_flags_scalar( non_cyclic_l, non_cyclic_n, non_cyclic_r, &
800	non_cyclic_s, advc_flag, extensive_degrad )
801
802
803	INTEGER(iwp) :: i !< index variable along x
804	INTEGER(iwp) :: j !< index variable along y
805	INTEGER(iwp) :: k !< index variable along z
806	INTEGER(iwp) :: k_mm !< dummy index along z
807	INTEGER(iwp) :: k_pp !< dummy index along z
808	INTEGER(iwp) :: k_ppp !< dummy index along z
809
810	INTEGER(iwp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::&
811	advc_flag !< flag array to control order of scalar advection
812
813	LOGICAL :: flag_set !< steering variable for advection flags
814	LOGICAL :: non_cyclic_l !< flag that indicates non-cyclic boundary on the left
815	LOGICAL :: non_cyclic_n !< flag that indicates non-cyclic boundary on the north
816	LOGICAL :: non_cyclic_r !< flag that indicates non-cyclic boundary on the right
817	LOGICAL :: non_cyclic_s !< flag that indicates non-cyclic boundary on the south
818
819	LOGICAL, OPTIONAL :: extensive_degrad !< flag indicating that extensive degradation is required, e.g. for
820	!< passive scalars nearby topography along the horizontal directions,
821	!< as no monotonic limiter can be applied there
822	!
823	!-- Set flags to steer the degradation of the advection scheme in advec_ws
824	!-- near topography, inflow- and outflow boundaries as well as bottom and
825	!-- top of model domain. advc_flags_m remains zero for all non-prognostic
826	!-- grid points.
827	DO i = nxl, nxr
828	DO j = nys, nyn
829	DO k = nzb+1, nzt
830	IF ( .NOT. BTEST(wall_flags_0(k,j,i),0) ) CYCLE
831	!
832	!-- scalar - x-direction
833	!-- WS1 (0), WS3 (1), WS5 (2)
834	IF ( ( .NOT. BTEST(wall_flags_0(k,j,i+1),0) &
835	.OR. .NOT. BTEST(wall_flags_0(k,j,i+2),0) &
836	.OR. .NOT. BTEST(wall_flags_0(k,j,i-1),0) ) &
837	.OR. ( non_cyclic_l .AND. i == 0 ) &
838	.OR. ( non_cyclic_r .AND. i == nx ) ) THEN
839	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
840	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+3),0) &
841	.AND. BTEST(wall_flags_0(k,j,i+1),0) &
842	.AND. BTEST(wall_flags_0(k,j,i+2),0) &
843	.AND. BTEST(wall_flags_0(k,j,i-1),0) &
844	) .OR. &
845	( .NOT. BTEST(wall_flags_0(k,j,i-2),0) &
846	.AND. BTEST(wall_flags_0(k,j,i+1),0) &
847	.AND. BTEST(wall_flags_0(k,j,i+2),0) &
848	.AND. BTEST(wall_flags_0(k,j,i-1),0) &
849	) &
850	.OR. &
851	( non_cyclic_r .AND. i == nx-1 ) .OR. &
852	( non_cyclic_l .AND. i == 1 ) ) THEN
853	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
854	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),0) &
855	.AND. BTEST(wall_flags_0(k,j,i+2),0) &
856	.AND. BTEST(wall_flags_0(k,j,i+3),0) &
857	.AND. BTEST(wall_flags_0(k,j,i-1),0) &
858	.AND. BTEST(wall_flags_0(k,j,i-2),0) ) &
859	THEN
860	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 )
861	ENDIF
862	!
863	!-- scalar - y-direction
864	!-- WS1 (3), WS3 (4), WS5 (5)
865	IF ( ( .NOT. BTEST(wall_flags_0(k,j+1,i),0) &
866	.OR. .NOT. BTEST(wall_flags_0(k,j+2,i),0) &
867	.OR. .NOT. BTEST(wall_flags_0(k,j-1,i),0)) &
868	.OR. ( non_cyclic_s .AND. j == 0 ) &
869	.OR. ( non_cyclic_n .AND. j == ny ) ) THEN
870	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
871	!
872	!-- WS3
873	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+3,i),0) &
874	.AND. BTEST(wall_flags_0(k,j+1,i),0) &
875	.AND. BTEST(wall_flags_0(k,j+2,i),0) &
876	.AND. BTEST(wall_flags_0(k,j-1,i),0) &
877	) .OR. &
878	( .NOT. BTEST(wall_flags_0(k,j-2,i),0) &
879	.AND. BTEST(wall_flags_0(k,j+1,i),0) &
880	.AND. BTEST(wall_flags_0(k,j+2,i),0) &
881	.AND. BTEST(wall_flags_0(k,j-1,i),0) &
882	) &
883	.OR. &
884	( non_cyclic_s .AND. j == 1 ) .OR. &
885	( non_cyclic_n .AND. j == ny-1 ) ) THEN
886	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
887	!
888	!-- WS5
889	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),0) &
890	.AND. BTEST(wall_flags_0(k,j+2,i),0) &
891	.AND. BTEST(wall_flags_0(k,j+3,i),0) &
892	.AND. BTEST(wall_flags_0(k,j-1,i),0) &
893	.AND. BTEST(wall_flags_0(k,j-2,i),0) ) &
894	THEN
895	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 )
896	ENDIF
897	!
898	!-- Near topography, set horizontal advection scheme to 1st order
899	!-- for passive scalars, even if only one direction may be
900	!-- blocked by topography. These locations will be identified
901	!-- by wall_flags_0 bit 31. Note, since several modules define
902	!-- advection flags but may apply different scalar boundary
903	!-- conditions, bit 31 is temporarily stored on advc_flags.
904	!-- Moreover, note that this extended degradtion for passive
905	!-- scalars is not required for the vertical direction as there
906	!-- the monotonic limiter can be applied.
907	IF ( PRESENT( extensive_degrad ) ) THEN
908	IF ( extensive_degrad ) THEN
909	!
910	!-- At all grid points that are within a three-grid point
911	!-- range to topography, set 1st-order scheme.
912	IF( BTEST( advc_flag(k,j,i), 31 ) ) THEN
913	!
914	!-- Clear flags that might indicate higher-order
915	!-- advection along x- and y-direction.
916	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
917	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
918	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
919	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
920	!
921	!-- Set flags that indicate 1st-order advection along
922	!-- x- and y-direction.
923	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
924	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
925	ENDIF
926	!
927	!-- Adjacent to this extended degradation zone, successively
928	!-- upgrade the order of the scheme if this grid point isn't
929	!-- flagged with bit 31 (indicating extended degradation
930	!-- zone).
931	IF ( .NOT. BTEST( advc_flag(k,j,i), 31 ) ) THEN
932	!
933	!-- x-direction. First, clear all previous settings, than
934	!-- set flag for 3rd-order scheme.
935	IF ( BTEST( advc_flag(k,j,i-1), 31 ) .AND. &
936	BTEST( advc_flag(k,j,i+1), 31 ) ) THEN
937	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
938	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
939	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
940
941	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
942	ENDIF
943	!
944	!-- x-direction. First, clear all previous settings, than
945	!-- set flag for 5rd-order scheme.
946	IF ( .NOT. BTEST( advc_flag(k,j,i-1), 31 ) .AND. &
947	BTEST( advc_flag(k,j,i-2), 31 ) .AND. &
948	.NOT. BTEST( advc_flag(k,j,i+1), 31 ) .AND. &
949	BTEST( advc_flag(k,j,i+2), 31 ) ) THEN
950	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
951	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
952	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
953
954	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 )
955	ENDIF
956	!
957	!-- y-direction. First, clear all previous settings, than
958	!-- set flag for 3rd-order scheme.
959	IF ( BTEST( advc_flag(k,j-1,i), 31 ) .AND. &
960	BTEST( advc_flag(k,j+1,i), 31 ) ) THEN
961	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
962	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
963	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
964
965	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
966	ENDIF
967	!
968	!-- y-direction. First, clear all previous settings, than
969	!-- set flag for 5rd-order scheme.
970	IF ( .NOT. BTEST( advc_flag(k,j-1,i), 31 ) .AND. &
971	BTEST( advc_flag(k,j-2,i), 31 ) .AND. &
972	.NOT. BTEST( advc_flag(k,j+1,i), 31 ) .AND. &
973	BTEST( advc_flag(k,j+2,i), 31 ) ) THEN
974	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
975	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
976	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
977
978	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 )
979	ENDIF
980	ENDIF
981
982	ENDIF
983
984	!
985	!-- Near lateral boundary flags might be overwritten. Set
986	!-- them again.
987	!-- x-direction
988	IF ( ( non_cyclic_l .AND. i == 0 ) .OR. &
989	( non_cyclic_r .AND. i == nx ) ) THEN
990	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
991	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
992	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
993
994	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
995	ENDIF
996
997	IF ( ( non_cyclic_l .AND. i == 1 ) .OR. &
998	( non_cyclic_r .AND. i == nx-1 ) ) THEN
999	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
1000	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
1001	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
1002
1003	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
1004	ENDIF
1005	!
1006	!-- y-direction
1007	IF ( ( non_cyclic_n .AND. j == 0 ) .OR. &
1008	( non_cyclic_s .AND. j == ny ) ) THEN
1009	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
1010	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
1011	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
1012
1013	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
1014	ENDIF
1015
1016	IF ( ( non_cyclic_n .AND. j == 1 ) .OR. &
1017	( non_cyclic_s .AND. j == ny-1 ) ) THEN
1018	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
1019	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
1020	advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
1021
1022	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
1023	ENDIF
1024
1025	ENDIF
1026
1027
1028	!
1029	!-- scalar - z-direction
1030	!-- WS1 (6), WS3 (7), WS5 (8)
1031	IF ( k == nzb+1 ) THEN
1032	k_mm = nzb
1033	ELSE
1034	k_mm = k - 2
1035	ENDIF
1036	IF ( k > nzt-1 ) THEN
1037	k_pp = nzt+1
1038	ELSE
1039	k_pp = k + 2
1040	ENDIF
1041	IF ( k > nzt-2 ) THEN
1042	k_ppp = nzt+1
1043	ELSE
1044	k_ppp = k + 3
1045	ENDIF
1046
1047	flag_set = .FALSE.
1048	IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),0) .AND. &
1049	BTEST(wall_flags_0(k,j,i),0) .OR. &
1050	.NOT. BTEST(wall_flags_0(k_pp,j,i),0) .AND. &
1051	BTEST(wall_flags_0(k,j,i),0) .OR. &
1052	k == nzt ) &
1053	THEN
1054	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 6 )
1055	flag_set = .TRUE.
1056	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),0) .OR. &
1057	.NOT. BTEST(wall_flags_0(k_ppp,j,i),0) ) .AND. &
1058	BTEST(wall_flags_0(k-1,j,i),0) .AND. &
1059	BTEST(wall_flags_0(k,j,i),0) .AND. &
1060	BTEST(wall_flags_0(k+1,j,i),0) .AND. &
1061	BTEST(wall_flags_0(k_pp,j,i),0) .OR. &
1062	k == nzt - 1 ) &
1063	THEN
1064	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 7 )
1065	flag_set = .TRUE.
1066	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),0) &
1067	.AND. BTEST(wall_flags_0(k-1,j,i),0) &
1068	.AND. BTEST(wall_flags_0(k,j,i),0) &
1069	.AND. BTEST(wall_flags_0(k+1,j,i),0) &
1070	.AND. BTEST(wall_flags_0(k_pp,j,i),0) &
1071	.AND. BTEST(wall_flags_0(k_ppp,j,i),0) &
1072	.AND. .NOT. flag_set ) &
1073	THEN
1074	advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 8 )
1075	ENDIF
1076
1077	ENDDO
1078	ENDDO
1079	ENDDO
1080	!
1081	!-- Exchange 3D integer wall_flags.
1082	!
1083	!-- Exchange ghost points for advection flags
1084	CALL exchange_horiz_int( advc_flag, nys, nyn, nxl, nxr, nzt, nbgp )
1085	!
1086	!-- Set boundary flags at inflow and outflow boundary in case of
1087	!-- non-cyclic boundary conditions.
1088	IF ( non_cyclic_l ) THEN
1089	advc_flag(:,:,nxl-1) = advc_flag(:,:,nxl)
1090	ENDIF
1091
1092	IF ( non_cyclic_r ) THEN
1093	advc_flag(:,:,nxr+1) = advc_flag(:,:,nxr)
1094	ENDIF
1095
1096	IF ( non_cyclic_n ) THEN
1097	advc_flag(:,nyn+1,:) = advc_flag(:,nyn,:)
1098	ENDIF
1099
1100	IF ( non_cyclic_s ) THEN
1101	advc_flag(:,nys-1,:) = advc_flag(:,nys,:)
1102	ENDIF
1103
1104
1105
1106	END SUBROUTINE ws_init_flags_scalar
1107
1108	!------------------------------------------------------------------------------!
1109	! Description:
1110	! ------------
1111	!> Initialize variables used for storing statistic quantities (fluxes, variances)
1112	!------------------------------------------------------------------------------!
1113	SUBROUTINE ws_statistics
1114
1115
1116	!
1117	!-- The arrays needed for statistical evaluation are set to to 0 at the
1118	!-- beginning of prognostic_equations.
1119	IF ( ws_scheme_mom ) THEN
1120	!$ACC KERNELS PRESENT(sums_wsus_ws_l, sums_wsvs_ws_l) &
1121	!$ACC PRESENT(sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l)
1122	sums_wsus_ws_l = 0.0_wp
1123	sums_wsvs_ws_l = 0.0_wp
1124	sums_us2_ws_l = 0.0_wp
1125	sums_vs2_ws_l = 0.0_wp
1126	sums_ws2_ws_l = 0.0_wp
1127	!$ACC END KERNELS
1128	ENDIF
1129
1130	IF ( ws_scheme_sca ) THEN
1131	!$ACC KERNELS PRESENT(sums_wspts_ws_l)
1132	sums_wspts_ws_l = 0.0_wp
1133	!$ACC END KERNELS
1134	IF ( humidity ) sums_wsqs_ws_l = 0.0_wp
1135	IF ( passive_scalar ) sums_wsss_ws_l = 0.0_wp
1136
1137	ENDIF
1138
1139	END SUBROUTINE ws_statistics
1140
1141
1142	!------------------------------------------------------------------------------!
1143	! Description:
1144	! ------------
1145	!> Scalar advection - Call for grid point i,j
1146	!------------------------------------------------------------------------------!
1147	SUBROUTINE advec_s_ws_ij( advc_flag, i, j, sk, sk_char, swap_flux_y_local, &
1148	swap_diss_y_local, swap_flux_x_local, &
1149	swap_diss_x_local, i_omp, tn, &
1150	non_cyclic_l, non_cyclic_n, &
1151	non_cyclic_r, non_cyclic_s, &
1152	flux_limitation )
1153
1154
1155	CHARACTER (LEN = *), INTENT(IN) :: sk_char !< string identifier, used for assign fluxes to the correct dimension in the analysis array
1156
1157	INTEGER(iwp) :: i !< grid index along x-direction
1158	INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread
1159	INTEGER(iwp) :: j !< grid index along y-direction
1160	INTEGER(iwp) :: k !< grid index along z-direction
1161	INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults
1162	INTEGER(iwp) :: k_mmm !< k-3 index in disretization, can be modified to avoid segmentation faults
1163	INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults
1164	INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults
1165	INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
1166	INTEGER(iwp) :: tn !< number of OpenMP thread
1167
1168	INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: &
1169	advc_flag !< flag array to control order of scalar advection
1170
1171	LOGICAL :: non_cyclic_l !< flag that indicates non-cyclic boundary on the left
1172	LOGICAL :: non_cyclic_n !< flag that indicates non-cyclic boundary on the north
1173	LOGICAL :: non_cyclic_r !< flag that indicates non-cyclic boundary on the right
1174	LOGICAL :: non_cyclic_s !< flag that indicates non-cyclic boundary on the south
1175	LOGICAL, OPTIONAL :: flux_limitation !< flag indicating flux limitation of the vertical advection
1176	LOGICAL :: limiter !< control flag indicating the application of flux limitation
1177
1178	REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom
1179	REAL(wp) :: div !< velocity diverence on scalar grid
1180	REAL(wp) :: div_in !< vertical flux divergence of ingoing fluxes
1181	REAL(wp) :: div_out !< vertical flux divergence of outgoing fluxes
1182	REAL(wp) :: f_corr_t !< correction flux at grid-cell top, i.e. the difference between high and low-order flux
1183	REAL(wp) :: f_corr_d !< correction flux at grid-cell bottom, i.e. the difference between high and low-order flux
1184	REAL(wp) :: f_corr_t_in !< correction flux of ingoing flux part at grid-cell top
1185	REAL(wp) :: f_corr_d_in !< correction flux of ingoing flux part at grid-cell bottom
1186	REAL(wp) :: f_corr_t_out !< correction flux of outgoing flux part at grid-cell top
1187	REAL(wp) :: f_corr_d_out !< correction flux of outgoing flux part at grid-cell bottom
1188	REAL(wp) :: fac_correction!< factor to limit the in- and outgoing fluxes
1189	REAL(wp) :: flux_d !< 6th-order flux at grid box bottom
1190	REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction
1191	REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction
1192	REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction
1193	REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction
1194	REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction
1195	REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction
1196	REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction
1197	REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction
1198	REAL(wp) :: ibit8 !< flag indicating 5th-order scheme along z-direction
1199	REAL(wp) :: max_val !< maximum value of the quanitity along the numerical stencil (in vertical direction)
1200	REAL(wp) :: min_val !< maximum value of the quanitity along the numerical stencil (in vertical direction)
1201	REAL(wp) :: mon !< monotone solution of the advection equation using 1st-order fluxes
1202	REAL(wp) :: u_comp !< advection velocity along x-direction
1203	REAL(wp) :: v_comp !< advection velocity along y-direction
1204	!
1205	!-- sk is an array from parameter list. It should not be a pointer, because
1206	!-- in that case the compiler can not assume a stride 1 and cannot perform
1207	!-- a strided one vector load. Adding the CONTIGUOUS keyword makes things
1208	!-- even worse, because the compiler cannot assume strided one in the
1209	!-- caller side.
1210	REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !< advected scalar
1211
1212	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box
1213	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box
1214	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !< discretized artificial dissipation at top of the grid box
1215	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !< discretized 6th-order flux at northward-side of the grid box
1216	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box
1217	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6th-order flux at top of the grid box
1218	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t_1st !< discretized 1st-order flux at top of the grid box
1219
1220	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local !< discretized artificial dissipation at southward-side of the grid box
1221	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local !< discretized 6th-order flux at northward-side of the grid box
1222
1223	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_diss_x_local !< discretized artificial dissipation at leftward-side of the grid box
1224	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_flux_x_local !< discretized 6th-order flux at leftward-side of the grid box
1225	!
1226	!-- Used local modified copy of nzb_max (used to degrade order of
1227	!-- discretization) at non-cyclic boundaries. Modify only at relevant points
1228	!-- instead of the entire subdomain. This should lead to better
1229	!-- load balance between boundary and non-boundary PEs.
1230	IF( non_cyclic_l .AND. i <= nxl + 2 .OR. &
1231	non_cyclic_r .AND. i >= nxr - 2 .OR. &
1232	non_cyclic_s .AND. j <= nys + 2 .OR. &
1233	non_cyclic_n .AND. j >= nyn - 2 ) THEN
1234	nzb_max_l = nzt
1235	ELSE
1236	nzb_max_l = nzb_max
1237	END IF
1238	!
1239	!-- Set control flag for flux limiter
1240	limiter = .FALSE.
1241	IF ( PRESENT( flux_limitation) ) limiter = flux_limitation
1242	!
1243	!-- Compute southside fluxes of the respective PE bounds.
1244	IF ( j == nys ) THEN
1245	!
1246	!-- Up to the top of the highest topography.
1247	DO k = nzb+1, nzb_max_l
1248
1249	ibit5 = REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp )
1250	ibit4 = REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp )
1251	ibit3 = REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp )
1252
1253	v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k)
1254	swap_flux_y_local(k,tn) = v_comp * ( &
1255	( 37.0_wp * ibit5 * adv_sca_5 &
1256	+ 7.0_wp * ibit4 * adv_sca_3 &
1257	+ ibit3 * adv_sca_1 &
1258	) * &
1259	( sk(k,j,i) + sk(k,j-1,i) ) &
1260	- ( 8.0_wp * ibit5 * adv_sca_5 &
1261	+ ibit4 * adv_sca_3 &
1262	) * &
1263	( sk(k,j+1,i) + sk(k,j-2,i) ) &
1264	+ ( ibit5 * adv_sca_5 &
1265	) * &
1266	( sk(k,j+2,i) + sk(k,j-3,i) ) &
1267	)
1268
1269	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
1270	( 10.0_wp * ibit5 * adv_sca_5 &
1271	+ 3.0_wp * ibit4 * adv_sca_3 &
1272	+ ibit3 * adv_sca_1 &
1273	) * &
1274	( sk(k,j,i) - sk(k,j-1,i) ) &
1275	- ( 5.0_wp * ibit5 * adv_sca_5 &
1276	+ ibit4 * adv_sca_3 &
1277	) * &
1278	( sk(k,j+1,i) - sk(k,j-2,i) ) &
1279	+ ( ibit5 * adv_sca_5 &
1280	) * &
1281	( sk(k,j+2,i) - sk(k,j-3,i) ) &
1282	)
1283
1284	ENDDO
1285	!
1286	!-- Above to the top of the highest topography. No degradation necessary.
1287	DO k = nzb_max_l+1, nzt
1288
1289	v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k)
1290	swap_flux_y_local(k,tn) = v_comp * ( &
1291	37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) &
1292	- 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
1293	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) &
1294	) * adv_sca_5
1295	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
1296	10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) &
1297	- 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
1298	+ sk(k,j+2,i) - sk(k,j-3,i) &
1299	) * adv_sca_5
1300
1301	ENDDO
1302
1303	ENDIF
1304	!
1305	!-- Compute leftside fluxes of the respective PE bounds.
1306	IF ( i == i_omp ) THEN
1307
1308	DO k = nzb+1, nzb_max_l
1309
1310	ibit2 = REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp )
1311	ibit1 = REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp )
1312	ibit0 = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )
1313
1314	u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k)
1315	swap_flux_x_local(k,j,tn) = u_comp * ( &
1316	( 37.0_wp * ibit2 * adv_sca_5 &
1317	+ 7.0_wp * ibit1 * adv_sca_3 &
1318	+ ibit0 * adv_sca_1 &
1319	) * &
1320	( sk(k,j,i) + sk(k,j,i-1) ) &
1321	- ( 8.0_wp * ibit2 * adv_sca_5 &
1322	+ ibit1 * adv_sca_3 &
1323	) * &
1324	( sk(k,j,i+1) + sk(k,j,i-2) ) &
1325	+ ( ibit2 * adv_sca_5 &
1326	) * &
1327	( sk(k,j,i+2) + sk(k,j,i-3) ) &
1328	)
1329
1330	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
1331	( 10.0_wp * ibit2 * adv_sca_5 &
1332	+ 3.0_wp * ibit1 * adv_sca_3 &
1333	+ ibit0 * adv_sca_1 &
1334	) * &
1335	( sk(k,j,i) - sk(k,j,i-1) ) &
1336	- ( 5.0_wp * ibit2 * adv_sca_5 &
1337	+ ibit1 * adv_sca_3 &
1338	) * &
1339	( sk(k,j,i+1) - sk(k,j,i-2) ) &
1340	+ ( ibit2 * adv_sca_5 &
1341	) * &
1342	( sk(k,j,i+2) - sk(k,j,i-3) ) &
1343	)
1344
1345	ENDDO
1346
1347	DO k = nzb_max_l+1, nzt
1348
1349	u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k)
1350	swap_flux_x_local(k,j,tn) = u_comp * ( &
1351	37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) &
1352	- 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
1353	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) &
1354	) * adv_sca_5
1355
1356	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
1357	10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) &
1358	- 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
1359	+ ( sk(k,j,i+2) - sk(k,j,i-3) ) &
1360	) * adv_sca_5
1361
1362	ENDDO
1363
1364	ENDIF
1365	!
1366	!-- Now compute the fluxes and tendency terms for the horizontal and
1367	!-- vertical parts up to the top of the highest topography.
1368	DO k = nzb+1, nzb_max_l
1369	!
1370	!-- Note: It is faster to conduct all multiplications explicitly, e.g.
1371	!-- * adv_sca_5 ... than to determine a factor and multiplicate the
1372	!-- flux at the end.
1373
1374	ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
1375	ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp )
1376	ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp )
1377
1378	u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
1379	flux_r(k) = u_comp * ( &
1380	( 37.0_wp * ibit2 * adv_sca_5 &
1381	+ 7.0_wp * ibit1 * adv_sca_3 &
1382	+ ibit0 * adv_sca_1 &
1383	) * &
1384	( sk(k,j,i+1) + sk(k,j,i) ) &
1385	- ( 8.0_wp * ibit2 * adv_sca_5 &
1386	+ ibit1 * adv_sca_3 &
1387	) * &
1388	( sk(k,j,i+2) + sk(k,j,i-1) ) &
1389	+ ( ibit2 * adv_sca_5 &
1390	) * &
1391	( sk(k,j,i+3) + sk(k,j,i-2) ) &
1392	)
1393
1394	diss_r(k) = -ABS( u_comp ) * ( &
1395	( 10.0_wp * ibit2 * adv_sca_5 &
1396	+ 3.0_wp * ibit1 * adv_sca_3 &
1397	+ ibit0 * adv_sca_1 &
1398	) * &
1399	( sk(k,j,i+1) - sk(k,j,i) ) &
1400	- ( 5.0_wp * ibit2 * adv_sca_5 &
1401	+ ibit1 * adv_sca_3 &
1402	) * &
1403	( sk(k,j,i+2) - sk(k,j,i-1) ) &
1404	+ ( ibit2 * adv_sca_5 &
1405	) * &
1406	( sk(k,j,i+3) - sk(k,j,i-2) ) &
1407	)
1408
1409	ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp )
1410	ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp )
1411	ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp )
1412
1413	v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
1414	flux_n(k) = v_comp * ( &
1415	( 37.0_wp * ibit5 * adv_sca_5 &
1416	+ 7.0_wp * ibit4 * adv_sca_3 &
1417	+ ibit3 * adv_sca_1 &
1418	) * &
1419	( sk(k,j+1,i) + sk(k,j,i) ) &
1420	- ( 8.0_wp * ibit5 * adv_sca_5 &
1421	+ ibit4 * adv_sca_3 &
1422	) * &
1423	( sk(k,j+2,i) + sk(k,j-1,i) ) &
1424	+ ( ibit5 * adv_sca_5 &
1425	) * &
1426	( sk(k,j+3,i) + sk(k,j-2,i) ) &
1427	)
1428
1429	diss_n(k) = -ABS( v_comp ) * ( &
1430	( 10.0_wp * ibit5 * adv_sca_5 &
1431	+ 3.0_wp * ibit4 * adv_sca_3 &
1432	+ ibit3 * adv_sca_1 &
1433	) * &
1434	( sk(k,j+1,i) - sk(k,j,i) ) &
1435	- ( 5.0_wp * ibit5 * adv_sca_5 &
1436	+ ibit4 * adv_sca_3 &
1437	) * &
1438	( sk(k,j+2,i) - sk(k,j-1,i) ) &
1439	+ ( ibit5 * adv_sca_5 &
1440	) * &
1441	( sk(k,j+3,i) - sk(k,j-2,i) ) &
1442	)
1443	ENDDO
1444	!
1445	!-- Now compute the fluxes and tendency terms for the horizontal and
1446	!-- vertical parts above the top of the highest topography. No degradation
1447	!-- for the horizontal parts, but for the vertical it is stell needed.
1448	DO k = nzb_max_l+1, nzt
1449
1450	u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
1451	flux_r(k) = u_comp * ( &
1452	37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) &
1453	- 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
1454	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
1455	diss_r(k) = -ABS( u_comp ) * ( &
1456	10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) &
1457	- 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
1458	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
1459
1460	v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
1461	flux_n(k) = v_comp * ( &
1462	37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) &
1463	- 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
1464	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
1465	diss_n(k) = -ABS( v_comp ) * ( &
1466	10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) &
1467	- 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
1468	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
1469
1470	ENDDO
1471	!
1472	!-- Now, compute vertical fluxes. Split loop into a part treating the
1473	!-- lowest 2 grid points with indirect indexing, a main loop without
1474	!-- indirect indexing, and a loop for the uppermost 2 grip points with
1475	!-- indirect indexing. This allows better vectorization for the main loop.
1476	!-- First, compute the flux at model surface, which need has to be
1477	!-- calculated explicetely for the tendency at
1478	!-- the first w-level. For topography wall this is done implicitely by
1479	!-- advc_flag.
1480	flux_t(nzb) = 0.0_wp
1481	diss_t(nzb) = 0.0_wp
1482	DO k = nzb+1, nzb+2
1483	ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
1484	ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
1485	ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
1486	!
1487	!-- k index has to be modified near bottom and top, else array
1488	!-- subscripts will be exceeded.
1489	k_ppp = k + 3 * ibit8
1490	k_pp = k + 2 * ( 1 - ibit6 )
1491	k_mm = k - 2 * ibit8
1492
1493
1494	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
1495	( 37.0_wp * ibit8 * adv_sca_5 &
1496	+ 7.0_wp * ibit7 * adv_sca_3 &
1497	+ ibit6 * adv_sca_1 &
1498	) * &
1499	( sk(k+1,j,i) + sk(k,j,i) ) &
1500	- ( 8.0_wp * ibit8 * adv_sca_5 &
1501	+ ibit7 * adv_sca_3 &
1502	) * &
1503	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
1504	+ ( ibit8 * adv_sca_5 &
1505	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
1506	)
1507
1508	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
1509	( 10.0_wp * ibit8 * adv_sca_5 &
1510	+ 3.0_wp * ibit7 * adv_sca_3 &
1511	+ ibit6 * adv_sca_1 &
1512	) * &
1513	( sk(k+1,j,i) - sk(k,j,i) ) &
1514	- ( 5.0_wp * ibit8 * adv_sca_5 &
1515	+ ibit7 * adv_sca_3 &
1516	) * &
1517	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
1518	+ ( ibit8 * adv_sca_5 &
1519	) * &
1520	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
1521	)
1522	ENDDO
1523
1524	DO k = nzb+3, nzt-2
1525	ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
1526	ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
1527	ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
1528
1529	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
1530	( 37.0_wp * ibit8 * adv_sca_5 &
1531	+ 7.0_wp * ibit7 * adv_sca_3 &
1532	+ ibit6 * adv_sca_1 &
1533	) * &
1534	( sk(k+1,j,i) + sk(k,j,i) ) &
1535	- ( 8.0_wp * ibit8 * adv_sca_5 &
1536	+ ibit7 * adv_sca_3 &
1537	) * &
1538	( sk(k+2,j,i) + sk(k-1,j,i) ) &
1539	+ ( ibit8 * adv_sca_5 &
1540	) * ( sk(k+3,j,i)+ sk(k-2,j,i) ) &
1541	)
1542
1543	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
1544	( 10.0_wp * ibit8 * adv_sca_5 &
1545	+ 3.0_wp * ibit7 * adv_sca_3 &
1546	+ ibit6 * adv_sca_1 &
1547	) * &
1548	( sk(k+1,j,i) - sk(k,j,i) ) &
1549	- ( 5.0_wp * ibit8 * adv_sca_5 &
1550	+ ibit7 * adv_sca_3 &
1551	) * &
1552	( sk(k+2,j,i) - sk(k-1,j,i) ) &
1553	+ ( ibit8 * adv_sca_5 &
1554	) * &
1555	( sk(k+3,j,i) - sk(k-2,j,i) ) &
1556	)
1557	ENDDO
1558
1559	DO k = nzt-1, nzt
1560	ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
1561	ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
1562	ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
1563	!
1564	!-- k index has to be modified near bottom and top, else array
1565	!-- subscripts will be exceeded.
1566	k_ppp = k + 3 * ibit8
1567	k_pp = k + 2 * ( 1 - ibit6 )
1568	k_mm = k - 2 * ibit8
1569
1570
1571	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
1572	( 37.0_wp * ibit8 * adv_sca_5 &
1573	+ 7.0_wp * ibit7 * adv_sca_3 &
1574	+ ibit6 * adv_sca_1 &
1575	) * &
1576	( sk(k+1,j,i) + sk(k,j,i) ) &
1577	- ( 8.0_wp * ibit8 * adv_sca_5 &
1578	+ ibit7 * adv_sca_3 &
1579	) * &
1580	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
1581	+ ( ibit8 * adv_sca_5 &
1582	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
1583	)
1584
1585	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
1586	( 10.0_wp * ibit8 * adv_sca_5 &
1587	+ 3.0_wp * ibit7 * adv_sca_3 &
1588	+ ibit6 * adv_sca_1 &
1589	) * &
1590	( sk(k+1,j,i) - sk(k,j,i) ) &
1591	- ( 5.0_wp * ibit8 * adv_sca_5 &
1592	+ ibit7 * adv_sca_3 &
1593	) * &
1594	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
1595	+ ( ibit8 * adv_sca_5 &
1596	) * &
1597	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
1598	)
1599	ENDDO
1600
1601	IF ( limiter ) THEN
1602	!
1603	!-- Compute monotone first-order fluxes which are required for mononte
1604	!-- flux limitation.
1605	flux_t_1st(nzb) = 0.0_wp
1606	DO k = nzb+1, nzb_max_l
1607	flux_t_1st(k) = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) &
1608	- ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) ) &
1609	* rho_air_zw(k) * adv_sca_1
1610	!
1611	!-- In flux limitation the total flux will be corrected. For the sake
1612	!-- of cleariness the higher-order advective and disspative fluxes
1613	!-- will be merged onto flux_t.
1614	flux_t(k) = flux_t(k) + diss_t(k)
1615	diss_t(k) = 0.0_wp
1616	ENDDO
1617	!
1618	!-- Flux limitation of vertical fluxes according to Skamarock (2006).
1619	!-- Please note, as flux limitation implies linear dependencies of fluxes,
1620	!-- flux limitation is only made for the vertical advection term. Limitation
1621	!-- of the horizontal terms cannot be parallelized.
1622	!-- Due to the linear dependency, the following loop will not be vectorized.
1623	!-- Further, note that the flux limiter is only applied within the urban
1624	!-- layer, i.e up to the topography top.
1625	DO k = nzb+1, nzb_max_l
1626	!
1627	!-- Compute one-dimensional divergence along the vertical direction,
1628	!-- which is used to correct the advection discretization. This is
1629	!-- necessary as in one-dimensional space the advection velocity
1630	!-- should actually be constant.
1631	div = ( w(k,j,i) * rho_air_zw(k) &
1632	- w(k-1,j,i) * rho_air_zw(k-1) &
1633	) * drho_air(k) * ddzw(k)
1634	!
1635	!-- Compute monotone solution of the advection equation from
1636	!-- 1st-order fluxes. Please note, the advection equation is corrected
1637	!-- by the divergence term (in 1D the advective flow should be divergence
1638	!-- free). Moreover, please note, as time-increment the full timestep
1639	!-- is used, even though a Runge-Kutta scheme will be used. However,
1640	!-- the length of the actual time increment is not important at all
1641	!-- since it cancels out later when the fluxes are limited.
1642	mon = sk(k,j,i) + ( - ( flux_t_1st(k) - flux_t_1st(k-1) ) &
1643	* drho_air(k) * ddzw(k) &
1644	+ div * sk(k,j,i) &
1645	) * dt_3d
1646	!
1647	!-- Determine minimum and maximum values along the numerical stencil.
1648	k_mmm = MAX( k - 3, nzb + 1 )
1649	k_ppp = MIN( k + 3, nzt + 1 )
1650
1651	min_val = MINVAL( sk(k_mmm:k_ppp,j,i) )
1652	max_val = MAXVAL( sk(k_mmm:k_ppp,j,i) )
1653	!
1654	!-- Compute difference between high- and low-order fluxes, which may
1655	!-- act as correction fluxes
1656	f_corr_t = flux_t(k) - flux_t_1st(k)
1657	f_corr_d = flux_t(k-1) - flux_t_1st(k-1)
1658	!
1659	!-- Determine outgoing fluxes, i.e. the part of the fluxes which can
1660	!-- decrease the value within the grid box
1661	f_corr_t_out = MAX( 0.0_wp, f_corr_t )
1662	f_corr_d_out = MIN( 0.0_wp, f_corr_d )
1663	!
1664	!-- Determine ingoing fluxes, i.e. the part of the fluxes which can
1665	!-- increase the value within the grid box
1666	f_corr_t_in = MIN( 0.0_wp, f_corr_t)
1667	f_corr_d_in = MAX( 0.0_wp, f_corr_d)
1668	!
1669	!-- Compute divergence of outgoing correction fluxes
1670	div_out = - ( f_corr_t_out - f_corr_d_out ) * drho_air(k) &
1671	* ddzw(k) * dt_3d
1672	!
1673	!-- Compute divergence of ingoing correction fluxes
1674	div_in = - ( f_corr_t_in - f_corr_d_in ) * drho_air(k) &
1675	* ddzw(k) * dt_3d
1676	!
1677	!-- Check if outgoing fluxes can lead to undershoots, i.e. values smaller
1678	!-- than the minimum value within the numerical stencil. If so, limit
1679	!-- them.
1680	IF ( mon - min_val < - div_out .AND. ABS( div_out ) > 0.0_wp ) &
1681	THEN
1682	fac_correction = ( mon - min_val ) / ( - div_out )
1683	f_corr_t_out = f_corr_t_out * fac_correction
1684	f_corr_d_out = f_corr_d_out * fac_correction
1685	ENDIF
1686	!
1687	!-- Check if ingoing fluxes can lead to overshoots, i.e. values larger
1688	!-- than the maximum value within the numerical stencil. If so, limit
1689	!-- them.
1690	IF ( mon - max_val > - div_in .AND. ABS( div_in ) > 0.0_wp ) &
1691	THEN
1692	fac_correction = ( mon - max_val ) / ( - div_in )
1693	f_corr_t_in = f_corr_t_in * fac_correction
1694	f_corr_d_in = f_corr_d_in * fac_correction
1695	ENDIF
1696	!
1697	!-- Finally add the limited fluxes to the original ones. If no
1698	!-- flux limitation was done, the fluxes equal the original ones.
1699	flux_t(k) = flux_t_1st(k) + f_corr_t_out + f_corr_t_in
1700	flux_t(k-1) = flux_t_1st(k-1) + f_corr_d_out + f_corr_d_in
1701	ENDDO
1702	ENDIF
1703
1704	DO k = nzb+1, nzb_max_l
1705
1706	flux_d = flux_t(k-1)
1707	diss_d = diss_t(k-1)
1708
1709	ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
1710	ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp )
1711	ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp )
1712
1713	ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp )
1714	ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp )
1715	ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp )
1716
1717	ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
1718	ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
1719	ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
1720	!
1721	!-- Calculate the divergence of the velocity field. A respective
1722	!-- correction is needed to overcome numerical instabilities introduced
1723	!-- by a not sufficient reduction of divergences near topography.
1724	div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) &
1725	- u(k,j,i) * ( &
1726	REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp ) &
1727	+ REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp ) &
1728	+ REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp ) &
1729	) &
1730	) * ddx &
1731	+ ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) &
1732	- v(k,j,i) * ( &
1733	REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp ) &
1734	+ REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp ) &
1735	+ REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp ) &
1736	) &
1737	) * ddy &
1738	+ ( w(k,j,i) * rho_air_zw(k) * &
1739	( ibit6 + ibit7 + ibit8 ) &
1740	- w(k-1,j,i) * rho_air_zw(k-1) * &
1741	( &
1742	REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp ) &
1743	+ REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp ) &
1744	+ REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp ) &
1745	) &
1746	) * drho_air(k) * ddzw(k)
1747
1748	tend(k,j,i) = tend(k,j,i) - ( &
1749	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
1750	swap_diss_x_local(k,j,tn) ) * ddx &
1751	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
1752	swap_diss_y_local(k,tn) ) * ddy &
1753	+ ( ( flux_t(k) + diss_t(k) ) - &
1754	( flux_d + diss_d ) &
1755	) * drho_air(k) * ddzw(k) &
1756	) + sk(k,j,i) * div
1757
1758
1759	swap_flux_y_local(k,tn) = flux_n(k)
1760	swap_diss_y_local(k,tn) = diss_n(k)
1761	swap_flux_x_local(k,j,tn) = flux_r(k)
1762	swap_diss_x_local(k,j,tn) = diss_r(k)
1763
1764	ENDDO
1765
1766	DO k = nzb_max_l+1, nzt
1767
1768	flux_d = flux_t(k-1)
1769	diss_d = diss_t(k-1)
1770	!
1771	!-- Calculate the divergence of the velocity field. A respective
1772	!-- correction is needed to overcome numerical instabilities introduced
1773	!-- by a not sufficient reduction of divergences near topography.
1774	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
1775	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
1776	+ ( w(k,j,i) * rho_air_zw(k) &
1777	- w(k-1,j,i) * rho_air_zw(k-1) &
1778	) * drho_air(k) * ddzw(k)
1779
1780	tend(k,j,i) = tend(k,j,i) - ( &
1781	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
1782	swap_diss_x_local(k,j,tn) ) * ddx &
1783	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
1784	swap_diss_y_local(k,tn) ) * ddy &
1785	+ ( ( flux_t(k) + diss_t(k) ) - &
1786	( flux_d + diss_d ) &
1787	) * drho_air(k) * ddzw(k) &
1788	) + sk(k,j,i) * div
1789
1790
1791	swap_flux_y_local(k,tn) = flux_n(k)
1792	swap_diss_y_local(k,tn) = diss_n(k)
1793	swap_flux_x_local(k,j,tn) = flux_r(k)
1794	swap_diss_x_local(k,j,tn) = diss_r(k)
1795
1796	ENDDO
1797
1798	!
1799	!-- Evaluation of statistics.
1800	SELECT CASE ( sk_char )
1801
1802	CASE ( 'pt' )
1803
1804	DO k = nzb, nzt
1805	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + &
1806	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1807	* ( w(k,j,i) - hom(k,1,3,0) ) &
1808	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1809	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1810	) * weight_substep(intermediate_timestep_count)
1811	ENDDO
1812
1813	CASE ( 'sa' )
1814
1815	DO k = nzb, nzt
1816	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) + &
1817	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1818	* ( w(k,j,i) - hom(k,1,3,0) ) &
1819	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1820	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1821	) * weight_substep(intermediate_timestep_count)
1822	ENDDO
1823
1824	CASE ( 'q' )
1825
1826	DO k = nzb, nzt
1827	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) + &
1828	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1829	* ( w(k,j,i) - hom(k,1,3,0) ) &
1830	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1831	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1832	) * weight_substep(intermediate_timestep_count)
1833	ENDDO
1834
1835	CASE ( 'qc' )
1836
1837	DO k = nzb, nzt
1838	sums_wsqcs_ws_l(k,tn) = sums_wsqcs_ws_l(k,tn) + &
1839	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1840	* ( w(k,j,i) - hom(k,1,3,0) ) &
1841	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1842	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1843	) * weight_substep(intermediate_timestep_count)
1844	ENDDO
1845
1846
1847	CASE ( 'qr' )
1848
1849	DO k = nzb, nzt
1850	sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) + &
1851	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1852	* ( w(k,j,i) - hom(k,1,3,0) ) &
1853	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1854	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1855	) * weight_substep(intermediate_timestep_count)
1856	ENDDO
1857
1858	CASE ( 'nc' )
1859
1860	DO k = nzb, nzt
1861	sums_wsncs_ws_l(k,tn) = sums_wsncs_ws_l(k,tn) + &
1862	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1863	* ( w(k,j,i) - hom(k,1,3,0) ) &
1864	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1865	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1866	) * weight_substep(intermediate_timestep_count)
1867	ENDDO
1868
1869	CASE ( 'nr' )
1870
1871	DO k = nzb, nzt
1872	sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) + &
1873	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1874	* ( w(k,j,i) - hom(k,1,3,0) ) &
1875	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1876	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1877	) * weight_substep(intermediate_timestep_count)
1878	ENDDO
1879
1880	CASE ( 's' )
1881
1882	DO k = nzb, nzt
1883	sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn) + &
1884	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1885	* ( w(k,j,i) - hom(k,1,3,0) ) &
1886	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1887	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1888	) * weight_substep(intermediate_timestep_count)
1889	ENDDO
1890
1891	CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' )
1892
1893	DO k = nzb, nzt
1894	sums_salsa_ws_l(k,tn) = sums_salsa_ws_l(k,tn) + &
1895	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1896	* ( w(k,j,i) - hom(k,1,3,0) ) &
1897	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1898	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1899	) * weight_substep(intermediate_timestep_count)
1900	ENDDO
1901
1902	! CASE ( 'kc' )
1903	!kk Has to be implemented for kpp chemistry
1904
1905
1906	END SELECT
1907
1908	END SUBROUTINE advec_s_ws_ij
1909
1910
1911
1912
1913	!------------------------------------------------------------------------------!
1914	! Description:
1915	! ------------
1916	!> Advection of u-component - Call for grid point i,j
1917	!------------------------------------------------------------------------------!
1918	SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
1919
1920
1921	INTEGER(iwp) :: i !< grid index along x-direction
1922	INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread
1923	INTEGER(iwp) :: j !< grid index along y-direction
1924	INTEGER(iwp) :: k !< grid index along z-direction
1925	INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults
1926	INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults
1927	INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults
1928	INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
1929	INTEGER(iwp) :: tn !< number of OpenMP thread
1930
1931	REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction
1932	REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction
1933	REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction
1934	REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction
1935	REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction
1936	REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction
1937	REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction
1938	REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction
1939	REAL(wp) :: ibit8 !< flag indicating 5th-order scheme along z-direction
1940	REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom
1941	REAL(wp) :: div !< diverence on u-grid
1942	REAL(wp) :: flux_d !< 6th-order flux at grid box bottom
1943	REAL(wp) :: gu !< Galilei-transformation velocity along x
1944	REAL(wp) :: gv !< Galilei-transformation velocity along y
1945	REAL(wp) :: u_comp_l !< advection velocity along x at leftmost grid point on subdomain
1946
1947	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box
1948	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box
1949	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !< discretized artificial dissipation at top of the grid box
1950	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !< discretized 6th-order flux at northward-side of the grid box
1951	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box
1952	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6th-order flux at top of the grid box
1953	REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !< advection velocity along x
1954	REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !< advection velocity along y
1955	REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< advection velocity along z
1956	!
1957	!-- Used local modified copy of nzb_max (used to degrade order of
1958	!-- discretization) at non-cyclic boundaries. Modify only at relevant points
1959	!-- instead of the entire subdomain. This should lead to better
1960	!-- load balance between boundary and non-boundary PEs.
1961	IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR. &
1962	( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i >= nxr - 2 .OR. &
1963	( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j <= nys + 2 .OR. &
1964	( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j >= nyn - 2 ) THEN
1965	nzb_max_l = nzt
1966	ELSE
1967	nzb_max_l = nzb_max
1968	END IF
1969
1970	gu = 2.0_wp * u_gtrans
1971	gv = 2.0_wp * v_gtrans
1972	!
1973	!-- Compute southside fluxes for the respective boundary of PE
1974	IF ( j == nys ) THEN
1975
1976	DO k = nzb+1, nzb_max_l
1977
1978	ibit5 = REAL( IBITS(advc_flags_m(k,j-1,i),5,1), KIND = wp )
1979	ibit4 = REAL( IBITS(advc_flags_m(k,j-1,i),4,1), KIND = wp )
1980	ibit3 = REAL( IBITS(advc_flags_m(k,j-1,i),3,1), KIND = wp )
1981
1982	v_comp(k) = v(k,j,i) + v(k,j,i-1) - gv
1983	flux_s_u(k,tn) = v_comp(k) * ( &
1984	( 37.0_wp * ibit5 * adv_mom_5 &
1985	+ 7.0_wp * ibit4 * adv_mom_3 &
1986	+ ibit3 * adv_mom_1 &
1987	) * &
1988	( u(k,j,i) + u(k,j-1,i) ) &
1989	- ( 8.0_wp * ibit5 * adv_mom_5 &
1990	+ ibit4 * adv_mom_3 &
1991	) * &
1992	( u(k,j+1,i) + u(k,j-2,i) ) &
1993	+ ( ibit5 * adv_mom_5 &
1994	) * &
1995	( u(k,j+2,i) + u(k,j-3,i) ) &
1996	)
1997
1998	diss_s_u(k,tn) = - ABS ( v_comp(k) ) * ( &
1999	( 10.0_wp * ibit5 * adv_mom_5 &
2000	+ 3.0_wp * ibit4 * adv_mom_3 &
2001	+ ibit3 * adv_mom_1 &
2002	) * &
2003	( u(k,j,i) - u(k,j-1,i) ) &
2004	- ( 5.0_wp * ibit5 * adv_mom_5 &
2005	+ ibit4 * adv_mom_3 &
2006	) * &
2007	( u(k,j+1,i) - u(k,j-2,i) ) &
2008	+ ( ibit5 * adv_mom_5 &
2009	) * &
2010	( u(k,j+2,i) - u(k,j-3,i) ) &
2011	)
2012
2013	ENDDO
2014
2015	DO k = nzb_max_l+1, nzt
2016
2017	v_comp(k) = v(k,j,i) + v(k,j,i-1) - gv
2018	flux_s_u(k,tn) = v_comp(k) * ( &
2019	37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) &
2020	- 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) &
2021	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
2022	diss_s_u(k,tn) = - ABS(v_comp(k)) * ( &
2023	10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) &
2024	- 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) &
2025	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
2026
2027	ENDDO
2028
2029	ENDIF
2030	!
2031	!-- Compute leftside fluxes for the respective boundary of PE
2032	IF ( i == i_omp .OR. i == nxlu ) THEN
2033
2034	DO k = nzb+1, nzb_max_l
2035
2036	ibit2 = REAL( IBITS(advc_flags_m(k,j,i-1),2,1), KIND = wp )
2037	ibit1 = REAL( IBITS(advc_flags_m(k,j,i-1),1,1), KIND = wp )
2038	ibit0 = REAL( IBITS(advc_flags_m(k,j,i-1),0,1), KIND = wp )
2039
2040	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
2041	flux_l_u(k,j,tn) = u_comp_l * ( &
2042	( 37.0_wp * ibit2 * adv_mom_5 &
2043	+ 7.0_wp * ibit1 * adv_mom_3 &
2044	+ ibit0 * adv_mom_1 &
2045	) * &
2046	( u(k,j,i) + u(k,j,i-1) ) &
2047	- ( 8.0_wp * ibit2 * adv_mom_5 &
2048	+ ibit1 * adv_mom_3 &
2049	) * &
2050	( u(k,j,i+1) + u(k,j,i-2) ) &
2051	+ ( ibit2 * adv_mom_5 &
2052	) * &
2053	( u(k,j,i+2) + u(k,j,i-3) ) &
2054	)
2055
2056	diss_l_u(k,j,tn) = - ABS( u_comp_l ) * ( &
2057	( 10.0_wp * ibit2 * adv_mom_5 &
2058	+ 3.0_wp * ibit1 * adv_mom_3 &
2059	+ ibit0 * adv_mom_1 &
2060	) * &
2061	( u(k,j,i) - u(k,j,i-1) ) &
2062	- ( 5.0_wp * ibit2 * adv_mom_5 &
2063	+ ibit1 * adv_mom_3 &
2064	) * &
2065	( u(k,j,i+1) - u(k,j,i-2) ) &
2066	+ ( ibit2 * adv_mom_5 &
2067	) * &
2068	( u(k,j,i+2) - u(k,j,i-3) ) &
2069	)
2070
2071	ENDDO
2072
2073	DO k = nzb_max_l+1, nzt
2074
2075	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
2076	flux_l_u(k,j,tn) = u_comp_l * ( &
2077	37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) &
2078	- 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) &
2079	+ ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
2080	diss_l_u(k,j,tn) = - ABS(u_comp_l) * ( &
2081	10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) &
2082	- 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) &
2083	+ ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
2084
2085	ENDDO
2086
2087	ENDIF
2088	!
2089	!-- Now compute the fluxes tendency terms for the horizontal and
2090	!-- vertical parts.
2091	DO k = nzb+1, nzb_max_l
2092
2093	ibit2 = REAL( IBITS(advc_flags_m(k,j,i),2,1), KIND = wp )
2094	ibit1 = REAL( IBITS(advc_flags_m(k,j,i),1,1), KIND = wp )
2095	ibit0 = REAL( IBITS(advc_flags_m(k,j,i),0,1), KIND = wp )
2096
2097	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2098	flux_r(k) = ( u_comp(k) - gu ) * ( &
2099	( 37.0_wp * ibit2 * adv_mom_5 &
2100	+ 7.0_wp * ibit1 * adv_mom_3 &
2101	+ ibit0 * adv_mom_1 &
2102	) * &
2103	( u(k,j,i+1) + u(k,j,i) ) &
2104	- ( 8.0_wp * ibit2 * adv_mom_5 &
2105	+ ibit1 * adv_mom_3 &
2106	) * &
2107	( u(k,j,i+2) + u(k,j,i-1) ) &
2108	+ ( ibit2 * adv_mom_5 &
2109	) * &
2110	( u(k,j,i+3) + u(k,j,i-2) ) &
2111	)
2112
2113	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
2114	( 10.0_wp * ibit2 * adv_mom_5 &
2115	+ 3.0_wp * ibit1 * adv_mom_3 &
2116	+ ibit0 * adv_mom_1 &
2117	) * &
2118	( u(k,j,i+1) - u(k,j,i) ) &
2119	- ( 5.0_wp * ibit2 * adv_mom_5 &
2120	+ ibit1 * adv_mom_3 &
2121	) * &
2122	( u(k,j,i+2) - u(k,j,i-1) ) &
2123	+ ( ibit2 * adv_mom_5 &
2124	) * &
2125	( u(k,j,i+3) - u(k,j,i-2) ) &
2126	)
2127
2128	ibit5 = REAL( IBITS(advc_flags_m(k,j,i),5,1), KIND = wp )
2129	ibit4 = REAL( IBITS(advc_flags_m(k,j,i),4,1), KIND = wp )
2130	ibit3 = REAL( IBITS(advc_flags_m(k,j,i),3,1), KIND = wp )
2131
2132	v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
2133	flux_n(k) = v_comp(k) * ( &
2134	( 37.0_wp * ibit5 * adv_mom_5 &
2135	+ 7.0_wp * ibit4 * adv_mom_3 &
2136	+ ibit3 * adv_mom_1 &
2137	) * &
2138	( u(k,j+1,i) + u(k,j,i) ) &
2139	- ( 8.0_wp * ibit5 * adv_mom_5 &
2140	+ ibit4 * adv_mom_3 &
2141	) * &
2142	( u(k,j+2,i) + u(k,j-1,i) ) &
2143	+ ( ibit5 * adv_mom_5 &
2144	) * &
2145	( u(k,j+3,i) + u(k,j-2,i) ) &
2146	)
2147
2148	diss_n(k) = - ABS ( v_comp(k) ) * ( &
2149	( 10.0_wp * ibit5 * adv_mom_5 &
2150	+ 3.0_wp * ibit4 * adv_mom_3 &
2151	+ ibit3 * adv_mom_1 &
2152	) * &
2153	( u(k,j+1,i) - u(k,j,i) ) &
2154	- ( 5.0_wp * ibit5 * adv_mom_5 &
2155	+ ibit4 * adv_mom_3 &
2156	) * &
2157	( u(k,j+2,i) - u(k,j-1,i) ) &
2158	+ ( ibit5 * adv_mom_5 &
2159	) * &
2160	( u(k,j+3,i) - u(k,j-2,i) ) &
2161	)
2162	ENDDO
2163
2164	DO k = nzb_max_l+1, nzt
2165
2166	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2167	flux_r(k) = ( u_comp(k) - gu ) * ( &
2168	37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) &
2169	- 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) &
2170	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
2171	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
2172	10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) &
2173	- 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) &
2174	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
2175
2176	v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
2177	flux_n(k) = v_comp(k) * ( &
2178	37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) &
2179	- 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) &
2180	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
2181	diss_n(k) = - ABS( v_comp(k) ) * ( &
2182	10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) &
2183	- 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) &
2184	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
2185
2186	ENDDO
2187	!
2188	!-- Now, compute vertical fluxes. Split loop into a part treating the
2189	!-- lowest 2 grid points with indirect indexing, a main loop without
2190	!-- indirect indexing, and a loop for the uppermost 2 grip points with
2191	!-- indirect indexing. This allows better vectorization for the main loop.
2192	!-- First, compute the flux at model surface, which need has to be
2193	!-- calculated explicetely for the tendency at
2194	!-- the first w-level. For topography wall this is done implicitely by
2195	!-- advc_flags_m.
2196	flux_t(nzb) = 0.0_wp
2197	diss_t(nzb) = 0.0_wp
2198	w_comp(nzb) = 0.0_wp
2199	DO k = nzb+1, nzb+2
2200	!
2201	!-- k index has to be modified near bottom and top, else array
2202	!-- subscripts will be exceeded.
2203	ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
2204	ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
2205	ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp )
2206
2207	k_ppp = k + 3 * ibit8
2208	k_pp = k + 2 * ( 1 - ibit6 )
2209	k_mm = k - 2 * ibit8
2210
2211	w_comp(k) = w(k,j,i) + w(k,j,i-1)
2212	flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &
2213	( 37.0_wp * ibit8 * adv_mom_5 &
2214	+ 7.0_wp * ibit7 * adv_mom_3 &
2215	+ ibit6 * adv_mom_1 &
2216	) * &
2217	( u(k+1,j,i) + u(k,j,i) ) &
2218	- ( 8.0_wp * ibit8 * adv_mom_5 &
2219	+ ibit7 * adv_mom_3 &
2220	) * &
2221	( u(k_pp,j,i) + u(k-1,j,i) ) &
2222	+ ( ibit8 * adv_mom_5 &
2223	) * &
2224	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
2225	)
2226
2227	diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &
2228	( 10.0_wp * ibit8 * adv_mom_5 &
2229	+ 3.0_wp * ibit7 * adv_mom_3 &
2230	+ ibit6 * adv_mom_1 &
2231	) * &
2232	( u(k+1,j,i) - u(k,j,i) ) &
2233	- ( 5.0_wp * ibit8 * adv_mom_5 &
2234	+ ibit7 * adv_mom_3 &
2235	) * &
2236	( u(k_pp,j,i) - u(k-1,j,i) ) &
2237	+ ( ibit8 * adv_mom_5 &
2238	) * &
2239	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
2240	)
2241	ENDDO
2242
2243	DO k = nzb+3, nzt-2
2244
2245	ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
2246	ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
2247	ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp )
2248
2249	w_comp(k) = w(k,j,i) + w(k,j,i-1)
2250	flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &
2251	( 37.0_wp * ibit8 * adv_mom_5 &
2252	+ 7.0_wp * ibit7 * adv_mom_3 &
2253	+ ibit6 * adv_mom_1 &
2254	) * &
2255	( u(k+1,j,i) + u(k,j,i) ) &
2256	- ( 8.0_wp * ibit8 * adv_mom_5 &
2257	+ ibit7 * adv_mom_3 &
2258	) * &
2259	( u(k+2,j,i) + u(k-1,j,i) ) &
2260	+ ( ibit8 * adv_mom_5 &
2261	) * &
2262	( u(k+3,j,i) + u(k-2,j,i) ) &
2263	)
2264
2265	diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &
2266	( 10.0_wp * ibit8 * adv_mom_5 &
2267	+ 3.0_wp * ibit7 * adv_mom_3 &
2268	+ ibit6 * adv_mom_1 &
2269	) * &
2270	( u(k+1,j,i) - u(k,j,i) ) &
2271	- ( 5.0_wp * ibit8 * adv_mom_5 &
2272	+ ibit7 * adv_mom_3 &
2273	) * &
2274	( u(k+2,j,i) - u(k-1,j,i) ) &
2275	+ ( ibit8 * adv_mom_5 &
2276	) * &
2277	( u(k+3,j,i) - u(k-2,j,i) ) &
2278	)
2279	ENDDO
2280
2281	DO k = nzt-1, nzt
2282	!
2283	!-- k index has to be modified near bottom and top, else array
2284	!-- subscripts will be exceeded.
2285	ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
2286	ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
2287	ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp )
2288
2289	k_ppp = k + 3 * ibit8
2290	k_pp = k + 2 * ( 1 - ibit6 )
2291	k_mm = k - 2 * ibit8
2292
2293	w_comp(k) = w(k,j,i) + w(k,j,i-1)
2294	flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &
2295	( 37.0_wp * ibit8 * adv_mom_5 &
2296	+ 7.0_wp * ibit7 * adv_mom_3 &
2297	+ ibit6 * adv_mom_1 &
2298	) * &
2299	( u(k+1,j,i) + u(k,j,i) ) &
2300	- ( 8.0_wp * ibit8 * adv_mom_5 &
2301	+ ibit7 * adv_mom_3 &
2302	) * &
2303	( u(k_pp,j,i) + u(k-1,j,i) ) &
2304	+ ( ibit8 * adv_mom_5 &
2305	) * &
2306	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
2307	)
2308
2309	diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &
2310	( 10.0_wp * ibit8 * adv_mom_5 &
2311	+ 3.0_wp * ibit7 * adv_mom_3 &
2312	+ ibit6 * adv_mom_1 &
2313	) * &
2314	( u(k+1,j,i) - u(k,j,i) ) &
2315	- ( 5.0_wp * ibit8 * adv_mom_5 &
2316	+ ibit7 * adv_mom_3 &
2317	) * &
2318	( u(k_pp,j,i) - u(k-1,j,i) ) &
2319	+ ( ibit8 * adv_mom_5 &
2320	) * &
2321	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
2322	)
2323	ENDDO
2324
2325	DO k = nzb+1, nzb_max_l
2326
2327	flux_d = flux_t(k-1)
2328	diss_d = diss_t(k-1)
2329
2330	ibit2 = REAL( IBITS(advc_flags_m(k,j,i),2,1), KIND = wp )
2331	ibit1 = REAL( IBITS(advc_flags_m(k,j,i),1,1), KIND = wp )
2332	ibit0 = REAL( IBITS(advc_flags_m(k,j,i),0,1), KIND = wp )
2333
2334	ibit5 = REAL( IBITS(advc_flags_m(k,j,i),5,1), KIND = wp )
2335	ibit4 = REAL( IBITS(advc_flags_m(k,j,i),4,1), KIND = wp )
2336	ibit3 = REAL( IBITS(advc_flags_m(k,j,i),3,1), KIND = wp )
2337
2338	ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
2339	ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
2340	ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp )
2341	!
2342	!-- Calculate the divergence of the velocity field. A respective
2343	!-- correction is needed to overcome numerical instabilities introduced
2344	!-- by a not sufficient reduction of divergences near topography.
2345	div = ( ( u_comp(k) * ( ibit0 + ibit1 + ibit2 ) &
2346	- ( u(k,j,i) + u(k,j,i-1) ) &
2347	* ( &
2348	REAL( IBITS(advc_flags_m(k,j,i-1),0,1), KIND = wp ) &
2349	+ REAL( IBITS(advc_flags_m(k,j,i-1),1,1), KIND = wp ) &
2350	+ REAL( IBITS(advc_flags_m(k,j,i-1),2,1), KIND = wp ) &
2351	) &
2352	) * ddx &
2353	+ ( ( v_comp(k) + gv ) * ( ibit3 + ibit4 + ibit5 ) &
2354	- ( v(k,j,i) + v(k,j,i-1 ) ) &
2355	* ( &
2356	REAL( IBITS(advc_flags_m(k,j-1,i),3,1), KIND = wp ) &
2357	+ REAL( IBITS(advc_flags_m(k,j-1,i),4,1), KIND = wp ) &
2358	+ REAL( IBITS(advc_flags_m(k,j-1,i),5,1), KIND = wp ) &
2359	) &
2360	) * ddy &
2361	+ ( w_comp(k) * rho_air_zw(k) * ( ibit6 + ibit7 + ibit8 ) &
2362	- w_comp(k-1) * rho_air_zw(k-1) &
2363	* ( &
2364	REAL( IBITS(advc_flags_m(k-1,j,i),6,1), KIND = wp ) &
2365	+ REAL( IBITS(advc_flags_m(k-1,j,i),7,1), KIND = wp ) &
2366	+ REAL( IBITS(advc_flags_m(k-1,j,i),8,1), KIND = wp ) &
2367	) &
2368	) * drho_air(k) * ddzw(k) &
2369	) * 0.5_wp
2370
2371	tend(k,j,i) = tend(k,j,i) - ( &
2372	( flux_r(k) + diss_r(k) &
2373	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
2374	+ ( flux_n(k) + diss_n(k) &
2375	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &
2376	+ ( ( flux_t(k) + diss_t(k) ) &
2377	- ( flux_d + diss_d ) &
2378	) * drho_air(k) * ddzw(k) &
2379	) + div * u(k,j,i)
2380
2381	flux_l_u(k,j,tn) = flux_r(k)
2382	diss_l_u(k,j,tn) = diss_r(k)
2383	flux_s_u(k,tn) = flux_n(k)
2384	diss_s_u(k,tn) = diss_n(k)
2385	!
2386	!-- Statistical Evaluation of u'u'. The factor has to be applied for
2387	!-- right evaluation when gallilei_trans = .T. .
2388	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
2389	+ ( flux_r(k) &
2390	* ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2391	/ ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &
2392	+ diss_r(k) &
2393	* ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2394	/ ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &
2395	) * weight_substep(intermediate_timestep_count)
2396	!
2397	!-- Statistical Evaluation of w'u'.
2398	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
2399	+ ( flux_t(k) &
2400	* ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2401	/ ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) &
2402	+ diss_t(k) &
2403	* ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2404	/ ( ABS( w_comp(k) ) + 1.0E-20_wp ) &
2405	) * weight_substep(intermediate_timestep_count)
2406	ENDDO
2407
2408	DO k = nzb_max_l+1, nzt
2409
2410	flux_d = flux_t(k-1)
2411	diss_d = diss_t(k-1)
2412	!
2413	!-- Calculate the divergence of the velocity field. A respective
2414	!-- correction is needed to overcome numerical instabilities introduced
2415	!-- by a not sufficient reduction of divergences near topography.
2416	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
2417	+ ( v_comp(k) + gv - ( v(k,j,i) + v(k,j,i-1) ) ) * ddy &
2418	+ ( w_comp(k) * rho_air_zw(k) &
2419	- w_comp(k-1) * rho_air_zw(k-1) &
2420	) * drho_air(k) * ddzw(k) &
2421	) * 0.5_wp
2422
2423	tend(k,j,i) = tend(k,j,i) - ( &
2424	( flux_r(k) + diss_r(k) &
2425	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
2426	+ ( flux_n(k) + diss_n(k) &
2427	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &
2428	+ ( ( flux_t(k) + diss_t(k) ) &
2429	- ( flux_d + diss_d ) &
2430	) * drho_air(k) * ddzw(k) &
2431	) + div * u(k,j,i)
2432
2433	flux_l_u(k,j,tn) = flux_r(k)
2434	diss_l_u(k,j,tn) = diss_r(k)
2435	flux_s_u(k,tn) = flux_n(k)
2436	diss_s_u(k,tn) = diss_n(k)
2437	!
2438	!-- Statistical Evaluation of u'u'. The factor has to be applied for
2439	!-- right evaluation when gallilei_trans = .T. .
2440	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
2441	+ ( flux_r(k) &
2442	* ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2443	/ ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &
2444	+ diss_r(k) &
2445	* ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2446	/ ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &
2447	) * weight_substep(intermediate_timestep_count)
2448	!
2449	!-- Statistical Evaluation of w'u'.
2450	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
2451	+ ( flux_t(k) &
2452	* ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2453	/ ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) &
2454	+ diss_t(k) &
2455	* ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2456	/ ( ABS( w_comp(k) ) + 1.0E-20_wp ) &
2457	) * weight_substep(intermediate_timestep_count)
2458	ENDDO
2459
2460
2461
2462	END SUBROUTINE advec_u_ws_ij
2463
2464
2465
2466	!-----------------------------------------------------------------------------!
2467	! Description:
2468	! ------------
2469	!> Advection of v-component - Call for grid point i,j
2470	!-----------------------------------------------------------------------------!
2471	SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn )
2472
2473
2474	INTEGER(iwp) :: i !< grid index along x-direction
2475	INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread
2476	INTEGER(iwp) :: j !< grid index along y-direction
2477	INTEGER(iwp) :: k !< grid index along z-direction
2478	INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults
2479	INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults
2480	INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults
2481	INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
2482	INTEGER(iwp) :: tn !< number of OpenMP thread
2483
2484	REAL(wp) :: ibit9 !< flag indicating 1st-order scheme along x-direction
2485	REAL(wp) :: ibit10 !< flag indicating 3rd-order scheme along x-direction
2486	REAL(wp) :: ibit11 !< flag indicating 5th-order scheme along x-direction
2487	REAL(wp) :: ibit12 !< flag indicating 1st-order scheme along y-direction
2488	REAL(wp) :: ibit13 !< flag indicating 3rd-order scheme along y-direction
2489	REAL(wp) :: ibit14 !< flag indicating 3rd-order scheme along y-direction
2490	REAL(wp) :: ibit15 !< flag indicating 1st-order scheme along z-direction
2491	REAL(wp) :: ibit16 !< flag indicating 3rd-order scheme along z-direction
2492	REAL(wp) :: ibit17 !< flag indicating 3rd-order scheme along z-direction
2493	REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom
2494	REAL(wp) :: div !< divergence on v-grid
2495	REAL(wp) :: flux_d !< 6th-order flux at grid box bottom
2496	REAL(wp) :: gu !< Galilei-transformation velocity along x
2497	REAL(wp) :: gv !< Galilei-transformation velocity along y
2498	REAL(wp) :: v_comp_l !< advection velocity along y on leftmost grid point on subdomain
2499
2500	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box
2501	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box
2502	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !< discretized artificial dissipation at top of the grid box
2503	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !< discretized 6th-order flux at northward-side of the grid box
2504	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box
2505	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6th-order flux at top of the grid box
2506	REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !< advection velocity along x
2507	REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !< advection velocity along y
2508	REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< advection velocity along z
2509	!
2510	!-- Used local modified copy of nzb_max (used to degrade order of
2511	!-- discretization) at non-cyclic boundaries. Modify only at relevant points
2512	!-- instead of the entire subdomain. This should lead to better
2513	!-- load balance between boundary and non-boundary PEs.
2514	IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR. &
2515	( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i >= nxr - 2 .OR. &
2516	( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j <= nys + 2 .OR. &
2517	( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j >= nyn - 2 ) THEN
2518	nzb_max_l = nzt
2519	ELSE
2520	nzb_max_l = nzb_max
2521	END IF
2522
2523	gu = 2.0_wp * u_gtrans
2524	gv = 2.0_wp * v_gtrans
2525
2526	!
2527	!-- Compute leftside fluxes for the respective boundary.
2528	IF ( i == i_omp ) THEN
2529
2530	DO k = nzb+1, nzb_max_l
2531
2532	ibit11 = REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp )
2533	ibit10 = REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp )
2534	ibit9 = REAL( IBITS(advc_flags_m(k,j,i-1),9,1), KIND = wp )
2535
2536	u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu
2537	flux_l_v(k,j,tn) = u_comp(k) * ( &
2538	( 37.0_wp * ibit11 * adv_mom_5 &
2539	+ 7.0_wp * ibit10 * adv_mom_3 &
2540	+ ibit9 * adv_mom_1 &
2541	) * &
2542	( v(k,j,i) + v(k,j,i-1) ) &
2543	- ( 8.0_wp * ibit11 * adv_mom_5 &
2544	+ ibit10 * adv_mom_3 &
2545	) * &
2546	( v(k,j,i+1) + v(k,j,i-2) ) &
2547	+ ( ibit11 * adv_mom_5 &
2548	) * &
2549	( v(k,j,i+2) + v(k,j,i-3) ) &
2550	)
2551
2552	diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( &
2553	( 10.0_wp * ibit11 * adv_mom_5 &
2554	+ 3.0_wp * ibit10 * adv_mom_3 &
2555	+ ibit9 * adv_mom_1 &
2556	) * &
2557	( v(k,j,i) - v(k,j,i-1) ) &
2558	- ( 5.0_wp * ibit11 * adv_mom_5 &
2559	+ ibit10 * adv_mom_3 &
2560	) * &
2561	( v(k,j,i+1) - v(k,j,i-2) ) &
2562	+ ( ibit11 * adv_mom_5 &
2563	) * &
2564	( v(k,j,i+2) - v(k,j,i-3) ) &
2565	)
2566
2567	ENDDO
2568
2569	DO k = nzb_max_l+1, nzt
2570
2571	u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu
2572	flux_l_v(k,j,tn) = u_comp(k) * ( &
2573	37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) &
2574	- 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) &
2575	+ ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
2576	diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( &
2577	10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) &
2578	- 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) &
2579	+ ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
2580
2581	ENDDO
2582
2583	ENDIF
2584	!
2585	!-- Compute southside fluxes for the respective boundary.
2586	IF ( j == nysv ) THEN
2587
2588	DO k = nzb+1, nzb_max_l
2589
2590	ibit14 = REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp )
2591	ibit13 = REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp )
2592	ibit12 = REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp )
2593
2594	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
2595	flux_s_v(k,tn) = v_comp_l * ( &
2596	( 37.0_wp * ibit14 * adv_mom_5 &
2597	+ 7.0_wp * ibit13 * adv_mom_3 &
2598	+ ibit12 * adv_mom_1 &
2599	) * &
2600	( v(k,j,i) + v(k,j-1,i) ) &
2601	- ( 8.0_wp * ibit14 * adv_mom_5 &
2602	+ ibit13 * adv_mom_3 &
2603	) * &
2604	( v(k,j+1,i) + v(k,j-2,i) ) &
2605	+ ( ibit14 * adv_mom_5 &
2606	) * &
2607	( v(k,j+2,i) + v(k,j-3,i) ) &
2608	)
2609
2610	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
2611	( 10.0_wp * ibit14 * adv_mom_5 &
2612	+ 3.0_wp * ibit13 * adv_mom_3 &
2613	+ ibit12 * adv_mom_1 &
2614	) * &
2615	( v(k,j,i) - v(k,j-1,i) ) &
2616	- ( 5.0_wp * ibit14 * adv_mom_5 &
2617	+ ibit13 * adv_mom_3 &
2618	) * &
2619	( v(k,j+1,i) - v(k,j-2,i) ) &
2620	+ ( ibit14 * adv_mom_5 &
2621	) * &
2622	( v(k,j+2,i) - v(k,j-3,i) ) &
2623	)
2624
2625	ENDDO
2626
2627	DO k = nzb_max_l+1, nzt
2628
2629	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
2630	flux_s_v(k,tn) = v_comp_l * ( &
2631	37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) &
2632	- 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) &
2633	+ ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
2634	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
2635	10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) &
2636	- 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) &
2637	+ ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
2638
2639	ENDDO
2640
2641	ENDIF
2642	!
2643	!-- Now compute the fluxes and tendency terms for the horizontal and
2644	!-- verical parts.
2645	DO k = nzb+1, nzb_max_l
2646
2647	ibit11 = REAL( IBITS(advc_flags_m(k,j,i),11,1), KIND = wp )
2648	ibit10 = REAL( IBITS(advc_flags_m(k,j,i),10,1), KIND = wp )
2649	ibit9 = REAL( IBITS(advc_flags_m(k,j,i),9,1), KIND = wp )
2650
2651	u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
2652	flux_r(k) = u_comp(k) * ( &
2653	( 37.0_wp * ibit11 * adv_mom_5 &
2654	+ 7.0_wp * ibit10 * adv_mom_3 &
2655	+ ibit9 * adv_mom_1 &
2656	) * &
2657	( v(k,j,i+1) + v(k,j,i) ) &
2658	- ( 8.0_wp * ibit11 * adv_mom_5 &
2659	+ ibit10 * adv_mom_3 &
2660	) * &
2661	( v(k,j,i+2) + v(k,j,i-1) ) &
2662	+ ( ibit11 * adv_mom_5 &
2663	) * &
2664	( v(k,j,i+3) + v(k,j,i-2) ) &
2665	)
2666
2667	diss_r(k) = - ABS( u_comp(k) ) * ( &
2668	( 10.0_wp * ibit11 * adv_mom_5 &
2669	+ 3.0_wp * ibit10 * adv_mom_3 &
2670	+ ibit9 * adv_mom_1 &
2671	) * &
2672	( v(k,j,i+1) - v(k,j,i) ) &
2673	- ( 5.0_wp * ibit11 * adv_mom_5 &
2674	+ ibit10 * adv_mom_3 &
2675	) * &
2676	( v(k,j,i+2) - v(k,j,i-1) ) &
2677	+ ( ibit11 * adv_mom_5 &
2678	) * &
2679	( v(k,j,i+3) - v(k,j,i-2) ) &
2680	)
2681
2682	ibit14 = REAL( IBITS(advc_flags_m(k,j,i),14,1), KIND = wp )
2683	ibit13 = REAL( IBITS(advc_flags_m(k,j,i),13,1), KIND = wp )
2684	ibit12 = REAL( IBITS(advc_flags_m(k,j,i),12,1), KIND = wp )
2685
2686
2687	v_comp(k) = v(k,j+1,i) + v(k,j,i)
2688	flux_n(k) = ( v_comp(k) - gv ) * ( &
2689	( 37.0_wp * ibit14 * adv_mom_5 &
2690	+ 7.0_wp * ibit13 * adv_mom_3 &
2691	+ ibit12 * adv_mom_1 &
2692	) * &
2693	( v(k,j+1,i) + v(k,j,i) ) &
2694	- ( 8.0_wp * ibit14 * adv_mom_5 &
2695	+ ibit13 * adv_mom_3 &
2696	) * &
2697	( v(k,j+2,i) + v(k,j-1,i) ) &
2698	+ ( ibit14 * adv_mom_5 &
2699	) * &
2700	( v(k,j+3,i) + v(k,j-2,i) ) &
2701	)
2702
2703	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
2704	( 10.0_wp * ibit14 * adv_mom_5 &
2705	+ 3.0_wp * ibit13 * adv_mom_3 &
2706	+ ibit12 * adv_mom_1 &
2707	) * &
2708	( v(k,j+1,i) - v(k,j,i) ) &
2709	- ( 5.0_wp * ibit14 * adv_mom_5 &
2710	+ ibit13 * adv_mom_3 &
2711	) * &
2712	( v(k,j+2,i) - v(k,j-1,i) ) &
2713	+ ( ibit14 * adv_mom_5 &
2714	) * &
2715	( v(k,j+3,i) - v(k,j-2,i) ) &
2716	)
2717	ENDDO
2718
2719	DO k = nzb_max_l+1, nzt
2720
2721	u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
2722	flux_r(k) = u_comp(k) * ( &
2723	37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) &
2724	- 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) &
2725	+ ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
2726
2727	diss_r(k) = - ABS( u_comp(k) ) * ( &
2728	10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) &
2729	- 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) &
2730	+ ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
2731
2732
2733	v_comp(k) = v(k,j+1,i) + v(k,j,i)
2734	flux_n(k) = ( v_comp(k) - gv ) * ( &
2735	37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) &
2736	- 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) &
2737	+ ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
2738
2739	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
2740	10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) &
2741	- 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) &
2742	+ ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
2743	ENDDO
2744	!
2745	!-- Now, compute vertical fluxes. Split loop into a part treating the
2746	!-- lowest 2 grid points with indirect indexing, a main loop without
2747	!-- indirect indexing, and a loop for the uppermost 2 grip points with
2748	!-- indirect indexing. This allows better vectorization for the main loop.
2749	!-- First, compute the flux at model surface, which need has to be
2750	!-- calculated explicetely for the tendency at
2751	!-- the first w-level. For topography wall this is done implicitely by
2752	!-- advc_flags_m.
2753	flux_t(nzb) = 0.0_wp
2754	diss_t(nzb) = 0.0_wp
2755	w_comp(nzb) = 0.0_wp
2756
2757	DO k = nzb+1, nzb+2
2758	!
2759	!-- k index has to be modified near bottom and top, else array
2760	!-- subscripts will be exceeded.
2761	ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp )
2762	ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp )
2763	ibit15 = REAL( IBITS(advc_flags_m(k,j,i),15,1), KIND = wp )
2764
2765	k_ppp = k + 3 * ibit17
2766	k_pp = k + 2 * ( 1 - ibit15 )
2767	k_mm = k - 2 * ibit17
2768
2769	w_comp(k) = w(k,j-1,i) + w(k,j,i)
2770	flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &
2771	( 37.0_wp * ibit17 * adv_mom_5 &
2772	+ 7.0_wp * ibit16 * adv_mom_3 &
2773	+ ibit15 * adv_mom_1 &
2774	) * &
2775	( v(k+1,j,i) + v(k,j,i) ) &
2776	- ( 8.0_wp * ibit17 * adv_mom_5 &
2777	+ ibit16 * adv_mom_3 &
2778	) * &
2779	( v(k_pp,j,i) + v(k-1,j,i) ) &
2780	+ ( ibit17 * adv_mom_5 &
2781	) * &
2782	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
2783	)
2784
2785	diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &
2786	( 10.0_wp * ibit17 * adv_mom_5 &
2787	+ 3.0_wp * ibit16 * adv_mom_3 &
2788	+ ibit15 * adv_mom_1 &
2789	) * &
2790	( v(k+1,j,i) - v(k,j,i) ) &
2791	- ( 5.0_wp * ibit17 * adv_mom_5 &
2792	+ ibit16 * adv_mom_3 &
2793	) * &
2794	( v(k_pp,j,i) - v(k-1,j,i) ) &
2795	+ ( ibit17 * adv_mom_5 &
2796	) * &
2797	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
2798	)
2799	ENDDO
2800
2801	DO k = nzb+3, nzt-2
2802
2803	ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp )
2804	ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp )
2805	ibit15 = REAL( IBITS(advc_flags_m(k,j,i),15,1), KIND = wp )
2806
2807	w_comp(k) = w(k,j-1,i) + w(k,j,i)
2808	flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &
2809	( 37.0_wp * ibit17 * adv_mom_5 &
2810	+ 7.0_wp * ibit16 * adv_mom_3 &
2811	+ ibit15 * adv_mom_1 &
2812	) * &
2813	( v(k+1,j,i) + v(k,j,i) ) &
2814	- ( 8.0_wp * ibit17 * adv_mom_5 &
2815	+ ibit16 * adv_mom_3 &
2816	) * &
2817	( v(k+2,j,i) + v(k-1,j,i) ) &
2818	+ ( ibit17 * adv_mom_5 &
2819	) * &
2820	( v(k+3,j,i) + v(k-2,j,i) ) &
2821	)
2822
2823	diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &
2824	( 10.0_wp * ibit17 * adv_mom_5 &
2825	+ 3.0_wp * ibit16 * adv_mom_3 &
2826	+ ibit15 * adv_mom_1 &
2827	) * &
2828	( v(k+1,j,i) - v(k,j,i) ) &
2829	- ( 5.0_wp * ibit17 * adv_mom_5 &
2830	+ ibit16 * adv_mom_3 &
2831	) * &
2832	( v(k+2,j,i) - v(k-1,j,i) ) &
2833	+ ( ibit17 * adv_mom_5 &
2834	) * &
2835	( v(k+3,j,i) - v(k-2,j,i) ) &
2836	)
2837	ENDDO
2838
2839	DO k = nzt-1, nzt
2840	!
2841	!-- k index has to be modified near bottom and top, else array
2842	!-- subscripts will be exceeded.
2843	ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp )
2844	ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp )
2845	ibit15 = REAL( IBITS(advc_flags_m(k,j,i),15,1), KIND = wp )
2846
2847	k_ppp = k + 3 * ibit17
2848	k_pp = k + 2 * ( 1 - ibit15 )
2849	k_mm = k - 2 * ibit17
2850
2851	w_comp(k) = w(k,j-1,i) + w(k,j,i)
2852	flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &
2853	( 37.0_wp * ibit17 * adv_mom_5 &
2854	+ 7.0_wp * ibit16 * adv_mom_3 &
2855	+ ibit15 * adv_mom_1 &
2856	) * &
2857	( v(k+1,j,i) + v(k,j,i) ) &
2858	- ( 8.0_wp * ibit17 * adv_mom_5 &
2859	+ ibit16 * adv_mom_3 &
2860	) * &
2861	( v(k_pp,j,i) + v(k-1,j,i) ) &
2862	+ ( ibit17 * adv_mom_5 &
2863	) * &
2864	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
2865	)
2866
2867	diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &
2868	( 10.0_wp * ibit17 * adv_mom_5 &
2869	+ 3.0_wp * ibit16 * adv_mom_3 &
2870	+ ibit15 * adv_mom_1 &
2871	) * &
2872	( v(k+1,j,i) - v(k,j,i) ) &
2873	- ( 5.0_wp * ibit17 * adv_mom_5 &
2874	+ ibit16 * adv_mom_3 &
2875	) * &
2876	( v(k_pp,j,i) - v(k-1,j,i) ) &
2877	+ ( ibit17 * adv_mom_5 &
2878	) * &
2879	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
2880	)
2881	ENDDO
2882
2883	DO k = nzb+1, nzb_max_l
2884
2885	flux_d = flux_t(k-1)
2886	diss_d = diss_t(k-1)
2887
2888	ibit11 = REAL( IBITS(advc_flags_m(k,j,i),11,1), KIND = wp )
2889	ibit10 = REAL( IBITS(advc_flags_m(k,j,i),10,1), KIND = wp )
2890	ibit9 = REAL( IBITS(advc_flags_m(k,j,i),9,1), KIND = wp )
2891
2892	ibit14 = REAL( IBITS(advc_flags_m(k,j,i),14,1), KIND = wp )
2893	ibit13 = REAL( IBITS(advc_flags_m(k,j,i),13,1), KIND = wp )
2894	ibit12 = REAL( IBITS(advc_flags_m(k,j,i),12,1), KIND = wp )
2895
2896	ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp )
2897	ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp )
2898	ibit15 = REAL( IBITS(advc_flags_m(k,j,i),15,1), KIND = wp )
2899	!
2900	!-- Calculate the divergence of the velocity field. A respective
2901	!-- correction is needed to overcome numerical instabilities introduced
2902	!-- by a not sufficient reduction of divergences near topography.
2903	div = ( ( ( u_comp(k) + gu ) &
2904	* ( ibit9 + ibit10 + ibit11 ) &
2905	- ( u(k,j-1,i) + u(k,j,i) ) &
2906	* ( &
2907	REAL( IBITS(advc_flags_m(k,j,i-1),9,1), KIND = wp ) &
2908	+ REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp ) &
2909	+ REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp ) &
2910	) &
2911	) * ddx &
2912	+ ( v_comp(k) &
2913	* ( ibit12 + ibit13 + ibit14 ) &
2914	- ( v(k,j,i) + v(k,j-1,i) ) &
2915	* ( &
2916	REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp ) &
2917	+ REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp ) &
2918	+ REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp ) &
2919	) &
2920	) * ddy &
2921	+ ( w_comp(k) * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )&
2922	- w_comp(k-1) * rho_air_zw(k-1) &
2923	* ( &
2924	REAL( IBITS(advc_flags_m(k-1,j,i),15,1), KIND = wp ) &
2925	+ REAL( IBITS(advc_flags_m(k-1,j,i),16,1), KIND = wp ) &
2926	+ REAL( IBITS(advc_flags_m(k-1,j,i),17,1), KIND = wp ) &
2927	) &
2928	) * drho_air(k) * ddzw(k) &
2929	) * 0.5_wp
2930
2931	tend(k,j,i) = tend(k,j,i) - ( &
2932	( flux_r(k) + diss_r(k) &
2933	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
2934	+ ( flux_n(k) + diss_n(k) &
2935	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &
2936	+ ( ( flux_t(k) + diss_t(k) ) &
2937	- ( flux_d + diss_d ) &
2938	) * drho_air(k) * ddzw(k) &
2939	) + v(k,j,i) * div
2940
2941	flux_l_v(k,j,tn) = flux_r(k)
2942	diss_l_v(k,j,tn) = diss_r(k)
2943	flux_s_v(k,tn) = flux_n(k)
2944	diss_s_v(k,tn) = diss_n(k)
2945	!
2946	!-- Statistical Evaluation of v'v'. The factor has to be applied for
2947	!-- right evaluation when gallilei_trans = .T. .
2948	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
2949	+ ( flux_n(k) &
2950	* ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2951	/ ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &
2952	+ diss_n(k) &
2953	* ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2954	/ ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &
2955	) * weight_substep(intermediate_timestep_count)
2956	!
2957	!-- Statistical Evaluation of w'u'.
2958	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
2959	+ ( flux_t(k) &
2960	* ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2961	/ ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) &
2962	+ diss_t(k) &
2963	* ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2964	/ ( ABS( w_comp(k) ) + 1.0E-20_wp ) &
2965	) * weight_substep(intermediate_timestep_count)
2966
2967	ENDDO
2968
2969	DO k = nzb_max_l+1, nzt
2970
2971	flux_d = flux_t(k-1)
2972	diss_d = diss_t(k-1)
2973	!
2974	!-- Calculate the divergence of the velocity field. A respective
2975	!-- correction is needed to overcome numerical instabilities introduced
2976	!-- by a not sufficient reduction of divergences near topography.
2977	div = ( ( u_comp(k) + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
2978	+ ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
2979	+ ( w_comp(k) * rho_air_zw(k) &
2980	- w_comp(k-1) * rho_air_zw(k-1) &
2981	) * drho_air(k) * ddzw(k) &
2982	) * 0.5_wp
2983
2984	tend(k,j,i) = tend(k,j,i) - ( &
2985	( flux_r(k) + diss_r(k) &
2986	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
2987	+ ( flux_n(k) + diss_n(k) &
2988	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &
2989	+ ( ( flux_t(k) + diss_t(k) ) &
2990	- ( flux_d + diss_d ) &
2991	) * drho_air(k) * ddzw(k) &
2992	) + v(k,j,i) * div
2993
2994	flux_l_v(k,j,tn) = flux_r(k)
2995	diss_l_v(k,j,tn) = diss_r(k)
2996	flux_s_v(k,tn) = flux_n(k)
2997	diss_s_v(k,tn) = diss_n(k)
2998	!
2999	!-- Statistical Evaluation of v'v'. The factor has to be applied for
3000	!-- right evaluation when gallilei_trans = .T. .
3001	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
3002	+ ( flux_n(k) &
3003	* ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
3004	/ ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &
3005	+ diss_n(k) &
3006	* ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
3007	/ ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &
3008	) * weight_substep(intermediate_timestep_count)
3009	!
3010	!-- Sta