Home

Context Navigation

advec_ws.f90 @ 4110

Last change on this file since 4110 was 4110, checked in by suehring, 4 years ago
last changes documented
Property svn:keywords set to `Id`
File size: 381.2 KB

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