Home

Context Navigation

advec_ws.f90 @ 4318

Last change on this file since 4318 was 4318, checked in by Giersch, 5 years ago
Indirect indexing for calculating vertical fluxes close to boundaries is only used for loop indices where it is really necessary
Property svn:keywords set to `Id`
File size: 375.7 KB

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