source: palm/trunk/SOURCE/advec_ws.f90 @ 4109

Last change on this file since 4109 was 4109, checked in by suehring, 2 years ago

Control discretization of advection term: separate initialization of WS advection flags for momentum and scalars. In this context, resort the bits and do some minor formatting; Make initialization of scalar-advection flags more flexible, i.e. introduce an arguemnt list to indicate non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols); Introduce extended 'degradation zones', where horizontal advection of passive scalars is discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3 grid points). Even though no building is within the numerical stencil, first-order scheme is used. At fourth and fifth grid point the order of the horizontal advection scheme is successively upgraded. These extended degradation zones are used to avoid stationary numerical oscillations, which are responsible for high concentration maxima that may appear under shear-free stable conditions. Therefore, an additional 3D interger array used to store control flags is introduced; Change interface for scalar advection routine; Bugfix, avoid uninitialized value sk_num in vector version of WS scalar advection; Chemistry: Decycling boundary conditions are only set at the ghost points not on the prognostic grid points; Land-surface model: Relax checks for non-consistent initialization in case static or dynamic input is provided. For example, soil_temperature or deep_soil_temperature is not mandatory any more if dynamic input is available. Also, improper settings of x_type in namelist are only checked if no static file is available.

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