Home

Context Navigation

advec_ws.f90 @ 4598

Last change on this file since 4598 was 4581, checked in by suehring, 10 months ago
mesoscale nesting: omit explicit pressure forcing via geostrophic wind components; chemistry: enable profile output of vertical fluxes; urban-surface: bugfix in initialization in case of cyclic_fill
Property svn:keywords set to `Id`
File size: 364.6 KB

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