Home

Context Navigation

advec_ws.f90 @ 4324

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