Home

Context Navigation

advec_ws.f90 @ 4317

Last change on this file since 4317 was 4317, checked in by Giersch, 5 years ago
Comments revised/added, formatting improved, fluxes for u,v, and scalars are explicitly set to zero at nzt+1, fluxes of w-component are now calculated only until nzt-1 (Prognostic equation for w-velocity component ends at nzt-1)
Property svn:keywords set to `Id`
File size: 375.5 KB

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