Home

Context Navigation

advec_ws.f90 @ 3661

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