Home

Context Navigation

advec_ws.f90 @ 4329

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