Home

Context Navigation

advec_ws.f90 @ 3870

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