Home

Context Navigation

advec_ws.f90 @ 4109

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