Home

Context Navigation

advec_ws.f90 @ 3871

Last change on this file since 3871 was 3871, checked in by knoop, 4 years ago
Moving prognostic equations of bcm into bulk_cloud_model_mod
Property svn:keywords set to `Id`
File size: 353.8 KB

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