Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Apr 26, 2020 3:57:55 PM (9 months ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/advec_ws.f90

    r4502 r4509  
    11!> @file advec_ws.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    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/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    … …  
    2524! -----------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4502 2020-04-17 16:14:16Z schwenkel
    2729! Implementation of ice microphysics
    2830!
    … …  
    3638!
    3739! 4466 2020-03-20 16:14:41Z suehring
    38 ! - vector branch further optimized (linear dependencies along z removed and
    39 !   loops are splitted)
    40 ! - topography closed channel flow with symmetric boundaries also implemented
    41 !   in vector branch
     40! - vector branch further optimized (linear dependencies along z removed and loops are splitted)
     41! - topography closed channel flow with symmetric boundaries also implemented in vector branch
    4242! - some formatting adjustments made and comments added
    4343! - cpu measures for vector branch added
    … …  
    5050!
    5151! 4360 2020-01-07 11:25:50Z suehring
    52 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    53 ! topography information used in wall_flags_static_0
     52! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     53! information used in wall_flags_static_0
    5454!
    5555! 4340 2019-12-16 08:17:03Z Giersch
    … …  
    6666!
    6767! 4327 2019-12-06 14:48:31Z Giersch
    68 ! Setting of advection flags for vertical fluxes of w revised, air density for
    69 ! vertical flux calculation of w at k=1 is considered now
     68! Setting of advection flags for vertical fluxes of w revised, air density for vertical flux
     69! calculation of w at k=1 is considered now
    7070!
    7171! 4325 2019-12-06 07:14:04Z Giersch
    72 ! Vertical fluxes of w are now set to zero at nzt and nzt+1, setting of
    73 ! advection flags for fluxes in z-direction revised, comments extended
     72! Vertical fluxes of w are now set to zero at nzt and nzt+1, setting of advection flags for fluxes
     73! in z-direction revised, comments extended
    7474!
    7575! 4324 2019-12-06 07:11:33Z Giersch
    76 ! Indirect indexing for calculating vertical fluxes close to boundaries is only
    77 ! used for loop indizes where it is really necessary
     76! Indirect indexing for calculating vertical fluxes close to boundaries is only used for loop
     77! indizes where it is really necessary
    7878!
    7979! 4317 2019-12-03 12:43:22Z Giersch
    80 ! Comments revised/added, formatting improved, fluxes for u,v, and scalars are
    81 ! explicitly set to zero at nzt+1, fluxes of w-component are now calculated only
    82 ! until nzt-1 (Prognostic equation for w-velocity component ends at nzt-1)
     80! Comments revised/added, formatting improved, fluxes for u,v, and scalars are explicitly set to
     81! zero at nzt+1, fluxes of w-component are now calculated only until nzt-1 (Prognostic equation for
     82! w-velocity component ends at nzt-1)
    8383!
    8484! 4204 2019-08-30 12:30:17Z knoop
    … …  
    8989!
    9090! 4110 2019-07-22 17:05:21Z suehring
    91 ! - Separate initialization of advection flags for momentum and scalars. In this
    92 !   context, resort the bits and do some minor formatting.
    93 ! - Make flag initialization for scalars more flexible, introduce an
    94 !   arguemnt list to indicate non-cyclic boundaries (required for decycled
    95 !   scalars such as chemical species or aerosols)
    96 ! - Introduce extended 'degradation zones', where horizontal advection of
    97 !   passive scalars is discretized by first-order scheme at all grid points
    98 !   that in the vicinity of buildings (<= 3 grid points). Even though no
    99 !   building is within the numerical stencil, first-order scheme is used.
    100 !   At fourth and fifth grid point the order of the horizontal advection scheme
    101 !   is successively upgraded.
    102 !   These extended degradation zones are used to avoid stationary numerical
    103 !   oscillations, which are responsible for high concentration maxima that may
    104 !   appear under shear-free stable conditions.
     91! - Separate initialization of advection flags for momentum and scalars. In this context, resort the
     92!   bits and do some minor formatting.
     93! - Make flag initialization for scalars more flexible, introduce an arguemnt list to indicate
     94!   non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols)
     95! - Introduce extended 'degradation zones', where horizontal advection of passive scalars is
     96!   discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3
     97!   grid points). Even though no building is within the numerical stencil, first-order scheme is
     98!   used. At fourth and fifth grid point the order of the horizontal advection scheme is
     99!   successively upgraded. These extended degradation zones are used to avoid stationary numerical
     100!   oscillations, which are responsible for high concentration maxima that may appear under
     101!   shear-free stable conditions.
    105102! - Change interface for scalar advection routine.
    106103! - Bugfix, avoid uninitialized value sk_num in vector version of scalar
    … …  
    108105!
    109106! 4109 2019-07-22 17:00:34Z suehring
    110 ! Implementation of a flux limiter according to Skamarock (2006) for the
    111 ! vertical scalar advection. Please note, this is only implemented for the
    112 ! cache-optimized version at the moment. Implementation for the vector-
    113 ! optimized version will follow after critical issues concerning
     107! Implementation of a flux limiter according to Skamarock (2006) for the vertical scalar advection.
     108! Please note, this is only implemented for the cache-optimized version at the moment.
     109! Implementation for the vector-optimized version will follow after critical issues concerning
    114110! vectorization are fixed.
    115111!
    … …  
    130126!
    131127! 3661 2019-01-08 18:22:50Z suehring
    132 ! - Minor bugfix in divergence correction (only has implications at
    133 !   downward-facing wall surfaces)
     128! - Minor bugfix in divergence correction (only has implications at downward-facing wall surfaces)
    134129! - Remove setting of Neumann condition for horizontal velocity variances
    135 ! - Split loops for tendency calculation and divergence correction in order to
    136 !   reduce bit queries
    137 ! - Introduce new parameter nzb_max_l to better control order degradation at
    138 !   non-cyclic boundaries
     130! - Split loops for tendency calculation and divergence correction in order to reduce bit queries
     131! - Introduce new parameter nzb_max_l to better control order degradation at non-cyclic boundaries
    139132!
    140133! 3655 2019-01-07 16:51:22Z knoop
    … …  
    151144! Description:
    152145! ------------
    153 !> Advection scheme for scalars and momentum using the flux formulation of
    154 !> Wicker and Skamarock 5th order. Additionally the module contains of a
    155 !> routine using for initialisation and steering of the statical evaluation.
    156 !> The computation of turbulent fluxes takes place inside the advection
    157 !> routines.
    158 !> Near non-cyclic boundaries the order of the applied advection scheme is
    159 !> degraded.
    160 !> A divergence correction is applied. It is necessary for topography, since
    161 !> the divergence is not sufficiently reduced, resulting in erroneous fluxes
    162 !> and could lead to numerical instabilities.
     146!> Advection scheme for scalars and momentum using the flux formulation of Wicker and Skamarock 5th
     147!> order. Additionally the module contains of a routine using for initialisation and steering of the
     148!> statical evaluation.
     149!> The computation of turbulent fluxes takes place inside the advection routines.
     150!> Near non-cyclic boundaries the order of the applied advection scheme is degraded.
     151!> A divergence correction is applied. It is necessary for topography, since the divergence is not
     152!> sufficiently reduced, resulting in erroneous fluxes and could lead to numerical instabilities.
    163153!>
    164154!> @todo Implement monotonic flux limiter also for vector version.
    165155!> @todo Move 3d arrays advc_flag, advc_flags_m from modules to advec_ws
    166156!> @todo Move arrays flux_l_x from modules to advec_ws
    167 !------------------------------------------------------------------------------!
     157!--------------------------------------------------------------------------------------------------!
    168158 MODULE advec_ws
    169159
    170     USE arrays_3d,                                                             &
    171         ONLY:  ddzu, ddzw, tend, u, v, w,                                      &
    172                drho_air, drho_air_zw, rho_air, rho_air_zw,                     &
    173                u_stokes_zu, v_stokes_zu,                                       &
    174                diss_l_diss, diss_l_e, diss_l_pt, diss_l_q,                     &
    175                diss_l_s, diss_l_u, diss_l_v, diss_l_w,                         &
    176                flux_l_diss, flux_l_e, flux_l_pt, flux_l_q, flux_l_s,           &
    177                flux_l_u, flux_l_v, flux_l_w,                                   &
    178                diss_s_diss, diss_s_e, diss_s_pt, diss_s_q, diss_s_s,           &
    179                diss_s_u, diss_s_v, diss_s_w,                                   &
    180                flux_s_diss, flux_s_e, flux_s_pt, flux_s_q, flux_s_s,           &
    181                flux_s_u, flux_s_v, flux_s_w
    182 
    183     USE control_parameters,                                                    &
    184         ONLY:  bc_dirichlet_l,                                                 &
    185                bc_dirichlet_n,                                                 &
    186                bc_dirichlet_r,                                                 &
    187                bc_dirichlet_s,                                                 &
    188                bc_radiation_l,                                                 &
    189                bc_radiation_n,                                                 &
    190                bc_radiation_r,                                                 &
    191                bc_radiation_s,                                                 &
    192                humidity,                                                       &
    193                loop_optimization,                                              &
    194                passive_scalar,                                                 &
    195                rans_tke_e,                                                     &
    196                symmetry_flag,                                                  &
    197                intermediate_timestep_count,                                    &
    198                u_gtrans,                                                       &
    199                v_gtrans,                                                       &
    200                ws_scheme_mom,                                                  &
    201                ws_scheme_sca,                                                  &
    202                dt_3d
    203 
    204     USE cpulog,                                                                &
    205         ONLY:  cpu_log,                                                        &
     160    USE arrays_3d,                                                                                 &
     161        ONLY:  ddzu, ddzw, tend, u, v, w,                                                          &
     162               diss_l_diss, diss_l_e, diss_l_pt, diss_l_q,                                         &
     163               diss_l_s, diss_l_u, diss_l_v, diss_l_w,                                             &
     164               diss_s_diss, diss_s_e, diss_s_pt, diss_s_q, diss_s_s,                               &
     165               diss_s_u, diss_s_v, diss_s_w,                                                       &
     166               drho_air, drho_air_zw, rho_air, rho_air_zw,                                         &
     167               flux_l_diss, flux_l_e, flux_l_pt, flux_l_q, flux_l_s,                               &
     168               flux_l_u, flux_l_v, flux_l_w,                                                       &
     169               flux_s_diss, flux_s_e, flux_s_pt, flux_s_q, flux_s_s,                               &
     170               flux_s_u, flux_s_v, flux_s_w,                                                       &
     171               u_stokes_zu, v_stokes_zu
     172
     173    USE control_parameters,                                                                        &
     174        ONLY:  bc_dirichlet_l,                                                                     &
     175               bc_dirichlet_n,                                                                     &
     176               bc_dirichlet_r,                                                                     &
     177               bc_dirichlet_s,                                                                     &
     178               bc_radiation_l,                                                                     &
     179               bc_radiation_n,                                                                     &
     180               bc_radiation_r,                                                                     &
     181               bc_radiation_s,                                                                     &
     182               dt_3d,                                                                              &
     183               humidity,                                                                           &
     184               intermediate_timestep_count,                                                        &
     185               loop_optimization,                                                                  &
     186               passive_scalar,                                                                     &
     187               rans_tke_e,                                                                         &
     188               symmetry_flag,                                                                      &
     189               u_gtrans,                                                                           &
     190               v_gtrans,                                                                           &
     191               ws_scheme_mom,                                                                      &
     192               ws_scheme_sca
     193
     194    USE cpulog,                                                                                    &
     195        ONLY:  cpu_log,                                                                            &
    206196               log_point_s
    207197
    208     USE exchange_horiz_mod,                                                    &
     198    USE exchange_horiz_mod,                                                                        &
    209199        ONLY:  exchange_horiz_int
    210200
    211     USE indices,                                                               &
    212         ONLY:  advc_flags_m,                                                   &
    213                advc_flags_s,                                                   &
    214                nbgp,                                                           &
    215                nx,                                                             &
    216                nxl,                                                            &
    217                nxlg,                                                           &
    218                nxlu,                                                           &
    219                nxr,                                                            &
    220                nxrg,                                                           &
    221                ny,                                                             &
    222                nyn,                                                            &
    223                nyng,                                                           &
    224                nys,                                                            &
    225                nysg,                                                           &
    226                nysv,                                                           &
    227                nzb,                                                            &
    228                nzb_max,                                                        &
    229                nzt,                                                            &
     201    USE indices,                                                                                   &
     202        ONLY:  advc_flags_m,                                                                       &
     203               advc_flags_s,                                                                       &
     204               nbgp,                                                                               &
     205               nx,                                                                                 &
     206               nxl,                                                                                &
     207               nxlg,                                                                               &
     208               nxlu,                                                                               &
     209               nxr,                                                                                &
     210               nxrg,                                                                               &
     211               ny,                                                                                 &
     212               nyn,                                                                                &
     213               nyng,                                                                               &
     214               nys,                                                                                &
     215               nysg,                                                                               &
     216               nysv,                                                                               &
     217               nzb,                                                                                &
     218               nzb_max,                                                                            &
     219               nzt,                                                                                &
    230220               wall_flags_total_0
    231221
    232     USE grid_variables,                                                        &
     222    USE grid_variables,                                                                            &
    233223        ONLY:  ddx, ddy
    234224
    235225    USE kinds
    236226
    237     USE pegrid,                                                                &
    238            ONLY:  threads_per_task
    239 
    240     USE statistics,                                                            &
    241         ONLY:  sums_salsa_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l,   &
    242                sums_wspts_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l,                &
    243                sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l,                &
    244                sums_wsqcs_ws_l, sums_wsqrs_ws_l,                               &
    245                sums_wsqis_ws_l, sums_wsnis_ws_l,                               &
    246                sums_wsncs_ws_l, sums_wsnrs_ws_l,                               &
    247                hom, weight_substep
     227    USE pegrid,                                                                                    &
     228        ONLY:  threads_per_task
     229
     230    USE statistics,                                                                                &
     231        ONLY:  hom, weight_substep,                                                                &
     232               sums_salsa_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l,                       &
     233               sums_wsncs_ws_l, sums_wsnrs_ws_l,                                                   &
     234               sums_wspts_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l,                                    &
     235               sums_wsqcs_ws_l, sums_wsqrs_ws_l,                                                   &
     236               sums_wsqis_ws_l, sums_wsnis_ws_l,                                                   &
     237               sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l
     238
     239
    248240
    249241    IMPLICIT NONE
    250 
     242   
    251243    REAL(wp) ::  adv_mom_1            !< 1/4 - constant used in 5th-order advection scheme for momentum advection (1st-order part)
    252244    REAL(wp) ::  adv_mom_3            !< 1/24 - constant used in 5th-order advection scheme for momentum advection (3rd-order part)
    … …  
    257249
    258250    PRIVATE
    259     PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, ws_init,          &
    260              ws_init_flags_momentum, ws_init_flags_scalar, ws_statistics
     251    PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, ws_init, ws_init_flags_momentum,      &
     252             ws_init_flags_scalar, ws_statistics
    261253
    262254    INTERFACE ws_init
    … …  
    299291
    300292
    301 !------------------------------------------------------------------------------!
     293!--------------------------------------------------------------------------------------------------!
    302294! Description:
    303295! ------------
    304296!> Initialization of WS-scheme
    305 !------------------------------------------------------------------------------!
    306     SUBROUTINE ws_init
    307 
    308 !
    309 !--    Set factors for scalar and momentum advection.
    310        adv_sca_5 = 1.0_wp /  60.0_wp
    311        adv_sca_3 = 1.0_wp /  12.0_wp
    312        adv_sca_1 = 1.0_wp /   2.0_wp
    313        adv_mom_5 = 1.0_wp / 120.0_wp
    314        adv_mom_3 = 1.0_wp /  24.0_wp
    315        adv_mom_1 = 1.0_wp /   4.0_wp
    316 !
    317 !--    Arrays needed for statical evaluation of fluxes.
    318        IF ( ws_scheme_mom )  THEN
    319 
    320           ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1),            &
    321                     sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1),            &
    322                     sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1),             &
    323                     sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1),             &
    324                     sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    325 
    326           sums_wsus_ws_l = 0.0_wp
    327           sums_wsvs_ws_l = 0.0_wp
    328           sums_us2_ws_l  = 0.0_wp
    329           sums_vs2_ws_l  = 0.0_wp
    330           sums_ws2_ws_l  = 0.0_wp
    331 
     297!--------------------------------------------------------------------------------------------------!
     298 SUBROUTINE ws_init
     299
     300!
     301!-- Set factors for scalar and momentum advection.
     302    adv_sca_5 = 1.0_wp /  60.0_wp
     303    adv_sca_3 = 1.0_wp /  12.0_wp
     304    adv_sca_1 = 1.0_wp /   2.0_wp
     305    adv_mom_5 = 1.0_wp / 120.0_wp
     306    adv_mom_3 = 1.0_wp /  24.0_wp
     307    adv_mom_1 = 1.0_wp /   4.0_wp
     308!
     309!-- Arrays needed for statical evaluation of fluxes.
     310    IF ( ws_scheme_mom )  THEN
     311
     312       ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1),                                   &
     313                 sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1),                                   &
     314                 sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1),                                    &
     315                 sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1),                                    &
     316                 sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     317
     318       sums_wsus_ws_l = 0.0_wp
     319       sums_wsvs_ws_l = 0.0_wp
     320       sums_us2_ws_l  = 0.0_wp
     321       sums_vs2_ws_l  = 0.0_wp
     322       sums_ws2_ws_l  = 0.0_wp
     323
     324    ENDIF
     325
     326    IF ( ws_scheme_sca )  THEN
     327
     328       ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     329       sums_wspts_ws_l = 0.0_wp
     330
     331       IF ( humidity  )  THEN
     332          ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     333          sums_wsqs_ws_l = 0.0_wp
    332334       ENDIF
    333335
     336       IF ( passive_scalar )  THEN
     337          ALLOCATE( sums_wsss_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     338          sums_wsss_ws_l = 0.0_wp
     339       ENDIF
     340
     341    ENDIF
     342
     343!
     344!-- Arrays needed for reasons of speed optimization
     345    IF ( ws_scheme_mom )  THEN
     346
     347       ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1),                                         &
     348                 flux_s_v(nzb+1:nzt,0:threads_per_task-1),                                         &
     349                 flux_s_w(nzb+1:nzt,0:threads_per_task-1),                                         &
     350                 diss_s_u(nzb+1:nzt,0:threads_per_task-1),                                         &
     351                 diss_s_v(nzb+1:nzt,0:threads_per_task-1),                                         &
     352                 diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
     353       ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                                 &
     354                 flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                                 &
     355                 flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                                 &
     356                 diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                                 &
     357                 diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                                 &
     358                 diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
     359
     360    ENDIF
     361!
     362!-- For the vector version the buffer arrays for scalars are not necessary, since internal arrays
     363!-- are used in the vector version.
     364    IF ( loop_optimization /= 'vector' )  THEN
    334365       IF ( ws_scheme_sca )  THEN
    335366
    336           ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    337           sums_wspts_ws_l = 0.0_wp
    338 
    339           IF ( humidity  )  THEN
    340              ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    341              sums_wsqs_ws_l = 0.0_wp
     367          ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1),                                     &
     368                    flux_s_e(nzb+1:nzt,0:threads_per_task-1),                                      &
     369                    diss_s_pt(nzb+1:nzt,0:threads_per_task-1),                                     &
     370                    diss_s_e(nzb+1:nzt,0:threads_per_task-1) )
     371          ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                             &
     372                    flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                              &
     373                    diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                             &
     374                    diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
     375
     376          IF ( rans_tke_e )  THEN
     377             ALLOCATE( flux_s_diss(nzb+1:nzt,0:threads_per_task-1),                                &
     378                       diss_s_diss(nzb+1:nzt,0:threads_per_task-1) )
     379             ALLOCATE( flux_l_diss(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                        &
     380                       diss_l_diss(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    342381          ENDIF
    343382
     383          IF ( humidity )  THEN
     384             ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1),                                   &
     385                       diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
     386             ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                           &
     387                       diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
     388          ENDIF
     389
    344390          IF ( passive_scalar )  THEN
    345              ALLOCATE( sums_wsss_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    346              sums_wsss_ws_l = 0.0_wp
     391             ALLOCATE( flux_s_s(nzb+1:nzt,0:threads_per_task-1),                                   &
     392                       diss_s_s(nzb+1:nzt,0:threads_per_task-1) )
     393             ALLOCATE( flux_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1),                           &
     394                       diss_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    347395          ENDIF
    348396
    349397       ENDIF
    350 
    351 !
    352 !--    Arrays needed for reasons of speed optimization
    353        IF ( ws_scheme_mom )  THEN
    354 
    355           ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1),               &
    356                     flux_s_v(nzb+1:nzt,0:threads_per_task-1),               &
    357                     flux_s_w(nzb+1:nzt,0:threads_per_task-1),               &
    358                     diss_s_u(nzb+1:nzt,0:threads_per_task-1),               &
    359                     diss_s_v(nzb+1:nzt,0:threads_per_task-1),               &
    360                     diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
    361           ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
    362                     flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
    363                     flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
    364                     diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
    365                     diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
    366                     diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    367 
    368        ENDIF
    369 !
    370 !--    For the vector version the buffer arrays for scalars are not necessary,
    371 !--    since internal arrays are used in the vector version.
    372        IF ( loop_optimization /= 'vector' )  THEN
    373           IF ( ws_scheme_sca )  THEN
    374 
    375              ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1),              &
    376                        flux_s_e(nzb+1:nzt,0:threads_per_task-1),               &
    377                        diss_s_pt(nzb+1:nzt,0:threads_per_task-1),              &
    378                        diss_s_e(nzb+1:nzt,0:threads_per_task-1) )
    379              ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),      &
    380                        flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
    381                        diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),      &
    382                        diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    383 
    384              IF ( rans_tke_e )  THEN
    385                 ALLOCATE( flux_s_diss(nzb+1:nzt,0:threads_per_task-1),         &
    386                           diss_s_diss(nzb+1:nzt,0:threads_per_task-1) )
    387                 ALLOCATE( flux_l_diss(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
    388                           diss_l_diss(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    389              ENDIF
    390 
    391              IF ( humidity )  THEN
    392                 ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1),            &
    393                           diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
    394                 ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
    395                           diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    396              ENDIF
    397 
    398              IF ( passive_scalar )  THEN
    399                 ALLOCATE( flux_s_s(nzb+1:nzt,0:threads_per_task-1),            &
    400                           diss_s_s(nzb+1:nzt,0:threads_per_task-1) )
    401                 ALLOCATE( flux_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
    402                           diss_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    403              ENDIF
    404 
    405           ENDIF
    406        ENDIF
    407 !
    408 !--    Initialize the flag arrays controlling degradation near walls, i.e.
    409 !--    to decrease the numerical stencil appropriately. The order of the scheme
    410 !--    is degraded near solid walls as well as near non-cyclic inflow and outflow
    411 !--    boundaries. Do this separately for momentum and scalars.
    412        IF ( ws_scheme_mom )  CALL ws_init_flags_momentum
    413 
    414        IF ( ws_scheme_sca )  THEN
    415           ALLOCATE( advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    416           advc_flags_s = 0
    417           CALL ws_init_flags_scalar( bc_dirichlet_l  .OR.  bc_radiation_l,     &
    418                                      bc_dirichlet_n  .OR.  bc_radiation_n,     &
    419                                      bc_dirichlet_r  .OR.  bc_radiation_r,     &
    420                                      bc_dirichlet_s  .OR.  bc_radiation_s,     &
    421                                      advc_flags_s )
    422        ENDIF
    423 
    424     END SUBROUTINE ws_init
    425 
    426 !------------------------------------------------------------------------------!
     398    ENDIF
     399!
     400!-- Initialize the flag arrays controlling degradation near walls, i.e. to decrease the numerical
     401!-- stencil appropriately. The order of the scheme is degraded near solid walls as well as near
     402!-- non-cyclic inflow and outflow boundaries. Do this separately for momentum and scalars.
     403    IF ( ws_scheme_mom )  CALL ws_init_flags_momentum
     404
     405    IF ( ws_scheme_sca )  THEN
     406       ALLOCATE( advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     407       advc_flags_s = 0
     408       CALL ws_init_flags_scalar( bc_dirichlet_l  .OR.  bc_radiation_l,                            &
     409                                  bc_dirichlet_n  .OR.  bc_radiation_n,                            &
     410                                  bc_dirichlet_r  .OR.  bc_radiation_r,                            &
     411                                  bc_dirichlet_s  .OR.  bc_radiation_s,                            &
     412                                  advc_flags_s )
     413    ENDIF
     414
     415 END SUBROUTINE ws_init
     416
     417!--------------------------------------------------------------------------------------------------!
    427418! Description:
    428419! ------------
    429 !> Initialization of flags to control the order of the advection scheme near
    430 !> solid walls and non-cyclic inflow boundaries, where the order is sucessively
    431 !> degraded.
    432 !------------------------------------------------------------------------------!
    433     SUBROUTINE ws_init_flags_momentum
    434 
    435 
    436        INTEGER(iwp) ::  i     !< index variable along x
    437        INTEGER(iwp) ::  j     !< index variable along y
    438        INTEGER(iwp) ::  k     !< index variable along z
    439        INTEGER(iwp) ::  k_mm  !< dummy index along z
    440        INTEGER(iwp) ::  k_pp  !< dummy index along z
    441        INTEGER(iwp) ::  k_ppp !< dummy index along z
    442 
    443        LOGICAL      ::  flag_set !< steering variable for advection flags
    444 
    445        ALLOCATE( advc_flags_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    446        advc_flags_m = 0
    447 !
    448 !--    Set advc_flags_m to steer the degradation of the advection scheme in advec_ws
    449 !--    near topography, inflow- and outflow boundaries as well as bottom and
    450 !--    top of model domain. advc_flags_m remains zero for all non-prognostic
    451 !--    grid points.
    452 !--    u-component
    453        DO  i = nxl, nxr
    454           DO  j = nys, nyn
    455              DO  k = nzb+1, nzt
    456 !
    457 !--             At first, set flags to WS1.
    458 !--             Since fluxes are swapped in advec_ws.f90, this is necessary to
    459 !--             in order to handle the left/south flux.
    460 !--             near vertical walls.
     420!> Initialization of flags to control the order of the advection scheme near solid walls and
     421!> non-cyclic inflow boundaries, where the order is sucessively degraded.
     422!--------------------------------------------------------------------------------------------------!
     423 SUBROUTINE ws_init_flags_momentum
     424
     425
     426    INTEGER(iwp) ::  i     !< index variable along x
     427    INTEGER(iwp) ::  j     !< index variable along y
     428    INTEGER(iwp) ::  k     !< index variable along z
     429    INTEGER(iwp) ::  k_mm  !< dummy index along z
     430    INTEGER(iwp) ::  k_pp  !< dummy index along z
     431    INTEGER(iwp) ::  k_ppp !< dummy index along z
     432
     433    LOGICAL      ::  flag_set !< steering variable for advection flags
     434
     435    ALLOCATE( advc_flags_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     436    advc_flags_m = 0
     437!
     438!-- Set advc_flags_m to steer the degradation of the advection scheme in advec_ws near
     439!-- topography, inflow- and outflow boundaries as well as bottom and top of model domain.
     440!-- advc_flags_m remains zero for all non-prognostic grid points.
     441!-- u-component
     442    DO  i = nxl, nxr
     443       DO  j = nys, nyn
     444          DO  k = nzb+1, nzt
     445!
     446!--          At first, set flags to WS1.
     447!--          Since fluxes are swapped in advec_ws.f90, this is necessary to
     448!--          in order to handle the left/south flux.
     449!--          near vertical walls.
     450             advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 )
     451             advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 )
     452!
     453!--          u component - x-direction
     454!--          WS1 (0), WS3 (1), WS5 (2)
     455             IF ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),1)  .OR.                                 &
     456                      ( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxlu  )  .OR.            &
     457                      ( ( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i == nxr   ) )                &
     458             THEN
    461459                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 )
     460             ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),1) .AND.                           &
     461                              BTEST(wall_flags_total_0(k,j,i+1),1) .OR.                            &
     462                        .NOT. BTEST(wall_flags_total_0(k,j,i-1),1) )                               &
     463                                              .OR.                                                 &
     464                      ( ( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i == nxr-1 )  .OR.            &
     465                      ( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i == nxlu+1) )                &
     466             THEN
     467                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 1 )
     468!
     469!--             Clear flag for WS1
     470                advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 )
     471             ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),1)  .AND.                                  &
     472                      BTEST(wall_flags_total_0(k,j,i+2),1)  .AND.                                  &
     473                      BTEST(wall_flags_total_0(k,j,i-1),1) )                                       &
     474             THEN
     475                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 2 )
     476!
     477!--             Clear flag for WS1
     478                advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 )
     479             ENDIF
     480!
     481!--          u component - y-direction
     482!--          WS1 (3), WS3 (4), WS5 (5)
     483             IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),1)  .OR.                                 &
     484                      ( ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j == nys   ) .OR.             &
     485                      ( ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j == nyn   ) )                &
     486             THEN
    462487                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 )
    463 !
    464 !--             u component - x-direction
    465 !--             WS1 (0), WS3 (1), WS5 (2)
    466                 IF ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),1)  .OR.          &
    467                          ( ( bc_dirichlet_l .OR. bc_radiation_l )              &
    468                            .AND. i <= nxlu  )    .OR.                          &
    469                          ( ( bc_dirichlet_r .OR. bc_radiation_r )              &
    470                            .AND. i == nxr   ) )                                &
    471                 THEN
    472                     advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 )
    473                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),1)  .AND.   &
    474                                  BTEST(wall_flags_total_0(k,j,i+1),1)  .OR.    &
    475                            .NOT. BTEST(wall_flags_total_0(k,j,i-1),1) )        &
    476                                                      .OR.                      &
    477                          ( ( bc_dirichlet_r .OR. bc_radiation_r )              &
    478                            .AND. i == nxr-1 )    .OR.                          &
    479                          ( ( bc_dirichlet_l .OR. bc_radiation_l )              &
    480                            .AND. i == nxlu+1) )                                &
    481                 THEN
    482                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 1 )
    483 !
    484 !--                Clear flag for WS1
    485                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 )
    486                 ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),1)  .AND.           &
    487                          BTEST(wall_flags_total_0(k,j,i+2),1)  .AND.           &
    488                          BTEST(wall_flags_total_0(k,j,i-1),1) )                &
    489                 THEN
    490                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 2 )
    491 !
    492 !--                Clear flag for WS1
    493                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 )
    494                 ENDIF
    495 !
    496 !--             u component - y-direction
    497 !--             WS1 (3), WS3 (4), WS5 (5)
    498                 IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),1)   .OR.         &
    499                          ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
    500                            .AND. j == nys   )    .OR.                          &
    501                          ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
    502                            .AND. j == nyn   ) )                                &
    503                 THEN
    504                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 )
    505                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),1)  .AND.   &
    506                                  BTEST(wall_flags_total_0(k,j+1,i),1)  .OR.    &
    507                            .NOT. BTEST(wall_flags_total_0(k,j-1,i),1) )        &
    508                                                      .OR.                      &
    509                          ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
    510                            .AND. j == nysv  )    .OR.                          &
    511                          ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
    512                            .AND. j == nyn-1 ) )                                &
    513                 THEN
    514                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 4 )
    515 !
    516 !--                Clear flag for WS1
    517                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 )
    518                 ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),1)  .AND.          &
    519                          BTEST(wall_flags_total_0(k,j+2,i),1)  .AND.          &
    520                          BTEST(wall_flags_total_0(k,j-1,i),1) )               &
    521                 THEN
    522                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 5 )
    523 !
    524 !--                Clear flag for WS1
    525                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 )
    526                 ENDIF
    527 !
    528 !--             u component - z-direction. Fluxes are calculated on w-grid
    529 !--             level. Boundary u-values at/within walls aren't used.
    530 !--             WS1 (6), WS3 (7), WS5 (8)
    531                 IF ( k == nzb+1 )  THEN
    532                    k_mm = nzb
    533                 ELSE
    534                    k_mm = k - 2
    535                 ENDIF
    536                 IF ( k > nzt-1 )  THEN
    537                    k_pp = nzt+1
    538                 ELSE
    539                    k_pp = k + 2
    540                 ENDIF
    541                 IF ( k > nzt-2 )  THEN
    542                    k_ppp = nzt+1
    543                 ELSE
    544                    k_ppp = k + 3
    545                 ENDIF
    546 
    547                 flag_set = .FALSE.
    548                 IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),1)       .AND.  &
    549                              BTEST(wall_flags_total_0(k,j,i),1)         .AND.  &
    550                              BTEST(wall_flags_total_0(k+1,j,i),1) )     .OR.   &
    551                      ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),1)      .AND.  &
    552                              BTEST(wall_flags_total_0(k+1,j,i),1)       .AND.  &
    553                              BTEST(wall_flags_total_0(k,j,i),1) )       .OR.   &
    554                      ( k == nzt .AND. symmetry_flag == 0 ) )                   &
    555                 THEN
    556                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 6 )
    557                    flag_set = .TRUE.
    558                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),1)    .OR. &
    559                            .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),1) ) .AND.&
    560                                  BTEST(wall_flags_total_0(k-1,j,i),1)     .AND.&
    561                                  BTEST(wall_flags_total_0(k,j,i),1)       .AND.&
    562                                  BTEST(wall_flags_total_0(k+1,j,i),1)     .AND.&
    563                                  BTEST(wall_flags_total_0(k_pp,j,i),1)    .AND.&
    564                            .NOT. flag_set                                  .OR.&
    565                          ( k == nzt - 1 .AND. symmetry_flag == 0 ) )           &
    566                 THEN
    567                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 )
    568                    flag_set = .TRUE.
    569                 ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),1)            .AND.&
    570                          BTEST(wall_flags_total_0(k-1,j,i),1)             .AND.&
    571                          BTEST(wall_flags_total_0(k,j,i),1)               .AND.&
    572                          BTEST(wall_flags_total_0(k+1,j,i),1)             .AND.&
    573                          BTEST(wall_flags_total_0(k_pp,j,i),1)            .AND.&
    574                          BTEST(wall_flags_total_0(k_ppp,j,i),1)           .AND.&
    575                          .NOT. flag_set )                                      &
    576                 THEN
    577                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 8 )
    578                 ENDIF
    579 
    580              ENDDO
     488             ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),1) .AND.                           &
     489                              BTEST(wall_flags_total_0(k,j+1,i),1) .OR.                            &
     490                        .NOT. BTEST(wall_flags_total_0(k,j-1,i),1) )                               &
     491                                              .OR.                                                 &
     492                      ( ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j == nysv  ) .OR.             &
     493                      ( ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j == nyn-1 ) )                &
     494             THEN
     495                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 4 )
     496!
     497!--             Clear flag for WS1
     498                advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 )
     499             ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),1)  .AND.                                  &
     500                      BTEST(wall_flags_total_0(k,j+2,i),1)  .AND.                                  &
     501                      BTEST(wall_flags_total_0(k,j-1,i),1) )                                       &
     502             THEN
     503                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 5 )
     504!
     505!--             Clear flag for WS1
     506                advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 )
     507             ENDIF
     508!
     509!--          u component - z-direction. Fluxes are calculated on w-grid level. Boundary u-values
     510!--          at/within walls aren't used.
     511!--          WS1 (6), WS3 (7), WS5 (8)
     512             IF ( k == nzb+1 )  THEN
     513                k_mm = nzb
     514             ELSE
     515                k_mm = k - 2
     516             ENDIF
     517             IF ( k > nzt-1 )  THEN
     518                k_pp = nzt+1
     519             ELSE
     520                k_pp = k + 2
     521             ENDIF
     522             IF ( k > nzt-2 )  THEN
     523                k_ppp = nzt+1
     524             ELSE
     525                k_ppp = k + 3
     526             ENDIF
     527
     528             flag_set = .FALSE.
     529             IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),1)       .AND.                         &
     530                          BTEST(wall_flags_total_0(k,j,i),1)         .AND.                         &
     531                          BTEST(wall_flags_total_0(k+1,j,i),1) )     .OR.                          &
     532                  ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),1)      .AND.                         &
     533                          BTEST(wall_flags_total_0(k+1,j,i),1)       .AND.                         &
     534                          BTEST(wall_flags_total_0(k,j,i),1) )       .OR.                          &
     535                  ( k == nzt .AND. symmetry_flag == 0 ) )                                          &
     536             THEN
     537                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 6 )
     538                flag_set = .TRUE.
     539             ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),1)    .OR.                        &
     540                        .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),1) ) .AND.                       &
     541                              BTEST(wall_flags_total_0(k-1,j,i),1)     .AND.                       &
     542                              BTEST(wall_flags_total_0(k,j,i),1)       .AND.                       &
     543                              BTEST(wall_flags_total_0(k+1,j,i),1)     .AND.                       &
     544                              BTEST(wall_flags_total_0(k_pp,j,i),1)    .AND.                       &
     545                        .NOT. flag_set                                 .OR.                        &
     546                      ( k == nzt - 1 .AND. symmetry_flag == 0 ) )                                  &
     547             THEN
     548                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 )
     549                flag_set = .TRUE.
     550             ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),1)            .AND.                       &
     551                      BTEST(wall_flags_total_0(k-1,j,i),1)             .AND.                       &
     552                      BTEST(wall_flags_total_0(k,j,i),1)               .AND.                       &
     553                      BTEST(wall_flags_total_0(k+1,j,i),1)             .AND.                       &
     554                      BTEST(wall_flags_total_0(k_pp,j,i),1)            .AND.                       &
     555                      BTEST(wall_flags_total_0(k_ppp,j,i),1)           .AND.                       &
     556                      .NOT. flag_set )                                                             &
     557             THEN
     558                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 8 )
     559             ENDIF
     560
    581561          ENDDO
    582562       ENDDO
    583 !
    584 !--    v-component
    585        DO  i = nxl, nxr
    586           DO  j = nys, nyn
    587              DO  k = nzb+1, nzt
    588 !
    589 !--             At first, set flags to WS1.
    590 !--             Since fluxes are swapped in advec_ws.f90, this is necessary to
    591 !--             in order to handle the left/south flux.
    592                 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9  )
     563    ENDDO
     564!
     565!-- v-component
     566    DO  i = nxl, nxr
     567       DO  j = nys, nyn
     568          DO  k = nzb+1, nzt
     569!
     570!--          At first, set flags to WS1.
     571!--          Since fluxes are swapped in advec_ws.f90, this is necessary to in order to handle the
     572!--          left/south flux.
     573             advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9  )
     574             advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 )
     575!
     576!--          v component - x-direction
     577!--          WS1 (9), WS3 (10), WS5 (11)
     578             IF ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),2)                       .OR.            &
     579                      ( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i == nxl   )  .OR.            &
     580                      ( ( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i == nxr   ) )                &
     581            THEN
     582                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 )
     583!
     584!--          WS3
     585             ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),2)   .AND.                         &
     586                              BTEST(wall_flags_total_0(k,j,i+1),2) ) .OR.                          &
     587                        .NOT. BTEST(wall_flags_total_0(k,j,i-1),2)                                 &
     588                                              .OR.                                                 &
     589                      ( ( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i == nxr-1 )  .OR.            &
     590                      ( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i == nxlu  ) )                &
     591             THEN
     592                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 10 )
     593!
     594!--             Clear flag for WS1
     595                advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 )
     596             ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),2)  .AND.                                  &
     597                      BTEST(wall_flags_total_0(k,j,i+2),2)  .AND.                                  &
     598                      BTEST(wall_flags_total_0(k,j,i-1),2) )                                       &
     599             THEN
     600                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 11 )
     601!
     602!--             Clear flag for WS1
     603                advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 )
     604             ENDIF
     605!
     606!--          v component - y-direction
     607!--          WS1 (12), WS3 (13), WS5 (14)
     608             IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),2)                      .OR.             &
     609                      ( ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j <= nysv  ) .OR.             &
     610                      ( ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j == nyn   ) )                &
     611             THEN
    593612                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 )
    594 !
    595 !--             v component - x-direction
    596 !--             WS1 (9), WS3 (10), WS5 (11)
    597                 IF ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),2)  .OR.          &
    598                          ( ( bc_dirichlet_l .OR. bc_radiation_l )              &
    599                            .AND. i == nxl   )    .OR.                          &
    600                          ( ( bc_dirichlet_r .OR. bc_radiation_r )              &
    601                            .AND. i == nxr   ) )                                &
    602                THEN
    603                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 )
    604 !
    605 !--             WS3
    606                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),2)   .AND.  &
    607                                  BTEST(wall_flags_total_0(k,j,i+1),2) ) .OR.   &
    608                            .NOT. BTEST(wall_flags_total_0(k,j,i-1),2)          &
    609                                                  .OR.                          &
    610                          ( ( bc_dirichlet_r .OR. bc_radiation_r )              &
    611                            .AND. i == nxr-1 )    .OR.                          &
    612                          ( ( bc_dirichlet_l .OR. bc_radiation_l )              &
    613                            .AND. i == nxlu  ) )                                &
    614                 THEN
    615                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 10 )
    616 !
    617 !--                Clear flag for WS1
    618                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 )
    619                 ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),2)  .AND.           &
    620                          BTEST(wall_flags_total_0(k,j,i+2),2)  .AND.           &
    621                          BTEST(wall_flags_total_0(k,j,i-1),2) )                &
    622                 THEN
    623                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 11 )
    624 !
    625 !--                Clear flag for WS1
    626                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 )
    627                 ENDIF
    628 !
    629 !--             v component - y-direction
    630 !--             WS1 (12), WS3 (13), WS5 (14)
    631                 IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),2) .OR.           &
    632                          ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
    633                            .AND. j <= nysv  )    .OR.                          &
    634                          ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
    635                            .AND. j == nyn   ) )                                &
    636                 THEN
    637                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 )
    638                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),2)  .AND.   &
    639                                  BTEST(wall_flags_total_0(k,j+1,i),2)  .OR.    &
    640                            .NOT. BTEST(wall_flags_total_0(k,j-1,i),2) )        &
    641                                                      .OR.                      &
    642                          ( (  bc_dirichlet_s .OR. bc_radiation_s )             &
    643                            .AND. j == nysv+1)    .OR.                          &
    644                          ( (  bc_dirichlet_n .OR. bc_radiation_n )             &
    645                            .AND. j == nyn-1 ) )                                &
    646                 THEN
    647                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 13 )
    648 !
    649 !--                Clear flag for WS1
    650                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 )
    651                 ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),2)  .AND.           &
    652                          BTEST(wall_flags_total_0(k,j+2,i),2)  .AND.           &
    653                          BTEST(wall_flags_total_0(k,j-1,i),2) )                &
    654                 THEN
    655                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 14 )
    656 !
    657 !--                Clear flag for WS1
    658                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 )
    659                 ENDIF
    660 !
    661 !--             v component - z-direction. Fluxes are calculated on w-grid
    662 !--             level. Boundary v-values at/within walls aren't used.
    663 !--             WS1 (15), WS3 (16), WS5 (17)
    664                 IF ( k == nzb+1 )  THEN
    665                    k_mm = nzb
    666                 ELSE
    667                    k_mm = k - 2
    668                 ENDIF
    669                 IF ( k > nzt-1 )  THEN
    670                    k_pp = nzt+1
    671                 ELSE
    672                    k_pp = k + 2
    673                 ENDIF
    674                 IF ( k > nzt-2 )  THEN
    675                    k_ppp = nzt+1
    676                 ELSE
    677                    k_ppp = k + 3
    678                 ENDIF
    679 
    680                 flag_set = .FALSE.
    681                 IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),2)         .AND.&
    682                              BTEST(wall_flags_total_0(k,j,i),2)           .AND.&
    683                              BTEST(wall_flags_total_0(k+1,j,i),2) )       .OR. &
    684                      ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),2)        .AND.&
    685                              BTEST(wall_flags_total_0(k+1,j,i),2)         .AND.&
    686                              BTEST(wall_flags_total_0(k,j,i),2) )         .OR. &
    687                      ( k == nzt .AND. symmetry_flag == 0 ) )                   &
    688                 THEN
    689                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 15 )
    690                    flag_set = .TRUE.
    691                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),2)    .OR. &
    692                            .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),2) ) .AND.&
    693                                  BTEST(wall_flags_total_0(k-1,j,i),2)     .AND.&
    694                                  BTEST(wall_flags_total_0(k,j,i),2)       .AND.&
    695                                  BTEST(wall_flags_total_0(k+1,j,i),2)     .AND.&
    696                                  BTEST(wall_flags_total_0(k_pp,j,i),2)    .AND.&
    697                            .NOT. flag_set                                  .OR.&
    698                          ( k == nzt - 1 .AND. symmetry_flag == 0 ) )           &
    699                 THEN
    700                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 )
    701                    flag_set = .TRUE.
    702                 ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),2)            .AND.&
    703                          BTEST(wall_flags_total_0(k-1,j,i),2)             .AND.&
    704                          BTEST(wall_flags_total_0(k,j,i),2)               .AND.&
    705                          BTEST(wall_flags_total_0(k+1,j,i),2)             .AND.&
    706                          BTEST(wall_flags_total_0(k_pp,j,i),2)            .AND.&
    707                          BTEST(wall_flags_total_0(k_ppp,j,i),2)           .AND.&
    708                          .NOT. flag_set )                                      &
    709                 THEN
    710                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 17 )
    711                 ENDIF
    712 
    713              ENDDO
     613             ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),2)  .AND.                          &
     614                              BTEST(wall_flags_total_0(k,j+1,i),2)  .OR.                           &
     615                        .NOT. BTEST(wall_flags_total_0(k,j-1,i),2) )                               &
     616                                              .OR.                                                 &
     617                      ( (  bc_dirichlet_s .OR. bc_radiation_s ) .AND. j == nysv+1) .OR.            &
     618                      ( (  bc_dirichlet_n .OR. bc_radiation_n ) .AND. j == nyn-1 ) )               &
     619             THEN
     620                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 13 )
     621!
     622!--             Clear flag for WS1
     623                advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 )
     624             ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),2)  .AND.                                  &
     625                      BTEST(wall_flags_total_0(k,j+2,i),2)  .AND.                                  &
     626                      BTEST(wall_flags_total_0(k,j-1,i),2) )                                       &
     627             THEN
     628                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 14 )
     629!
     630!--             Clear flag for WS1
     631                advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 )
     632             ENDIF
     633!
     634!--          v component - z-direction. Fluxes are calculated on w-grid level. Boundary v-values
     635!--          at/within walls aren't used.
     636!--          WS1 (15), WS3 (16), WS5 (17)
     637             IF ( k == nzb+1 )  THEN
     638                k_mm = nzb
     639             ELSE
     640                k_mm = k - 2
     641             ENDIF
     642             IF ( k > nzt-1 )  THEN
     643                k_pp = nzt+1
     644             ELSE
     645                k_pp = k + 2
     646             ENDIF
     647             IF ( k > nzt-2 )  THEN
     648                k_ppp = nzt+1
     649             ELSE
     650                k_ppp = k + 3
     651             ENDIF
     652
     653             flag_set = .FALSE.
     654             IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),2)         .AND.                       &
     655                          BTEST(wall_flags_total_0(k,j,i),2)           .AND.                       &
     656                          BTEST(wall_flags_total_0(k+1,j,i),2) )       .OR.                        &
     657                  ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),2)        .AND.                       &
     658                          BTEST(wall_flags_total_0(k+1,j,i),2)         .AND.                       &
     659                          BTEST(wall_flags_total_0(k,j,i),2) )         .OR.                        &
     660                  ( k == nzt .AND. symmetry_flag == 0 ) )                                          &
     661             THEN
     662                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 15 )
     663                flag_set = .TRUE.
     664             ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),2)    .OR.                        &
     665                        .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),2) ) .AND.                       &
     666                              BTEST(wall_flags_total_0(k-1,j,i),2)     .AND.                       &
     667                              BTEST(wall_flags_total_0(k,j,i),2)       .AND.                       &
     668                              BTEST(wall_flags_total_0(k+1,j,i),2)     .AND.                       &
     669                              BTEST(wall_flags_total_0(k_pp,j,i),2)    .AND.                       &
     670                        .NOT. flag_set                                  .OR.                       &
     671                      ( k == nzt - 1 .AND. symmetry_flag == 0 ) )                                  &
     672             THEN
     673                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 )
     674                flag_set = .TRUE.
     675             ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),2)            .AND.                       &
     676                      BTEST(wall_flags_total_0(k-1,j,i),2)             .AND.                       &
     677                      BTEST(wall_flags_total_0(k,j,i),2)               .AND.                       &
     678                      BTEST(wall_flags_total_0(k+1,j,i),2)             .AND.                       &
     679                      BTEST(wall_flags_total_0(k_pp,j,i),2)            .AND.                       &
     680                      BTEST(wall_flags_total_0(k_ppp,j,i),2)           .AND.                       &
     681                      .NOT. flag_set )                                                             &
     682             THEN
     683                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 17 )
     684             ENDIF
     685
    714686          ENDDO
    715687       ENDDO
    716 !
    717 !--    w - component
    718        DO  i = nxl, nxr
    719           DO  j = nys, nyn
    720              DO  k = nzb+1, nzt
    721 !
    722 !--             At first, set flags to WS1.
    723 !--             Since fluxes are swapped in advec_ws.f90, this is necessary to
    724 !--             in order to handle the left/south flux.
     688    ENDDO
     689!
     690!-- w - component
     691    DO  i = nxl, nxr
     692       DO  j = nys, nyn
     693          DO  k = nzb+1, nzt
     694!
     695!--          At first, set flags to WS1.
     696!--          Since fluxes are swapped in advec_ws.f90, this is necessary to in order to handle the
     697!--          left/south flux.
     698             advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 )
     699             advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 )
     700!
     701!--          w component - x-direction
     702!--          WS1 (18), WS3 (19), WS5 (20)
     703             IF ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),3)                        .OR.           &
     704                      ( (  bc_dirichlet_l .OR. bc_radiation_l ) .AND. i == nxl   )  .OR.           &
     705                      ( (  bc_dirichlet_r .OR. bc_radiation_r ) .AND. i == nxr   ) )               &
     706             THEN
    725707                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 )
     708             ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),3)  .AND.                          &
     709                              BTEST(wall_flags_total_0(k,j,i+1),3)  .OR.                           &
     710                        .NOT. BTEST(wall_flags_total_0(k,j,i-1),3) )                               &
     711                                              .OR.                                                 &
     712                      ( ( bc_dirichlet_r .OR. bc_radiation_r )  .AND. i == nxr-1 )  .OR.           &
     713                      ( ( bc_dirichlet_l .OR.  bc_radiation_l ) .AND. i == nxlu  ) )               &
     714             THEN
     715                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 19 )
     716!
     717!--             Clear flag for WS1
     718                advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 )
     719             ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),3)  .AND.                                  &
     720                      BTEST(wall_flags_total_0(k,j,i+2),3)  .AND.                                  &
     721                      BTEST(wall_flags_total_0(k,j,i-1),3) )                                       &
     722             THEN
     723                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i),20 )
     724!
     725!--             Clear flag for WS1
     726                advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 )
     727             ENDIF
     728!
     729!--          w component - y-direction
     730!--          WS1 (21), WS3 (22), WS5 (23)
     731             IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),3)                       .OR.            &
     732                      ( ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j == nys   )  .OR.            &
     733                      ( ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j == nyn   ) )                &
     734             THEN
    726735                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 )
    727 !
    728 !--             w component - x-direction
    729 !--             WS1 (18), WS3 (19), WS5 (20)
    730                 IF ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),3) .OR.           &
    731                          ( (  bc_dirichlet_l .OR. bc_radiation_l )             &
    732                            .AND. i == nxl   )    .OR.                          &
    733                          ( (  bc_dirichlet_r .OR. bc_radiation_r )             &
    734                            .AND. i == nxr   ) )                                &
    735                 THEN
    736                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 )
    737                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),3)  .AND.   &
    738                                  BTEST(wall_flags_total_0(k,j,i+1),3)  .OR.    &
    739                            .NOT. BTEST(wall_flags_total_0(k,j,i-1),3) )        &
    740                                                      .OR.                      &
    741                          ( ( bc_dirichlet_r .OR. bc_radiation_r )              &
    742                            .AND. i == nxr-1 )    .OR.                          &
    743                          ( ( bc_dirichlet_l .OR.  bc_radiation_l )             &
    744                            .AND. i == nxlu  ) )                                &
    745                 THEN
    746                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 19 )
    747 !
    748 !--                Clear flag for WS1
    749                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 )
    750                 ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),3)  .AND.           &
    751                          BTEST(wall_flags_total_0(k,j,i+2),3)  .AND.           &
    752                          BTEST(wall_flags_total_0(k,j,i-1),3) )                &
    753                 THEN
    754                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i),20 )
    755 !
    756 !--                Clear flag for WS1
    757                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 )
    758                 ENDIF
    759 !
    760 !--             w component - y-direction
    761 !--             WS1 (21), WS3 (22), WS5 (23)
    762                 IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),3) .OR.           &
    763                          ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
    764                            .AND. j == nys   )    .OR.                          &
    765                          ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
    766                            .AND. j == nyn   ) )                                &
    767                 THEN
    768                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 )
    769                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),3)  .AND.   &
    770                                  BTEST(wall_flags_total_0(k,j+1,i),3)  .OR.    &
    771                            .NOT. BTEST(wall_flags_total_0(k,j-1,i),3) )        &
    772                                                      .OR.                      &
    773                          ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
    774                            .AND. j == nysv  )    .OR.                          &
    775                          ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
    776                            .AND. j == nyn-1 ) )                                &
    777                 THEN
    778                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 22 )
    779 !
    780 !--                Clear flag for WS1
    781                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 )
    782                 ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),3)  .AND.           &
    783                          BTEST(wall_flags_total_0(k,j+2,i),3)  .AND.           &
    784                          BTEST(wall_flags_total_0(k,j-1,i),3) )                &
    785                 THEN
    786                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 23 )
    787 !
    788 !--                Clear flag for WS1
    789                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 )
    790                 ENDIF
    791 !
    792 !--             w component - z-direction. Fluxes are calculated on scalar grid
    793 !--             level. Boundary w-values at walls are used. Flux at k=i is
    794 !--             defined at scalar position k=i+1 with i being an integer.
    795 !--             WS1 (24), WS3 (25), WS5 (26)
    796                 IF ( k == nzb+1 )  THEN
    797                    k_mm = nzb
    798                 ELSE
    799                    k_mm = k - 2
    800                 ENDIF
    801                 IF ( k > nzt-1 )  THEN
    802                    k_pp = nzt+1
    803                 ELSE
    804                    k_pp = k + 2
    805                 ENDIF
    806                 IF ( k > nzt-2 )  THEN
    807                    k_ppp = nzt+1
    808                 ELSE
    809                    k_ppp = k + 3
    810                 ENDIF
    811 
    812                 flag_set = .FALSE.
    813                 IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i),3)          .AND. &
    814                              BTEST(wall_flags_total_0(k+1,j,i),3) )      .OR.  &
    815                      ( .NOT. BTEST(wall_flags_total_0(k+1,j,i),3)        .AND. &
    816                              BTEST(wall_flags_total_0(k,j,i),3) )        .OR.  &
    817                      k == nzt -1 )                                             &
    818                 THEN
    819 !
    820 !--                Please note, at k == nzb_w_inner(j,i) a flag is explicitly
    821 !--                set, although this is not a prognostic level. However,
    822 !--                contrary to the advection of u,v and s this is necessary
    823 !--                because flux_t(nzb_w_inner(j,i)) is used for the tendency
    824 !--                at k == nzb_w_inner(j,i)+1.
    825                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 24 )
    826                    flag_set = .TRUE.
    827                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),3)    .AND. &
    828                                  BTEST(wall_flags_total_0(k,j,i),3)      .AND. &
    829                                  BTEST(wall_flags_total_0(k+1,j,i),3)    .AND. &
    830                                  BTEST(wall_flags_total_0(k_pp,j,i),3) ) .OR.  &
    831                          ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),3)   .AND. &
    832                                  BTEST(wall_flags_total_0(k+1,j,i),3)    .AND. &
    833                                  BTEST(wall_flags_total_0(k,j,i),3)      .AND. &
    834                                  BTEST(wall_flags_total_0(k-1,j,i),3) )  .AND. &
    835                            .NOT. flag_set                          .OR.        &
    836                          k == nzt - 2 )                                        &
    837                 THEN
    838                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 25 )
    839                    flag_set = .TRUE.
    840                 ELSEIF ( BTEST(wall_flags_total_0(k-1,j,i),3)            .AND. &
    841                          BTEST(wall_flags_total_0(k,j,i),3)              .AND. &
    842                          BTEST(wall_flags_total_0(k+1,j,i),3)            .AND. &
    843                          BTEST(wall_flags_total_0(k_pp,j,i),3)           .AND. &
    844                          .NOT. flag_set )                                      &
    845                 THEN
    846                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 26 )
    847                 ENDIF
    848 
    849              ENDDO
     736             ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),3)  .AND.                          &
     737                              BTEST(wall_flags_total_0(k,j+1,i),3)  .OR.                           &
     738                        .NOT. BTEST(wall_flags_total_0(k,j-1,i),3) )                               &
     739                                              .OR.                                                 &
     740                      ( ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j == nysv  )  .OR.            &
     741                      ( ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j == nyn-1 ) )                &
     742             THEN
     743                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 22 )
     744!
     745!--             Clear flag for WS1
     746                advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 )
     747             ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),3)  .AND.                                  &
     748                      BTEST(wall_flags_total_0(k,j+2,i),3)  .AND.                                  &
     749                      BTEST(wall_flags_total_0(k,j-1,i),3) )                                       &
     750             THEN
     751                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 23 )
     752!
     753!--             Clear flag for WS1
     754                advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 )
     755             ENDIF
     756!
     757!--          w component - z-direction. Fluxes are calculated on scalar grid level. Boundary
     758!--          w-values at walls are used. Flux at k=i is defined at scalar position k=i+1 with i
     759!--          being an integer.
     760!--          WS1 (24), WS3 (25), WS5 (26)
     761             IF ( k == nzb+1 )  THEN
     762                k_mm = nzb
     763             ELSE
     764                k_mm = k - 2
     765             ENDIF
     766             IF ( k > nzt-1 )  THEN
     767                k_pp = nzt+1
     768             ELSE
     769                k_pp = k + 2
     770             ENDIF
     771             IF ( k > nzt-2 )  THEN
     772                k_ppp = nzt+1
     773             ELSE
     774                k_ppp = k + 3
     775             ENDIF
     776
     777             flag_set = .FALSE.
     778             IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i),3)          .AND.                        &
     779                          BTEST(wall_flags_total_0(k+1,j,i),3) )      .OR.                         &
     780                  ( .NOT. BTEST(wall_flags_total_0(k+1,j,i),3)        .AND.                        &
     781                          BTEST(wall_flags_total_0(k,j,i),3) )        .OR.                         &
     782                  k == nzt -1 )                                                                    &
     783             THEN
     784!
     785!--             Please note, at k == nzb_w_inner(j,i) a flag is explicitly set, although this is not
     786!--             a prognostic level. However, contrary to the advection of u,v and s this is
     787!--             necessary because flux_t(nzb_w_inner(j,i)) is used for the tendency at k ==
     788!--             0nzb_w_inner(j,i)+1.
     789                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 24 )
     790                flag_set = .TRUE.
     791             ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),3)    .AND.                        &
     792                              BTEST(wall_flags_total_0(k,j,i),3)      .AND.                        &
     793                              BTEST(wall_flags_total_0(k+1,j,i),3)    .AND.                        &
     794                              BTEST(wall_flags_total_0(k_pp,j,i),3) ) .OR.                         &
     795                      ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),3)   .AND.                        &
     796                              BTEST(wall_flags_total_0(k+1,j,i),3)    .AND.                        &
     797                              BTEST(wall_flags_total_0(k,j,i),3)      .AND.                        &
     798                              BTEST(wall_flags_total_0(k-1,j,i),3) )  .AND.                        &
     799                        .NOT. flag_set                                .OR.                         &
     800                      k == nzt - 2 )                                                               &
     801             THEN
     802                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 25 )
     803                flag_set = .TRUE.
     804             ELSEIF ( BTEST(wall_flags_total_0(k-1,j,i),3)            .AND.                        &
     805                      BTEST(wall_flags_total_0(k,j,i),3)              .AND.                        &
     806                      BTEST(wall_flags_total_0(k+1,j,i),3)            .AND.                        &
     807                      BTEST(wall_flags_total_0(k_pp,j,i),3)           .AND.                        &
     808                      .NOT. flag_set )                                                             &
     809             THEN
     810                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 26 )
     811             ENDIF
     812
    850813          ENDDO
    851814       ENDDO
    852 !
    853 !--    Exchange ghost points for advection flags
    854        CALL exchange_horiz_int( advc_flags_m, nys, nyn, nxl, nxr, nzt, nbgp )
    855 !
    856 !--    Set boundary flags at inflow and outflow boundary in case of
    857 !--    non-cyclic boundary conditions.
    858        IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
    859           advc_flags_m(:,:,nxl-1) = advc_flags_m(:,:,nxl)
    860        ENDIF
    861 
    862        IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
    863           advc_flags_m(:,:,nxr+1) = advc_flags_m(:,:,nxr)
    864        ENDIF
    865 
    866        IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
    867           advc_flags_m(:,nyn+1,:) = advc_flags_m(:,nyn,:)
    868        ENDIF
    869 
    870        IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
    871           advc_flags_m(:,nys-1,:) = advc_flags_m(:,nys,:)
    872        ENDIF
    873 
    874     END SUBROUTINE ws_init_flags_momentum
    875 
    876 
    877 !------------------------------------------------------------------------------!
     815    ENDDO
     816!
     817!-- Exchange ghost points for advection flags
     818    CALL exchange_horiz_int( advc_flags_m, nys, nyn, nxl, nxr, nzt, nbgp )
     819!
     820!-- Set boundary flags at inflow and outflow boundary in case of
     821!-- non-cyclic boundary conditions.
     822    IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
     823       advc_flags_m(:,:,nxl-1) = advc_flags_m(:,:,nxl)
     824    ENDIF
     825
     826    IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
     827       advc_flags_m(:,:,nxr+1) = advc_flags_m(:,:,nxr)
     828    ENDIF
     829
     830    IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
     831       advc_flags_m(:,nyn+1,:) = advc_flags_m(:,nyn,:)
     832    ENDIF
     833
     834    IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
     835       advc_flags_m(:,nys-1,:) = advc_flags_m(:,nys,:)
     836    ENDIF
     837
     838 END SUBROUTINE ws_init_flags_momentum
     839
     840
     841!--------------------------------------------------------------------------------------------------!
    878842! Description:
    879843! ------------
    880 !> Initialization of flags to control the order of the advection scheme near
    881 !> solid walls and non-cyclic inflow boundaries, where the order is sucessively
    882 !> degraded.
    883 !------------------------------------------------------------------------------!
    884     SUBROUTINE ws_init_flags_scalar( non_cyclic_l, non_cyclic_n, non_cyclic_r, &
    885                                      non_cyclic_s, advc_flag, extensive_degrad )
    886 
    887 
    888        INTEGER(iwp) ::  i     !< index variable along x
    889        INTEGER(iwp) ::  j     !< index variable along y
    890        INTEGER(iwp) ::  k     !< index variable along z
    891        INTEGER(iwp) ::  k_mm  !< dummy index along z
    892        INTEGER(iwp) ::  k_pp  !< dummy index along z
    893        INTEGER(iwp) ::  k_ppp !< dummy index along z
    894 
    895        INTEGER(iwp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::&
    896                                                   advc_flag !< flag array to control order of scalar advection
    897 
    898        LOGICAL ::  flag_set     !< steering variable for advection flags
    899        LOGICAL ::  non_cyclic_l !< flag that indicates non-cyclic boundary on the left
    900        LOGICAL ::  non_cyclic_n !< flag that indicates non-cyclic boundary on the north
    901        LOGICAL ::  non_cyclic_r !< flag that indicates non-cyclic boundary on the right
    902        LOGICAL ::  non_cyclic_s !< flag that indicates non-cyclic boundary on the south
    903 
    904        LOGICAL, OPTIONAL ::  extensive_degrad !< flag indicating that extensive degradation is required, e.g. for
    905                                               !< passive scalars nearby topography along the horizontal directions,
    906                                               !< as no monotonic limiter can be applied there
    907 !
    908 !--    Set flags to steer the degradation of the advection scheme in advec_ws
    909 !--    near topography, inflow- and outflow boundaries as well as bottom and
    910 !--    top of model domain. advc_flags_m remains zero for all non-prognostic
    911 !--    grid points.
    912        DO  i = nxl, nxr
    913           DO  j = nys, nyn
    914              DO  k = nzb+1, nzt
    915                 IF ( .NOT.  BTEST(wall_flags_total_0(k,j,i),0) )  CYCLE
    916 !
    917 !--             scalar - x-direction
    918 !--             WS1 (0), WS3 (1), WS5 (2)
    919                 IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),0)      .OR.    &
    920                        .NOT. BTEST(wall_flags_total_0(k,j,i+2),0)      .OR.    &
    921                        .NOT. BTEST(wall_flags_total_0(k,j,i-1),0) )    .OR.    &
    922                        ( non_cyclic_l  .AND.  i == 0  )                .OR.    &
    923                        ( non_cyclic_r  .AND.  i == nx ) )  THEN
    924                  advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
    925                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+3),0)  .AND.   &
    926                                  BTEST(wall_flags_total_0(k,j,i+1),0)  .AND.   &
    927                                  BTEST(wall_flags_total_0(k,j,i+2),0)  .AND.   &
    928                                  BTEST(wall_flags_total_0(k,j,i-1),0)          &
    929                          )                       .OR.                          &
    930                          ( .NOT. BTEST(wall_flags_total_0(k,j,i-2),0)  .AND.   &
    931                                  BTEST(wall_flags_total_0(k,j,i+1),0)  .AND.   &
    932                                  BTEST(wall_flags_total_0(k,j,i+2),0)  .AND.   &
    933                                  BTEST(wall_flags_total_0(k,j,i-1),0)          &
    934                          )                                                     &
    935                                                  .OR.                          &
    936                          ( non_cyclic_r  .AND.  i == nx-1 )  .OR.              &
    937                          ( non_cyclic_l  .AND.  i == 1    ) )  THEN
    938                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
    939                 ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),0)           .AND.  &
    940                          BTEST(wall_flags_total_0(k,j,i+2),0)           .AND.  &
    941                          BTEST(wall_flags_total_0(k,j,i+3),0)           .AND.  &
    942                          BTEST(wall_flags_total_0(k,j,i-1),0)           .AND.  &
    943                          BTEST(wall_flags_total_0(k,j,i-2),0) )                &
    944                 THEN
    945                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 )
    946                 ENDIF
    947 !
    948 !--             scalar - y-direction
    949 !--             WS1 (3), WS3 (4), WS5 (5)
    950                 IF ( ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),0)        .OR.  &
    951                        .NOT. BTEST(wall_flags_total_0(k,j+2,i),0)        .OR.  &
    952                        .NOT. BTEST(wall_flags_total_0(k,j-1,i),0))       .OR.  &
    953                      ( non_cyclic_s  .AND.  j == 0  )                    .OR.  &
    954                      ( non_cyclic_n  .AND.  j == ny ) )  THEN
    955                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
    956 !
    957 !--             WS3
    958                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+3,i),0)    .AND. &
    959                                  BTEST(wall_flags_total_0(k,j+1,i),0)    .AND. &
    960                                  BTEST(wall_flags_total_0(k,j+2,i),0)    .AND. &
    961                                  BTEST(wall_flags_total_0(k,j-1,i),0)          &
    962                          )                       .OR.                          &
    963                          ( .NOT. BTEST(wall_flags_total_0(k,j-2,i),0)    .AND. &
    964                                  BTEST(wall_flags_total_0(k,j+1,i),0)    .AND. &
    965                                  BTEST(wall_flags_total_0(k,j+2,i),0)    .AND. &
    966                                  BTEST(wall_flags_total_0(k,j-1,i),0)          &
    967                          )                                                     &
    968                                                  .OR.                          &
    969                          ( non_cyclic_s  .AND.  j == 1    )  .OR.              &
    970                          ( non_cyclic_n  .AND.  j == ny-1 ) )  THEN
    971                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
    972 !
    973 !--             WS5
    974                 ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),0)           .AND.  &
    975                          BTEST(wall_flags_total_0(k,j+2,i),0)           .AND.  &
    976                          BTEST(wall_flags_total_0(k,j+3,i),0)           .AND.  &
    977                          BTEST(wall_flags_total_0(k,j-1,i),0)           .AND.  &
    978                          BTEST(wall_flags_total_0(k,j-2,i),0) )                &
    979                 THEN
    980                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 )
    981                 ENDIF
    982 !
    983 !--             Near topography, set horizontal advection scheme to 1st order
    984 !--             for passive scalars, even if only one direction may be
    985 !--             blocked by topography. These locations will be identified
    986 !--             by wall_flags_total_0 bit 31. Note, since several modules define
    987 !--             advection flags but may apply different scalar boundary
    988 !--             conditions, bit 31 is temporarily stored on advc_flags.
    989 !--             Moreover, note that this extended degradtion for passive
    990 !--             scalars is not required for the vertical direction as there
    991 !--             the monotonic limiter can be applied.
    992                 IF ( PRESENT( extensive_degrad ) )  THEN
    993                    IF ( extensive_degrad )  THEN
    994 !
    995 !--                   At all grid points that are within a three-grid point
    996 !--                   range to topography, set 1st-order scheme.
    997                       IF( BTEST( advc_flag(k,j,i), 31 ) )  THEN
    998 !
    999 !--                      Clear flags that might indicate higher-order
    1000 !--                      advection along x- and y-direction.
     844!> Initialization of flags to control the order of the advection scheme near solid walls and
     845!> non-cyclic inflow boundaries, where the order is sucessively degraded.
     846!--------------------------------------------------------------------------------------------------!
     847 SUBROUTINE ws_init_flags_scalar( non_cyclic_l, non_cyclic_n, non_cyclic_r, non_cyclic_s,          &
     848                                  advc_flag, extensive_degrad )
     849
     850
     851    INTEGER(iwp) ::  i     !< index variable along x
     852    INTEGER(iwp) ::  j     !< index variable along y
     853    INTEGER(iwp) ::  k     !< index variable along z
     854    INTEGER(iwp) ::  k_mm  !< dummy index along z
     855    INTEGER(iwp) ::  k_pp  !< dummy index along z
     856    INTEGER(iwp) ::  k_ppp !< dummy index along z
     857
     858    INTEGER(iwp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::                       &
     859                                           advc_flag !< flag array to control order of scalar advection
     860
     861    LOGICAL ::  flag_set     !< steering variable for advection flags
     862    LOGICAL ::  non_cyclic_l !< flag that indicates non-cyclic boundary on the left
     863    LOGICAL ::  non_cyclic_n !< flag that indicates non-cyclic boundary on the north
     864    LOGICAL ::  non_cyclic_r !< flag that indicates non-cyclic boundary on the right
     865    LOGICAL ::  non_cyclic_s !< flag that indicates non-cyclic boundary on the south
     866
     867    LOGICAL, OPTIONAL ::  extensive_degrad !< flag indicating that extensive degradation is required, e.g. for
     868                                           !< passive scalars nearby topography along the horizontal directions,
     869                                           !< as no monotonic limiter can be applied there
     870!
     871!-- Set flags to steer the degradation of the advection scheme in advec_ws near topography, inflow-
     872!-- and outflow boundaries as well as bottom and top of model domain. advc_flags_m remains zero for
     873!-- all non-prognostic grid points.
     874    DO  i = nxl, nxr
     875       DO  j = nys, nyn
     876          DO  k = nzb+1, nzt
     877             IF ( .NOT.  BTEST(wall_flags_total_0(k,j,i),0) )  CYCLE
     878!
     879!--          scalar - x-direction
     880!--          WS1 (0), WS3 (1), WS5 (2)
     881             IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),0)      .OR.                           &
     882                    .NOT. BTEST(wall_flags_total_0(k,j,i+2),0)      .OR.                           &
     883                    .NOT. BTEST(wall_flags_total_0(k,j,i-1),0) )    .OR.                           &
     884                    ( non_cyclic_l  .AND.  i == 0  )                .OR.                           &
     885                    ( non_cyclic_r  .AND.  i == nx ) )                                             &
     886             THEN
     887                advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
     888             ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+3),0)  .AND.                          &
     889                              BTEST(wall_flags_total_0(k,j,i+1),0)  .AND.                          &
     890                              BTEST(wall_flags_total_0(k,j,i+2),0)  .AND.                          &
     891                              BTEST(wall_flags_total_0(k,j,i-1),0)                                 &
     892                      )                       .OR.                                                 &
     893                      ( .NOT. BTEST(wall_flags_total_0(k,j,i-2),0)  .AND.                          &
     894                              BTEST(wall_flags_total_0(k,j,i+1),0)  .AND.                          &
     895                              BTEST(wall_flags_total_0(k,j,i+2),0)  .AND.                          &
     896                              BTEST(wall_flags_total_0(k,j,i-1),0)                                 &
     897                      )                       .OR.                                                 &
     898                      ( non_cyclic_r  .AND.  i == nx-1 )            .OR.                           &
     899                      ( non_cyclic_l  .AND.  i == 1    ) )                                         &
     900             THEN
     901                advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
     902             ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),0)           .AND.                         &
     903                      BTEST(wall_flags_total_0(k,j,i+2),0)           .AND.                         &
     904                      BTEST(wall_flags_total_0(k,j,i+3),0)           .AND.                         &
     905                      BTEST(wall_flags_total_0(k,j,i-1),0)           .AND.                         &
     906                      BTEST(wall_flags_total_0(k,j,i-2),0) )                                       &
     907             THEN
     908                advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 )
     909             ENDIF
     910!
     911!--          scalar - y-direction
     912!--          WS1 (3), WS3 (4), WS5 (5)
     913             IF ( ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),0)        .OR.                         &
     914                    .NOT. BTEST(wall_flags_total_0(k,j+2,i),0)        .OR.                         &
     915                    .NOT. BTEST(wall_flags_total_0(k,j-1,i),0))       .OR.                         &
     916                  ( non_cyclic_s  .AND.  j == 0  )                    .OR.                         &
     917                  ( non_cyclic_n  .AND.  j == ny ) )                                               &
     918             THEN
     919                advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
     920!
     921!--          WS3
     922             ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+3,i),0)    .AND.                        &
     923                              BTEST(wall_flags_total_0(k,j+1,i),0)    .AND.                        &
     924                              BTEST(wall_flags_total_0(k,j+2,i),0)    .AND.                        &
     925                              BTEST(wall_flags_total_0(k,j-1,i),0)                                 &
     926                      )                       .OR.                                                 &
     927                      ( .NOT. BTEST(wall_flags_total_0(k,j-2,i),0)    .AND.                        &
     928                              BTEST(wall_flags_total_0(k,j+1,i),0)    .AND.                        &
     929                              BTEST(wall_flags_total_0(k,j+2,i),0)    .AND.                        &
     930                              BTEST(wall_flags_total_0(k,j-1,i),0)                                 &
     931                      )                       .OR.                                                 &
     932                      ( non_cyclic_s  .AND.  j == 1    )  .OR.                                     &
     933                      ( non_cyclic_n  .AND.  j == ny-1 ) )                                         &
     934             THEN
     935                advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
     936!
     937!--          WS5
     938             ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),0)           .AND.                         &
     939                      BTEST(wall_flags_total_0(k,j+2,i),0)           .AND.                         &
     940                      BTEST(wall_flags_total_0(k,j+3,i),0)           .AND.                         &
     941                      BTEST(wall_flags_total_0(k,j-1,i),0)           .AND.                         &
     942                      BTEST(wall_flags_total_0(k,j-2,i),0) )                                       &
     943             THEN
     944                advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 )
     945             ENDIF
     946!
     947!--          Near topography, set horizontal advection scheme to 1st order for passive scalars, even
     948!--          if only one direction may be blocked by topography. These locations will be identified
     949!--          by wall_flags_total_0 bit 31. Note, since several modules define advection flags but
     950!--          may apply different scalar boundary conditions, bit 31 is temporarily stored on
     951!--          advc_flags.
     952!--          Moreover, note that this extended degradtion for passive scalars is not required for
     953!--          the vertical direction as there the monotonic limiter can be applied.
     954             IF ( PRESENT( extensive_degrad ) )  THEN
     955                IF ( extensive_degrad )  THEN
     956!
     957!--                At all grid points that are within a three-grid point range to topography, set
     958!--                1st-order scheme.
     959                   IF( BTEST( advc_flag(k,j,i), 31 ) )  THEN
     960!
     961!--                   Clear flags that might indicate higher-order advection along x- and
     962!--                   y-direction.
     963                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
     964                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
     965                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
     966                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
     967!
     968!--                   Set flags that indicate 1st-order advection along x- and y-direction.
     969                      advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
     970                      advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
     971                   ENDIF
     972!
     973!--                Adjacent to this extended degradation zone, successively upgrade the order of the
     974!--                scheme if this grid point isn't flagged with bit 31 (indicating extended
     975!--                degradation zone).
     976                   IF ( .NOT. BTEST( advc_flag(k,j,i), 31 ) )  THEN
     977!
     978!--                   x-direction. First, clear all previous settings, then set flag for 3rd-order
     979!--                   scheme.
     980                      IF ( BTEST( advc_flag(k,j,i-1), 31 )  .AND.                                  &
     981                           BTEST( advc_flag(k,j,i+1), 31 ) )  THEN
     982                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
    1001983                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
    1002984                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
     985
     986                         advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
     987                      ENDIF
     988!
     989!--                   x-direction. First, clear all previous settings, then set flag for 5rd-order
     990!--                   scheme.
     991                      IF ( .NOT. BTEST( advc_flag(k,j,i-1), 31 )  .AND.                            &
     992                                 BTEST( advc_flag(k,j,i-2), 31 )  .AND.                            &
     993                           .NOT. BTEST( advc_flag(k,j,i+1), 31 )  .AND.                            &
     994                                 BTEST( advc_flag(k,j,i+2), 31 ) )  THEN
     995                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
     996                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
     997                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
     998
     999                         advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 )
     1000                      ENDIF
     1001!
     1002!--                   y-direction. First, clear all previous settings, then set flag for 3rd-order
     1003!--                   scheme.
     1004                      IF ( BTEST( advc_flag(k,j-1,i), 31 )  .AND.                                  &
     1005                           BTEST( advc_flag(k,j+1,i), 31 ) )  THEN
     1006                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
    10031007                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
    10041008                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
    1005 !
    1006 !--                      Set flags that indicate 1st-order advection along
    1007 !--                      x- and y-direction.
    1008                          advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
    1009                          advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
     1009
     1010                         advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
    10101011                      ENDIF
    10111012!
    1012 !--                   Adjacent to this extended degradation zone, successively
    1013 !--                   upgrade the order of the scheme if this grid point isn't
    1014 !--                   flagged with bit 31 (indicating extended degradation
    1015 !--                   zone).
    1016                       IF ( .NOT. BTEST( advc_flag(k,j,i), 31 ) )  THEN
    1017 !
    1018 !--                      x-direction. First, clear all previous settings, than
    1019 !--                      set flag for 3rd-order scheme.
    1020                          IF ( BTEST( advc_flag(k,j,i-1), 31 )  .AND.           &
    1021                               BTEST( advc_flag(k,j,i+1), 31 ) )  THEN
    1022                             advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
    1023                             advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
    1024                             advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
    1025 
    1026                             advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
    1027                          ENDIF
    1028 !
    1029 !--                      x-direction. First, clear all previous settings, than
    1030 !--                      set flag for 5rd-order scheme.
    1031                          IF ( .NOT. BTEST( advc_flag(k,j,i-1), 31 )  .AND.     &
    1032                                     BTEST( advc_flag(k,j,i-2), 31 )  .AND.     &
    1033                               .NOT. BTEST( advc_flag(k,j,i+1), 31 )  .AND.     &
    1034                                     BTEST( advc_flag(k,j,i+2), 31 ) )  THEN
    1035                             advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
    1036                             advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
    1037                             advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
    1038 
    1039                             advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 )
    1040                          ENDIF
    1041 !
    1042 !--                      y-direction. First, clear all previous settings, than
    1043 !--                      set flag for 3rd-order scheme.
    1044                          IF ( BTEST( advc_flag(k,j-1,i), 31 )  .AND.           &
    1045                               BTEST( advc_flag(k,j+1,i), 31 ) )  THEN
    1046                             advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
    1047                             advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
    1048                             advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
    1049 
    1050                             advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
    1051                          ENDIF
    1052 !
    1053 !--                      y-direction. First, clear all previous settings, than
    1054 !--                      set flag for 5rd-order scheme.
    1055                          IF ( .NOT. BTEST( advc_flag(k,j-1,i), 31 )  .AND.     &
    1056                                     BTEST( advc_flag(k,j-2,i), 31 )  .AND.     &
    1057                               .NOT. BTEST( advc_flag(k,j+1,i), 31 )  .AND.     &
    1058                                     BTEST( advc_flag(k,j+2,i), 31 ) )  THEN
    1059                             advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
    1060                             advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
    1061                             advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
    1062 
    1063                             advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 )
    1064                          ENDIF
     1013!--                   y-direction. First, clear all previous settings, then set flag for 5rd-order
     1014!--                   scheme.
     1015                      IF ( .NOT. BTEST( advc_flag(k,j-1,i), 31 )  .AND.                            &
     1016                                 BTEST( advc_flag(k,j-2,i), 31 )  .AND.                            &
     1017                           .NOT. BTEST( advc_flag(k,j+1,i), 31 )  .AND.                            &
     1018                                 BTEST( advc_flag(k,j+2,i), 31 ) )  THEN
     1019                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
     1020                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
     1021                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
     1022
     1023                         advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 )
    10651024                      ENDIF
    1066 
    10671025                   ENDIF
    10681026
    1069 !
    1070 !--                Near lateral boundary flags might be overwritten. Set
    1071 !--                them again.
    1072 !--                x-direction
    1073                    IF ( ( non_cyclic_l  .AND.  i == 0  )  .OR.                 &
    1074                         ( non_cyclic_r  .AND.  i == nx ) )  THEN
    1075                       advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
    1076                       advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
    1077                       advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
    1078 
    1079                       advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
    1080                    ENDIF
    1081 
    1082                    IF ( ( non_cyclic_l  .AND.  i == 1    )  .OR.               &
    1083                         ( non_cyclic_r  .AND.  i == nx-1 ) )  THEN
    1084                       advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
    1085                       advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
    1086                       advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
    1087 
    1088                       advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
    1089                    ENDIF
    1090 !
    1091 !--                y-direction
    1092                    IF ( ( non_cyclic_n  .AND.  j == 0  )  .OR.                 &
    1093                         ( non_cyclic_s  .AND.  j == ny ) )  THEN
    1094                       advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
    1095                       advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
    1096                       advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
    1097 
    1098                       advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
    1099                    ENDIF
    1100 
    1101                    IF ( ( non_cyclic_n  .AND.  j == 1    )  .OR.               &
    1102                         ( non_cyclic_s  .AND.  j == ny-1 ) )  THEN
    1103                       advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
    1104                       advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
    1105                       advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
    1106 
    1107                       advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
    1108                    ENDIF
    1109 
    11101027                ENDIF
    11111028
    1112 
    1113 !
    1114 !--             scalar - z-direction. Fluxes are calculated on w-grid
    1115 !--             level. Boundary values at/within walls aren't used.
    1116 !--             WS1 (6), WS3 (7), WS5 (8)
    1117                 IF ( k == nzb+1 )  THEN
    1118                    k_mm = nzb
    1119                 ELSE
    1120                    k_mm = k - 2
     1029!
     1030!--             Near lateral boundary flags might be overwritten. Set them again.
     1031!--             x-direction
     1032                IF ( ( non_cyclic_l  .AND.  i == 0  )  .OR.                                        &
     1033                     ( non_cyclic_r  .AND.  i == nx ) )  THEN
     1034                   advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
     1035                   advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
     1036                   advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
     1037
     1038                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
    11211039                ENDIF
    1122                 IF ( k > nzt-1 )  THEN
    1123                    k_pp = nzt+1
    1124                 ELSE
    1125                    k_pp = k + 2
     1040
     1041                IF ( ( non_cyclic_l  .AND.  i == 1    )  .OR.                                      &
     1042                     ( non_cyclic_r  .AND.  i == nx-1 ) )  THEN
     1043                   advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
     1044                   advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
     1045                   advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
     1046
     1047                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
    11261048                ENDIF
    1127                 IF ( k > nzt-2 )  THEN
    1128                    k_ppp = nzt+1
    1129                 ELSE
    1130                    k_ppp = k + 3
     1049!
     1050!--             y-direction
     1051                IF ( ( non_cyclic_n  .AND.  j == 0  )  .OR.                                        &
     1052                     ( non_cyclic_s  .AND.  j == ny ) )  THEN
     1053                   advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
     1054                   advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
     1055                   advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
     1056
     1057                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
    11311058                ENDIF
    11321059
    1133                 flag_set = .FALSE.
    1134                 IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),0)       .AND.  &
    1135                              BTEST(wall_flags_total_0(k,j,i),0)         .AND.  &
    1136                              BTEST(wall_flags_total_0(k+1,j,i),0) )     .OR.   &
    1137                      ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),0)      .AND.  &
    1138                              BTEST(wall_flags_total_0(k+1,j,i),0)       .AND.  &
    1139                              BTEST(wall_flags_total_0(k,j,i),0) )       .OR.   &
    1140                      ( k == nzt .AND. symmetry_flag == 0 ) )                   &
    1141                 THEN
    1142                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 6 )
    1143                    flag_set = .TRUE.
    1144                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),0)    .OR. &
    1145                            .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),0) ) .AND.&
    1146                                  BTEST(wall_flags_total_0(k-1,j,i),0)     .AND.&
    1147                                  BTEST(wall_flags_total_0(k,j,i),0)       .AND.&
    1148                                  BTEST(wall_flags_total_0(k+1,j,i),0)     .AND.&
    1149                                  BTEST(wall_flags_total_0(k_pp,j,i),0)    .AND.&
    1150                            .NOT. flag_set                                  .OR.&
    1151                          ( k == nzt - 1 .AND. symmetry_flag == 0 ) )           &
    1152                 THEN
    1153                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 7 )
    1154                    flag_set = .TRUE.
    1155                 ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),0)         .AND.   &
    1156                          BTEST(wall_flags_total_0(k-1,j,i),0)          .AND.   &
    1157                          BTEST(wall_flags_total_0(k,j,i),0)            .AND.   &
    1158                          BTEST(wall_flags_total_0(k+1,j,i),0)          .AND.   &
    1159                          BTEST(wall_flags_total_0(k_pp,j,i),0)         .AND.   &
    1160                          BTEST(wall_flags_total_0(k_ppp,j,i),0)        .AND.   &
    1161                         .NOT. flag_set )                                       &
    1162                 THEN
    1163                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 8 )
     1060                IF ( ( non_cyclic_n  .AND.  j == 1    )  .OR.                                      &
     1061                     ( non_cyclic_s  .AND.  j == ny-1 ) )  THEN
     1062                   advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
     1063                   advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
     1064                   advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
     1065
     1066                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
    11641067                ENDIF
    11651068
    1166              ENDDO
     1069             ENDIF
     1070
     1071
     1072!
     1073!--          scalar - z-direction. Fluxes are calculated on w-grid level. Boundary values at/within
     1074!--          walls aren't used.
     1075!--          WS1 (6), WS3 (7), WS5 (8)
     1076             IF ( k == nzb+1 )  THEN
     1077                k_mm = nzb
     1078             ELSE
     1079                k_mm = k - 2
     1080             ENDIF
     1081             IF ( k > nzt-1 )  THEN
     1082                k_pp = nzt+1
     1083             ELSE
     1084                k_pp = k + 2
     1085             ENDIF
     1086             IF ( k > nzt-2 )  THEN
     1087                k_ppp = nzt+1
     1088             ELSE
     1089                k_ppp = k + 3
     1090             ENDIF
     1091
     1092             flag_set = .FALSE.
     1093             IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),0)       .AND.                         &
     1094                          BTEST(wall_flags_total_0(k,j,i),0)         .AND.                         &
     1095                          BTEST(wall_flags_total_0(k+1,j,i),0) )     .OR.                          &
     1096                  ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),0)      .AND.                         &
     1097                          BTEST(wall_flags_total_0(k+1,j,i),0)       .AND.                         &
     1098                          BTEST(wall_flags_total_0(k,j,i),0) )       .OR.                          &
     1099                  ( k == nzt .AND. symmetry_flag == 0 ) )                                          &
     1100             THEN
     1101                advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 6 )
     1102                flag_set = .TRUE.
     1103             ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),0)    .OR.                        &
     1104                        .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),0) ) .AND.                       &
     1105                              BTEST(wall_flags_total_0(k-1,j,i),0)     .AND.                       &
     1106                              BTEST(wall_flags_total_0(k,j,i),0)       .AND.                       &
     1107                              BTEST(wall_flags_total_0(k+1,j,i),0)     .AND.                       &
     1108                              BTEST(wall_flags_total_0(k_pp,j,i),0)    .AND.                       &
     1109                        .NOT. flag_set                                  .OR.                       &
     1110                      ( k == nzt - 1 .AND. symmetry_flag == 0 ) )                                  &
     1111             THEN
     1112                advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 7 )
     1113                flag_set = .TRUE.
     1114             ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),0)         .AND.                          &
     1115                      BTEST(wall_flags_total_0(k-1,j,i),0)          .AND.                          &
     1116                      BTEST(wall_flags_total_0(k,j,i),0)            .AND.                          &
     1117                      BTEST(wall_flags_total_0(k+1,j,i),0)          .AND.                          &
     1118                      BTEST(wall_flags_total_0(k_pp,j,i),0)         .AND.                          &
     1119                      BTEST(wall_flags_total_0(k_ppp,j,i),0)        .AND.                          &
     1120                     .NOT. flag_set )                                                              &
     1121             THEN
     1122                advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 8 )
     1123             ENDIF
     1124
    11671125          ENDDO
    11681126       ENDDO
    1169 !
    1170 !--    Exchange 3D integer wall_flags.
    1171 !
    1172 !--    Exchange ghost points for advection flags
    1173        CALL exchange_horiz_int( advc_flag, nys, nyn, nxl, nxr, nzt, nbgp )
    1174 !
    1175 !--    Set boundary flags at inflow and outflow boundary in case of
    1176 !--   non-cyclic boundary conditions.
    1177        IF ( non_cyclic_l )  THEN
    1178           advc_flag(:,:,nxl-1) = advc_flag(:,:,nxl)
    1179        ENDIF
    1180 
    1181        IF ( non_cyclic_r )  THEN
    1182          advc_flag(:,:,nxr+1) = advc_flag(:,:,nxr)
    1183        ENDIF
    1184 
    1185        IF ( non_cyclic_n )  THEN
    1186           advc_flag(:,nyn+1,:) = advc_flag(:,nyn,:)
    1187        ENDIF
    1188 
    1189        IF ( non_cyclic_s )  THEN
    1190           advc_flag(:,nys-1,:) = advc_flag(:,nys,:)
    1191        ENDIF
    1192 
    1193 
    1194 
    1195     END SUBROUTINE ws_init_flags_scalar
    1196 
    1197 !------------------------------------------------------------------------------!
     1127    ENDDO
     1128!
     1129!-- Exchange 3D integer wall_flags.
     1130!
     1131!-- Exchange ghost points for advection flags
     1132    CALL exchange_horiz_int( advc_flag, nys, nyn, nxl, nxr, nzt, nbgp )
     1133!
     1134!-- Set boundary flags at inflow and outflow boundary in case of non-cyclic boundary conditions.
     1135    IF ( non_cyclic_l )  THEN
     1136       advc_flag(:,:,nxl-1) = advc_flag(:,:,nxl)
     1137    ENDIF
     1138
     1139    IF ( non_cyclic_r )  THEN
     1140      advc_flag(:,:,nxr+1) = advc_flag(:,:,nxr)
     1141    ENDIF
     1142
     1143    IF ( non_cyclic_n )  THEN
     1144       advc_flag(:,nyn+1,:) = advc_flag(:,nyn,:)
     1145    ENDIF
     1146
     1147    IF ( non_cyclic_s )  THEN
     1148       advc_flag(:,nys-1,:) = advc_flag(:,nys,:)
     1149    ENDIF
     1150
     1151
     1152
     1153 END SUBROUTINE ws_init_flags_scalar
     1154
     1155!--------------------------------------------------------------------------------------------------!
    11981156! Description:
    11991157! ------------
    12001158!> Initialize variables used for storing statistic quantities (fluxes, variances)
    1201 !------------------------------------------------------------------------------!
    1202     SUBROUTINE ws_statistics
    1203 
    1204 
    1205 !
    1206 !--    The arrays needed for statistical evaluation are set to to 0 at the
    1207 !--    beginning of prognostic_equations.
    1208        IF ( ws_scheme_mom )  THEN
    1209           !$ACC KERNELS PRESENT(sums_wsus_ws_l, sums_wsvs_ws_l) &
    1210           !$ACC PRESENT(sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l)
    1211           sums_wsus_ws_l = 0.0_wp
    1212           sums_wsvs_ws_l = 0.0_wp
    1213           sums_us2_ws_l  = 0.0_wp
    1214           sums_vs2_ws_l  = 0.0_wp
    1215           sums_ws2_ws_l  = 0.0_wp
    1216           !$ACC END KERNELS
    1217        ENDIF
    1218 
    1219        IF ( ws_scheme_sca )  THEN
    1220           !$ACC KERNELS PRESENT(sums_wspts_ws_l)
    1221           sums_wspts_ws_l = 0.0_wp
    1222           !$ACC END KERNELS
    1223           IF ( humidity       )  sums_wsqs_ws_l = 0.0_wp
    1224           IF ( passive_scalar )  sums_wsss_ws_l = 0.0_wp
    1225 
    1226        ENDIF
    1227 
    1228     END SUBROUTINE ws_statistics
    1229 
    1230 
    1231 !------------------------------------------------------------------------------!
     1159!--------------------------------------------------------------------------------------------------!
     1160 SUBROUTINE ws_statistics
     1161
     1162
     1163!
     1164!-- The arrays needed for statistical evaluation are set to to 0 at the beginning of
     1165!-- prognostic_equations.
     1166    IF ( ws_scheme_mom )  THEN
     1167       !$ACC KERNELS PRESENT(sums_wsus_ws_l, sums_wsvs_ws_l) &
     1168       !$ACC PRESENT(sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l)
     1169       sums_wsus_ws_l = 0.0_wp
     1170       sums_wsvs_ws_l = 0.0_wp
     1171       sums_us2_ws_l  = 0.0_wp
     1172       sums_vs2_ws_l  = 0.0_wp
     1173       sums_ws2_ws_l  = 0.0_wp
     1174       !$ACC END KERNELS
     1175    ENDIF
     1176
     1177    IF ( ws_scheme_sca )  THEN
     1178       !$ACC KERNELS PRESENT(sums_wspts_ws_l)
     1179       sums_wspts_ws_l = 0.0_wp
     1180       !$ACC END KERNELS
     1181       IF ( humidity       )  sums_wsqs_ws_l = 0.0_wp
     1182       IF ( passive_scalar )  sums_wsss_ws_l = 0.0_wp
     1183
     1184    ENDIF
     1185
     1186 END SUBROUTINE ws_statistics
     1187
     1188
     1189!--------------------------------------------------------------------------------------------------!
    12321190! Description:
    12331191! ------------
    12341192!> Scalar advection - Call for grid point i,j
    1235 !------------------------------------------------------------------------------!
    1236     SUBROUTINE advec_s_ws_ij( advc_flag, i, j, sk, sk_char, swap_flux_y_local, &
    1237                               swap_diss_y_local, swap_flux_x_local,            &
    1238                               swap_diss_x_local, i_omp, tn,                    &
    1239                               non_cyclic_l, non_cyclic_n,                      &
    1240                               non_cyclic_r, non_cyclic_s,                      &
    1241                               flux_limitation )
    1242 
    1243 
    1244        CHARACTER (LEN = *), INTENT(IN) ::  sk_char !< string identifier, used for assign fluxes to the
    1245                                                    !<correct dimension in the analysis array
    1246 
    1247        INTEGER(iwp) ::  i         !< grid index along x-direction
    1248        INTEGER(iwp) ::  i_omp     !< leftmost index on subdomain, or in case of OpenMP, on thread
    1249        INTEGER(iwp) ::  j         !< grid index along y-direction
    1250        INTEGER(iwp) ::  k         !< grid index along z-direction
    1251        INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
    1252        INTEGER(iwp) ::  k_mmm     !< k-3 index in disretization, can be modified to avoid segmentation faults
    1253        INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
    1254        INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
    1255        INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
    1256        INTEGER(iwp) ::  tn        !< number of OpenMP thread
    1257 
    1258        INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::   &
    1259                                                   advc_flag !< flag array to control order of scalar advection
    1260 
    1261        LOGICAL           ::  non_cyclic_l    !< flag that indicates non-cyclic boundary on the left
    1262        LOGICAL           ::  non_cyclic_n    !< flag that indicates non-cyclic boundary on the north
    1263        LOGICAL           ::  non_cyclic_r    !< flag that indicates non-cyclic boundary on the right
    1264        LOGICAL           ::  non_cyclic_s    !< flag that indicates non-cyclic boundary on the south
    1265        LOGICAL, OPTIONAL ::  flux_limitation !< flag indicating flux limitation of the vertical advection
    1266        LOGICAL           ::  limiter         !< control flag indicating the application of flux limitation
    1267 
    1268        REAL(wp) ::  diss_d        !< artificial dissipation term at grid box bottom
    1269        REAL(wp) ::  div           !< velocity diverence on scalar grid
    1270        REAL(wp) ::  div_in        !< vertical flux divergence of ingoing fluxes
    1271        REAL(wp) ::  div_out       !< vertical flux divergence of outgoing fluxes
    1272        REAL(wp) ::  f_corr_t      !< correction flux at grid-cell top, i.e. the difference between high and low-order flux
    1273        REAL(wp) ::  f_corr_d      !< correction flux at grid-cell bottom, i.e. the difference between high and low-order flux
    1274        REAL(wp) ::  f_corr_t_in   !< correction flux of ingoing flux part at grid-cell top
    1275        REAL(wp) ::  f_corr_d_in   !< correction flux of ingoing flux part at grid-cell bottom
    1276        REAL(wp) ::  f_corr_t_out  !< correction flux of outgoing flux part at grid-cell top
    1277        REAL(wp) ::  f_corr_d_out  !< correction flux of outgoing flux part at grid-cell bottom
    1278        REAL(wp) ::  fac_correction!< factor to limit the in- and outgoing fluxes
    1279        REAL(wp) ::  flux_d        !< 6th-order flux at grid box bottom
    1280        REAL(wp) ::  ibit0         !< flag indicating 1st-order scheme along x-direction
    1281        REAL(wp) ::  ibit1         !< flag indicating 3rd-order scheme along x-direction
    1282        REAL(wp) ::  ibit2         !< flag indicating 5th-order scheme along x-direction
    1283        REAL(wp) ::  ibit3         !< flag indicating 1st-order scheme along y-direction
    1284        REAL(wp) ::  ibit4         !< flag indicating 3rd-order scheme along y-direction
    1285        REAL(wp) ::  ibit5         !< flag indicating 5th-order scheme along y-direction
    1286        REAL(wp) ::  ibit6         !< flag indicating 1st-order scheme along z-direction
    1287        REAL(wp) ::  ibit7         !< flag indicating 3rd-order scheme along z-direction
    1288        REAL(wp) ::  ibit8         !< flag indicating 5th-order scheme along z-direction
    1289        REAL(wp) ::  max_val       !< maximum value of the quanitity along the numerical stencil (in vertical direction)
    1290        REAL(wp) ::  min_val       !< maximum value of the quanitity along the numerical stencil (in vertical direction)
    1291        REAL(wp) ::  mon           !< monotone solution of the advection equation using 1st-order fluxes
    1292        REAL(wp) ::  u_comp        !< advection velocity along x-direction
    1293        REAL(wp) ::  v_comp        !< advection velocity along y-direction
    1294 !
    1295 !--    sk is an array from parameter list. It should not be a pointer, because
    1296 !--    in that case the compiler can not assume a stride 1 and cannot perform
    1297 !--    a strided one vector load. Adding the CONTIGUOUS keyword makes things
    1298 !--    even worse, because the compiler cannot assume strided one in the
    1299 !--    caller side.
    1300        REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !<  advected scalar
    1301 
    1302        REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_n     !< discretized artificial dissipation at northward-side
    1303        REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_r     !< discretized artificial dissipation at rightward-side
    1304        REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_t     !< discretized artificial dissipation at top
    1305        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_n     !< discretized 6th-order flux at northward-side
    1306        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_r     !< discretized 6th-order flux at rightward-side
    1307        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t     !< discretized 6th-order flux at top
    1308        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t_1st !< discretized 1st-order flux at top
    1309 
    1310        REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_diss_y_local !< discretized artificial dissipation at southward-side
    1311        REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_flux_y_local !< discretized 6th-order flux at northward-side
    1312        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_diss_x_local !< discretized artificial dissipation at leftward-side
    1313        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_flux_x_local !< discretized 6th-order flux at leftward-side
    1314 !
    1315 !--    Used local modified copy of nzb_max (used to degrade order of
    1316 !--    discretization) at non-cyclic boundaries. Modify only at relevant points
    1317 !--    instead of the entire subdomain. This should lead to better
    1318 !--    load balance between boundary and non-boundary PEs.
    1319        IF( non_cyclic_l  .AND.  i <= nxl + 2  .OR.                             &
    1320            non_cyclic_r  .AND.  i >= nxr - 2  .OR.                             &
    1321            non_cyclic_s  .AND.  j <= nys + 2  .OR.                             &
    1322            non_cyclic_n  .AND.  j >= nyn - 2 )  THEN
    1323           nzb_max_l = nzt
    1324        ELSE
    1325           nzb_max_l = nzb_max
    1326        END IF
    1327 !
    1328 !--    Set control flag for flux limiter
    1329        limiter = .FALSE.
    1330        IF ( PRESENT( flux_limitation) )  limiter = flux_limitation
    1331 !
    1332 !--    Compute southside fluxes of the respective PE bounds.
    1333        IF ( j == nys )  THEN
    1334 !
    1335 !--       Up to the top of the highest topography.
    1336           DO  k = nzb+1, nzb_max_l
    1337 
    1338              ibit5 = REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp )
    1339              ibit4 = REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp )
    1340              ibit3 = REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp )
    1341 
    1342              v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
    1343              swap_flux_y_local(k,tn) = v_comp *         (                      &
    1344                                                ( 37.0_wp * ibit5 * adv_sca_5   &
    1345                                             +     7.0_wp * ibit4 * adv_sca_3   &
    1346                                             +              ibit3 * adv_sca_1   &
    1347                                                ) *                             &
    1348                                            ( sk(k,j,i)  + sk(k,j-1,i)     )    &
    1349                                          -     (  8.0_wp * ibit5 * adv_sca_5   &
    1350                                             +              ibit4 * adv_sca_3   &
    1351                                                 ) *                            &
    1352                                            ( sk(k,j+1,i) + sk(k,j-2,i)    )    &
    1353                                          +     (           ibit5 * adv_sca_5   &
    1354                                                ) *                             &
    1355                                            ( sk(k,j+2,i) + sk(k,j-3,i)    )    &
    1356                                                         )
    1357 
    1358              swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                      &
    1359                                                ( 10.0_wp * ibit5 * adv_sca_5   &
    1360                                             +     3.0_wp * ibit4 * adv_sca_3   &
    1361                                             +              ibit3 * adv_sca_1   &
    1362                                                ) *                             &
    1363                                             ( sk(k,j,i)   - sk(k,j-1,i)  )     &
    1364                                         -      (  5.0_wp * ibit5 * adv_sca_5   &
    1365                                             +              ibit4 * adv_sca_3   &
    1366                                             ) *                                &
    1367                                             ( sk(k,j+1,i) - sk(k,j-2,i)  )     &
    1368                                         +      (           ibit5 * adv_sca_5   &
    1369                                                ) *                             &
    1370                                             ( sk(k,j+2,i) - sk(k,j-3,i)  )     &
    1371                                                         )
    1372 
     1193!--------------------------------------------------------------------------------------------------!
     1194 SUBROUTINE advec_s_ws_ij( advc_flag, i, j, sk, sk_char, swap_flux_y_local, swap_diss_y_local,     &
     1195                           swap_flux_x_local, swap_diss_x_local, i_omp, tn, non_cyclic_l,          &
     1196                           non_cyclic_n, non_cyclic_r, non_cyclic_s, flux_limitation )
     1197
     1198
     1199    CHARACTER (LEN = *), INTENT(IN) ::  sk_char !< string identifier, used for assign fluxes to the
     1200                                                !<correct dimension in the analysis array
     1201
     1202    INTEGER(iwp) ::  i         !< grid index along x-direction
     1203    INTEGER(iwp) ::  i_omp     !< leftmost index on subdomain, or in case of OpenMP, on thread
     1204    INTEGER(iwp) ::  j         !< grid index along y-direction
     1205    INTEGER(iwp) ::  k         !< grid index along z-direction
     1206    INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
     1207    INTEGER(iwp) ::  k_mmm     !< k-3 index in disretization, can be modified to avoid segmentation faults
     1208    INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
     1209    INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
     1210    INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
     1211    INTEGER(iwp) ::  tn        !< number of OpenMP thread
     1212
     1213    INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::                          &
     1214                                        advc_flag !< flag array to control order of scalar advection
     1215
     1216    LOGICAL           ::  limiter         !< control flag indicating the application of flux limitation
     1217    LOGICAL           ::  non_cyclic_l    !< flag that indicates non-cyclic boundary on the left
     1218    LOGICAL           ::  non_cyclic_n    !< flag that indicates non-cyclic boundary on the north
     1219    LOGICAL           ::  non_cyclic_r    !< flag that indicates non-cyclic boundary on the right
     1220    LOGICAL           ::  non_cyclic_s    !< flag that indicates non-cyclic boundary on the south
     1221    LOGICAL, OPTIONAL ::  flux_limitation !< flag indicating flux limitation of the vertical advection
     1222
     1223    REAL(wp) ::  diss_d        !< artificial dissipation term at grid box bottom
     1224    REAL(wp) ::  div           !< velocity diverence on scalar grid
     1225    REAL(wp) ::  div_in        !< vertical flux divergence of ingoing fluxes
     1226    REAL(wp) ::  div_out       !< vertical flux divergence of outgoing fluxes
     1227    REAL(wp) ::  f_corr_d      !< correction flux at grid-cell bottom, i.e. the difference between high and low-order flux
     1228    REAL(wp) ::  f_corr_t      !< correction flux at grid-cell top, i.e. the difference between high and low-order flux
     1229    REAL(wp) ::  f_corr_d_in   !< correction flux of ingoing flux part at grid-cell bottom
     1230    REAL(wp) ::  f_corr_t_in   !< correction flux of ingoing flux part at grid-cell top
     1231    REAL(wp) ::  f_corr_d_out  !< correction flux of outgoing flux part at grid-cell bottom
     1232    REAL(wp) ::  f_corr_t_out  !< correction flux of outgoing flux part at grid-cell top
     1233    REAL(wp) ::  fac_correction!< factor to limit the in- and outgoing fluxes
     1234    REAL(wp) ::  flux_d        !< 6th-order flux at grid box bottom
     1235    REAL(wp) ::  ibit0         !< flag indicating 1st-order scheme along x-direction
     1236    REAL(wp) ::  ibit1         !< flag indicating 3rd-order scheme along x-direction
     1237    REAL(wp) ::  ibit2         !< flag indicating 5th-order scheme along x-direction
     1238    REAL(wp) ::  ibit3         !< flag indicating 1st-order scheme along y-direction
     1239    REAL(wp) ::  ibit4         !< flag indicating 3rd-order scheme along y-direction
     1240    REAL(wp) ::  ibit5         !< flag indicating 5th-order scheme along y-direction
     1241    REAL(wp) ::  ibit6         !< flag indicating 1st-order scheme along z-direction
     1242    REAL(wp) ::  ibit7         !< flag indicating 3rd-order scheme along z-direction
     1243    REAL(wp) ::  ibit8         !< flag indicating 5th-order scheme along z-direction
     1244    REAL(wp) ::  max_val       !< maximum value of the quanitity along the numerical stencil (in vertical direction)
     1245    REAL(wp) ::  min_val       !< maximum value of the quanitity along the numerical stencil (in vertical direction)
     1246    REAL(wp) ::  mon           !< monotone solution of the advection equation using 1st-order fluxes
     1247    REAL(wp) ::  u_comp        !< advection velocity along x-direction
     1248    REAL(wp) ::  v_comp        !< advection velocity along y-direction
     1249
     1250    REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_n     !< discretized artificial dissipation at northward-side
     1251    REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_r     !< discretized artificial dissipation at rightward-side
     1252    REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_t     !< discretized artificial dissipation at top
     1253    REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_n     !< discretized 6th-order flux at northward-side
     1254    REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_r     !< discretized 6th-order flux at rightward-side
     1255    REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t     !< discretized 6th-order flux at top
     1256    REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t_1st !< discretized 1st-order flux at top
     1257
     1258    REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_diss_y_local !< discretized artificial dissipation at southward-side
     1259    REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_flux_y_local !< discretized 6th-order flux at northward-side
     1260    REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_diss_x_local !< discretized artificial dissipation at leftward-side
     1261    REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_flux_x_local !< discretized 6th-order flux at leftward-side
     1262
     1263!
     1264!-- sk is an array from parameter list. It should not be a pointer, because in that case the
     1265!-- compiler can not assume a stride 1 and cannot perform a strided one vector load. Adding the
     1266!-- CONTIGUOUS keyword makes things even worse, because the compiler cannot assume strided one in
     1267!-- the caller side.
     1268    REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !< advected scalar
     1269
     1270!
     1271!-- Used local modified copy of nzb_max (used to degrade order of discretization) at non-cyclic
     1272!-- boundaries. Modify only at relevant points instead of the entire subdomain. This should lead to
     1273!-- betterload balance between boundary and non-boundary PEs.
     1274    IF( non_cyclic_l  .AND.  i <= nxl + 2  .OR.                                                    &
     1275        non_cyclic_r  .AND.  i >= nxr - 2  .OR.                                                    &
     1276        non_cyclic_s  .AND.  j <= nys + 2  .OR.                                                    &
     1277        non_cyclic_n  .AND.  j >= nyn - 2 )  THEN
     1278       nzb_max_l = nzt
     1279    ELSE
     1280       nzb_max_l = nzb_max
     1281    END IF
     1282!
     1283!-- Set control flag for flux limiter
     1284    limiter = .FALSE.
     1285    IF ( PRESENT( flux_limitation) )  limiter = flux_limitation
     1286!
     1287!-- Compute southside fluxes of the respective PE bounds.
     1288    IF ( j == nys )  THEN
     1289!
     1290!--    Up to the top of the highest topography.
     1291       DO  k = nzb+1, nzb_max_l
     1292
     1293          ibit5 = REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp )
     1294          ibit4 = REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp )
     1295          ibit3 = REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp )
     1296
     1297          v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
     1298          swap_flux_y_local(k,tn) = v_comp *         (                                             &
     1299                                            ( 37.0_wp * ibit5 * adv_sca_5                          &
     1300                                         +     7.0_wp * ibit4 * adv_sca_3                          &
     1301                                         +              ibit3 * adv_sca_1                          &
     1302                                            ) * ( sk(k,j,i)   + sk(k,j-1,i) )                      &
     1303                                         -  (  8.0_wp * ibit5 * adv_sca_5                          &
     1304                                         +              ibit4 * adv_sca_3                          &
     1305                                            ) * ( sk(k,j+1,i) + sk(k,j-2,i) )                      &
     1306                                         +  (           ibit5 * adv_sca_5 )                        &
     1307                                              * ( sk(k,j+2,i) + sk(k,j-3,i) )                      &
     1308                                                     )
     1309
     1310          swap_diss_y_local(k,tn) = - ABS( v_comp ) * (                                            &
     1311                                            ( 10.0_wp * ibit5 * adv_sca_5                          &
     1312                                         +     3.0_wp * ibit4 * adv_sca_3                          &
     1313                                         +              ibit3 * adv_sca_1                          &
     1314                                            ) * ( sk(k,j,i)   - sk(k,j-1,i) )                      &
     1315                                         -  ( 5.0_wp  * ibit5 * adv_sca_5                          &
     1316                                         +              ibit4 * adv_sca_3                          &
     1317                                            ) * ( sk(k,j+1,i) - sk(k,j-2,i) )                      &
     1318                                         +  (           ibit5 * adv_sca_5   )                      &
     1319                                              * ( sk(k,j+2,i) - sk(k,j-3,i)  )                     &
     1320                                                      )
     1321
     1322       ENDDO
     1323!
     1324!--    Above to the top of the highest topography. No degradation necessary.
     1325       DO  k = nzb_max_l+1, nzt
     1326
     1327          v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
     1328          swap_flux_y_local(k,tn) = v_comp * ( 37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )             &
     1329                                              - 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) )             &
     1330                                              + ( sk(k,j+2,i) + sk(k,j-3,i) )                      &
     1331                                             ) * adv_sca_5
     1332          swap_diss_y_local(k,tn) = - ABS( v_comp ) * (                                            &
     1333                                     10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )                       &
     1334                                    - 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )                       &
     1335                                    +            sk(k,j+2,i) - sk(k,j-3,i)                         &
     1336                                                      ) * adv_sca_5
     1337
     1338       ENDDO
     1339
     1340    ENDIF
     1341!
     1342!-- Compute leftside fluxes of the respective PE bounds.
     1343    IF ( i == i_omp )  THEN
     1344
     1345       DO  k = nzb+1, nzb_max_l
     1346
     1347          ibit2 = REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp )
     1348          ibit1 = REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp )
     1349          ibit0 = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )
     1350
     1351          u_comp                    = u(k,j,i) - u_gtrans + u_stokes_zu(k)
     1352          swap_flux_x_local(k,j,tn) = u_comp * (                                                   &
     1353                                            ( 37.0_wp * ibit2 * adv_sca_5                          &
     1354                                         +     7.0_wp * ibit1 * adv_sca_3                          &
     1355                                         +              ibit0 * adv_sca_1                          &
     1356                                            ) *   ( sk(k,j,i) + sk(k,j,i-1) )                      &
     1357                                         -  (  8.0_wp * ibit2 * adv_sca_5                          &
     1358                                         +              ibit1 * adv_sca_3                          &
     1359                                            ) * ( sk(k,j,i+1) + sk(k,j,i-2) )                      &
     1360                                         +  (           ibit2 * adv_sca_5                          &
     1361                                            ) * ( sk(k,j,i+2) + sk(k,j,i-3) )                      &
     1362                                               )
     1363
     1364          swap_diss_x_local(k,j,tn) = - ABS( u_comp ) * (                                          &
     1365                                            ( 10.0_wp * ibit2 * adv_sca_5                          &
     1366                                         +     3.0_wp * ibit1 * adv_sca_3                          &
     1367                                         +              ibit0 * adv_sca_1                          &
     1368                                            ) *   ( sk(k,j,i) - sk(k,j,i-1) )                      &
     1369                                         -  (  5.0_wp * ibit2 * adv_sca_5                          &
     1370                                         +              ibit1 * adv_sca_3                          &
     1371                                            ) * ( sk(k,j,i+1) - sk(k,j,i-2) )                      &
     1372                                         +  (           ibit2 * adv_sca_5                          &
     1373                                            ) * ( sk(k,j,i+2) - sk(k,j,i-3)    )                   &
     1374                                                       )
     1375
     1376       ENDDO
     1377
     1378       DO  k = nzb_max_l+1, nzt
     1379
     1380          u_comp                    = u(k,j,i) - u_gtrans + u_stokes_zu(k)
     1381          swap_flux_x_local(k,j,tn) = u_comp  * (                                                  &
     1382                                        37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) )                      &
     1383                                      -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) )                    &
     1384                                      +           ( sk(k,j,i+2) + sk(k,j,i-3) )                    &
     1385                                                ) * adv_sca_5
     1386
     1387          swap_diss_x_local(k,j,tn) = - ABS( u_comp ) * (                                          &
     1388                                        10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) )                    &
     1389                                       - 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) )                    &
     1390                                       +          ( sk(k,j,i+2) - sk(k,j,i-3) )                    &
     1391                                                        ) * adv_sca_5
     1392
     1393       ENDDO
     1394
     1395    ENDIF
     1396!
     1397!-- Now compute the fluxes for the horizontal termns up to the highest
     1398!-- topography.
     1399    DO  k = nzb+1, nzb_max_l
     1400
     1401       ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
     1402       ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp )
     1403       ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp )
     1404
     1405       u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
     1406       flux_r(k) = u_comp * (                                                                      &
     1407                   (  37.0_wp * ibit2 * adv_sca_5                                                  &
     1408                   +   7.0_wp * ibit1 * adv_sca_3                                                  &
     1409                   +            ibit0 * adv_sca_1 ) * ( sk(k,j,i+1) + sk(k,j,i)   )                &
     1410                   - ( 8.0_wp * ibit2 * adv_sca_5                                                  &
     1411                   +            ibit1 * adv_sca_3 ) * ( sk(k,j,i+2) + sk(k,j,i-1) )                &
     1412                   + (          ibit2 * adv_sca_5 ) * ( sk(k,j,i+3) + sk(k,j,i-2) )                &
     1413                            )
     1414
     1415       diss_r(k) = - ABS( u_comp ) * (                                                             &
     1416                   (  10.0_wp * ibit2 * adv_sca_5                                                  &
     1417                   +   3.0_wp * ibit1 * adv_sca_3                                                  &
     1418                   +            ibit0 * adv_sca_1 ) * ( sk(k,j,i+1) - sk(k,j,i) )                  &
     1419                   - ( 5.0_wp * ibit2 * adv_sca_5                                                  &
     1420                   +            ibit1 * adv_sca_3 ) * ( sk(k,j,i+2) - sk(k,j,i-1) )                &
     1421                   + (          ibit2 * adv_sca_5 ) * ( sk(k,j,i+3) - sk(k,j,i-2) )                &
     1422                                     )
     1423
     1424       ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp )
     1425       ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp )
     1426       ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp )
     1427
     1428       v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
     1429       flux_n(k) = v_comp * (                                                                      &
     1430                   (  37.0_wp * ibit5 * adv_sca_5                                                  &
     1431                   +   7.0_wp * ibit4 * adv_sca_3                                                  &
     1432                   +            ibit3 * adv_sca_1 ) * ( sk(k,j+1,i) + sk(k,j,i)   )                &
     1433                   - ( 8.0_wp * ibit5 * adv_sca_5                                                  &
     1434                   +            ibit4 * adv_sca_3 ) * ( sk(k,j+2,i) + sk(k,j-1,i) )                &
     1435                   + (          ibit5 * adv_sca_5 ) * ( sk(k,j+3,i) + sk(k,j-2,i) )                &
     1436                            )
     1437
     1438       diss_n(k) = - ABS( v_comp ) * (                                                             &
     1439                   (  10.0_wp * ibit5 * adv_sca_5                                                  &
     1440                   +   3.0_wp * ibit4 * adv_sca_3                                                  &
     1441                   +            ibit3 * adv_sca_1 ) * ( sk(k,j+1,i) - sk(k,j,i)   )                &
     1442                   - ( 5.0_wp * ibit5 * adv_sca_5                                                  &
     1443                   +            ibit4 * adv_sca_3 ) * ( sk(k,j+2,i) - sk(k,j-1,i) )                &
     1444                   + (          ibit5 * adv_sca_5 ) * ( sk(k,j+3,i) - sk(k,j-2,i) )                &
     1445                                     )
     1446    ENDDO
     1447!
     1448!-- Now compute the fluxes for the horizontal terms above the topography
     1449!-- where no degradation along the horizontal parts is necessary (except
     1450!-- for the non-cyclic lateral boundaries treated by nzb_max_l).
     1451    DO  k = nzb_max_l+1, nzt
     1452
     1453       u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
     1454       flux_r(k) = u_comp * (                                                                      &
     1455                     37.0_wp * ( sk(k,j,i+1) + sk(k,j,i)   )                                       &
     1456                   -  8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) )                                       &
     1457                   +           ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
     1458       diss_r(k) = - ABS( u_comp ) * (                                                             &
     1459                     10.0_wp * ( sk(k,j,i+1) - sk(k,j,i)   )                                       &
     1460                   -  5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) )                                       &
     1461                   +           ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
     1462
     1463       v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
     1464       flux_n(k) = v_comp * (                                                                      &
     1465                     37.0_wp * ( sk(k,j+1,i) + sk(k,j,i)   )                                       &
     1466                   -  8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) )                                       &
     1467                   +           ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
     1468       diss_n(k) = - ABS( v_comp ) * (                                                             &
     1469                     10.0_wp * ( sk(k,j+1,i) - sk(k,j,i)   )                                       &
     1470                   -  5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) )                                       &
     1471                   +           ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
     1472
     1473    ENDDO
     1474!
     1475!-- Now, compute vertical fluxes. Split loop into a part treating the lowest grid points with
     1476!-- indirect indexing, a main loop without indirect indexing, and a loop for the uppermost grid
     1477!-- points with indirect indexing. This allows better vectorization for the main loop.
     1478!-- First, compute the flux at model surface, which need has to be calculated explicetely for the
     1479!-- tendency at the first w-level. For topography wall this is done implicitely by advc_flag.
     1480    flux_t(nzb) = 0.0_wp
     1481    diss_t(nzb) = 0.0_wp
     1482
     1483    DO  k = nzb+1, nzb+1
     1484       ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     1485       ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     1486       ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
     1487!
     1488!--    k index has to be modified near bottom and top, else array subscripts will be exceeded.
     1489       k_ppp = k + 3 * ibit8
     1490       k_pp  = k + 2 * ( 1 - ibit6 )
     1491       k_mm  = k - 2 * ibit8
     1492
     1493       flux_t(k) = w(k,j,i)  * rho_air_zw(k) * (                                                   &
     1494                   (  37.0_wp * ibit8 * adv_sca_5                                                  &
     1495                   +   7.0_wp * ibit7 * adv_sca_3                                                  &
     1496                   +            ibit6 * adv_sca_1 ) * ( sk(k+1,j,i)  + sk(k,j,i)    )              &
     1497                   - ( 8.0_wp * ibit8 * adv_sca_5                                                  &
     1498                   +            ibit7 * adv_sca_3 ) * ( sk(k_pp,j,i) + sk(k-1,j,i)  )              &
     1499                   + (          ibit8 * adv_sca_5 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )              &
     1500                                               )
     1501
     1502       diss_t(k) = - ABS( w(k,j,i) ) * rho_air_zw(k) * (                                           &
     1503                   (  10.0_wp * ibit8 * adv_sca_5                                                  &
     1504                   +   3.0_wp * ibit7 * adv_sca_3                                                  &
     1505                   +            ibit6 * adv_sca_1 ) * ( sk(k+1,j,i)   - sk(k,j,i)    )             &
     1506                   - ( 5.0_wp * ibit8 * adv_sca_5                                                  &
     1507                   +            ibit7 * adv_sca_3 ) * ( sk(k_pp,j,i)  - sk(k-1,j,i)  )             &
     1508                   + (          ibit8 * adv_sca_5 ) * ( sk(k_ppp,j,i) - sk(k_mm,j,i) )             &
     1509                                                       )
     1510    ENDDO
     1511
     1512    DO  k = nzb+2, nzt-2
     1513       ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     1514       ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     1515       ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
     1516
     1517       flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                                                    &
     1518                   (  37.0_wp * ibit8 * adv_sca_5                                                  &
     1519                   +   7.0_wp * ibit7 * adv_sca_3                                                  &
     1520                   +            ibit6 * adv_sca_1 ) * ( sk(k+1,j,i)  + sk(k,j,i)   )               &
     1521                   - ( 8.0_wp * ibit8 * adv_sca_5                                                  &
     1522                   +            ibit7 * adv_sca_3 ) * ( sk(k+2,j,i)  + sk(k-1,j,i) )               &
     1523                   + (          ibit8 * adv_sca_5 ) * ( sk(k+3,j,i)  + sk(k-2,j,i) )               &
     1524                                              )
     1525
     1526       diss_t(k) = - ABS( w(k,j,i) ) * rho_air_zw(k) * (                                           &
     1527                   (  10.0_wp * ibit8 * adv_sca_5                                                  &
     1528                   +   3.0_wp * ibit7 * adv_sca_3                                                  &
     1529                   +            ibit6 * adv_sca_1 ) * ( sk(k+1,j,i) - sk(k,j,i)   )                &
     1530                   - ( 5.0_wp * ibit8 * adv_sca_5                                                  &
     1531                   +            ibit7 * adv_sca_3 ) * ( sk(k+2,j,i) - sk(k-1,j,i) )                &
     1532                   + (          ibit8 * adv_sca_5 ) * ( sk(k+3,j,i) - sk(k-2,j,i) )                &
     1533                                                       )
     1534    ENDDO
     1535
     1536    DO  k = nzt-1, nzt-symmetry_flag
     1537       ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     1538       ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     1539       ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
     1540!
     1541!--    k index has to be modified near bottom and top, else array subscripts will be exceeded.
     1542       k_ppp = k + 3 * ibit8
     1543       k_pp  = k + 2 * ( 1 - ibit6  )
     1544       k_mm  = k - 2 * ibit8
     1545
     1546
     1547       flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                                                    &
     1548                   (  37.0_wp * ibit8 * adv_sca_5                                                  &
     1549                   +   7.0_wp * ibit7 * adv_sca_3                                                  &
     1550                   +            ibit6 * adv_sca_1 ) * ( sk(k+1,j,i)  + sk(k,j,i)    )              &
     1551                   - ( 8.0_wp * ibit8 * adv_sca_5                                                  &
     1552                   +            ibit7 * adv_sca_3 ) * ( sk(k_pp,j,i) + sk(k-1,j,i)  )              &
     1553                   + (          ibit8 * adv_sca_5 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )              &
     1554                                              )
     1555
     1556       diss_t(k) = - ABS( w(k,j,i) ) * rho_air_zw(k) * (                                           &
     1557                   (  10.0_wp * ibit8 * adv_sca_5                                                  &
     1558                   +   3.0_wp * ibit7 * adv_sca_3                                                  &
     1559                   +            ibit6 * adv_sca_1 ) * ( sk(k+1,j,i)   - sk(k,j,i)    )             &
     1560                   - ( 5.0_wp * ibit8 * adv_sca_5                                                  &
     1561                   +            ibit7 * adv_sca_3 ) * ( sk(k_pp,j,i)  - sk(k-1,j,i)  )             &
     1562                   + (          ibit8 * adv_sca_5 ) * ( sk(k_ppp,j,i) - sk(k_mm,j,i) )             &
     1563                                                       )
     1564    ENDDO
     1565
     1566!
     1567!-- Set resolved/turbulent flux at model top to zero (w-level). In case that a symmetric behavior
     1568!-- between bottom and top shall be guaranteed (closed channel flow), the flux at nzt is also set to
     1569!-- zero.
     1570    IF ( symmetry_flag == 1 ) THEN
     1571       flux_t(nzt) = 0.0_wp
     1572       diss_t(nzt) = 0.0_wp
     1573    ENDIF
     1574    flux_t(nzt+1) = 0.0_wp
     1575    diss_t(nzt+1) = 0.0_wp
     1576
     1577
     1578    IF ( limiter )  THEN
     1579!
     1580!--    Compute monotone first-order fluxes which are required for mononteflux limitation.
     1581       flux_t_1st(nzb) = 0.0_wp
     1582       DO  k = nzb+1, nzb_max_l
     1583          flux_t_1st(k) = (      w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )                          &
     1584                          - ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) )                        &
     1585                          * rho_air_zw(k) * adv_sca_1
     1586!
     1587!--       In flux limitation the total flux will be corrected. For the sake of cleariness the
     1588!--       higher-order advective and disspative fluxes will be merged onto flux_t.
     1589          flux_t(k) = flux_t(k) + diss_t(k)
     1590          diss_t(k) = 0.0_wp
     1591       ENDDO
     1592!
     1593!--    Flux limitation of vertical fluxes according to Skamarock (2006).
     1594!--    Please note, as flux limitation implies linear dependencies of fluxes, flux limitation is
     1595!--    only made for the vertical advection term. Limitation of the horizontal terms cannot be
     1596!--    parallelized.
     1597!--    Due to the linear dependency, the following loop will not be vectorized.
     1598!--    Further, note that the flux limiter is only applied within the urban layer, i.e up to the
     1599!--    topography top.
     1600       DO  k = nzb+1, nzb_max_l
     1601!
     1602!--       Compute one-dimensional divergence along the vertical direction, which is used to correct
     1603!--       the advection discretization. This is necessary as in one-dimensional space the advection
     1604!--       velocity should actually be constant.
     1605          div = ( w(k,j,i)   * rho_air_zw(k)                                                       &
     1606                - w(k-1,j,i) * rho_air_zw(k-1)                                                     &
     1607                ) * drho_air(k) * ddzw(k)
     1608!
     1609!--       Compute monotone solution of the advection equation from 1st-order fluxes. Please note,
     1610!--       the advection equation is corrected by the divergence term (in 1D the advective flow
     1611!--       should be divergence free). Moreover, please note, as time-increment the full timestep is
     1612!--       used, even though a Runge-Kutta scheme will be used. However, the length of the actual
     1613!--       time increment is not important at all since it cancels out later when the fluxes are
     1614!--       limited.
     1615          mon = sk(k,j,i) + ( - ( flux_t_1st(k) - flux_t_1st(k-1) )                                &
     1616                          * drho_air(k) * ddzw(k)                                                  &
     1617                          + div * sk(k,j,i)                                                        &
     1618                            ) * dt_3d
     1619!
     1620!--       Determine minimum and maximum values along the numerical stencil.
     1621          k_mmm = MAX( k - 3, nzb + 1 )
     1622          k_ppp = MIN( k + 3, nzt + 1 )
     1623
     1624          min_val = MINVAL( sk(k_mmm:k_ppp,j,i) )
     1625          max_val = MAXVAL( sk(k_mmm:k_ppp,j,i) )
     1626!
     1627!--       Compute difference between high- and low-order fluxes, which may act as correction fluxes
     1628          f_corr_t = flux_t(k)   - flux_t_1st(k)
     1629          f_corr_d = flux_t(k-1) - flux_t_1st(k-1)
     1630!
     1631!--       Determine outgoing fluxes, i.e. the part of the fluxes which can decrease the value within
     1632!--       the grid box
     1633          f_corr_t_out = MAX( 0.0_wp, f_corr_t )
     1634          f_corr_d_out = MIN( 0.0_wp, f_corr_d )
     1635!
     1636!--       Determine ingoing fluxes, i.e. the part of the fluxes which can  increase the value within
     1637!--       the grid box
     1638          f_corr_t_in = MIN( 0.0_wp, f_corr_t)
     1639          f_corr_d_in = MAX( 0.0_wp, f_corr_d)
     1640!
     1641!--       Compute divergence of outgoing correction fluxes
     1642          div_out = - ( f_corr_t_out - f_corr_d_out ) * drho_air(k) * ddzw(k) * dt_3d
     1643!
     1644!--       Compute divergence of ingoing correction fluxes
     1645          div_in = - ( f_corr_t_in - f_corr_d_in ) * drho_air(k) * ddzw(k) * dt_3d
     1646!
     1647!--       Check if outgoing fluxes can lead to undershoots, i.e. values smaller than the minimum
     1648!--       value within the numerical stencil. If so, limit them.
     1649          IF ( mon - min_val < - div_out  .AND.  ABS( div_out ) > 0.0_wp )  THEN
     1650             fac_correction = ( mon - min_val ) / ( - div_out )
     1651             f_corr_t_out = f_corr_t_out * fac_correction
     1652             f_corr_d_out = f_corr_d_out * fac_correction
     1653          ENDIF
     1654!
     1655!--       Check if ingoing fluxes can lead to overshoots, i.e. values larger than the maximum value
     1656!--       within the numerical stencil. If so, limit them.
     1657          IF ( mon - max_val > - div_in  .AND.  ABS( div_in ) > 0.0_wp )  THEN
     1658             fac_correction = ( mon - max_val ) / ( - div_in )
     1659             f_corr_t_in = f_corr_t_in * fac_correction
     1660             f_corr_d_in = f_corr_d_in * fac_correction
     1661          ENDIF
     1662!
     1663!--       Finally add the limited fluxes to the original ones. If no flux limitation was done, the
     1664!--       fluxes equal the original ones.
     1665          flux_t(k)   = flux_t_1st(k)   + f_corr_t_out + f_corr_t_in
     1666          flux_t(k-1) = flux_t_1st(k-1) + f_corr_d_out + f_corr_d_in
     1667       ENDDO
     1668    ENDIF
     1669!
     1670!-- Now compute the tendency term including divergence correction.
     1671    DO  k = nzb+1, nzb_max_l
     1672
     1673       flux_d    = flux_t(k-1)
     1674       diss_d    = diss_t(k-1)
     1675
     1676       ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
     1677       ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp )
     1678       ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp )
     1679
     1680       ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp )
     1681       ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp )
     1682       ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp )
     1683
     1684       ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     1685       ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     1686       ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
     1687!
     1688!--    Calculate the divergence of the velocity field. A respective correction is needed to overcome
     1689!--    numerical instabilities introduced by an insufficient reduction of divergences near
     1690!--    topography.
     1691       div         =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )                                    &
     1692                       - u(k,j,i)   * (                                                            &
     1693                       REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )                            &
     1694                     + REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp )                            &
     1695                     + REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp )                            &
     1696                                      )                                                            &
     1697                       ) * ddx                                                                     &
     1698                     + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )                                    &
     1699                       - v(k,j,i)   * (                                                            &
     1700                       REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp )                            &
     1701                     + REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp )                            &
     1702                     + REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp )                            &
     1703                                      )                                                            &
     1704                       ) * ddy                                                                     &
     1705                     + ( w(k,j,i)   * rho_air_zw(k) * ( ibit6 + ibit7 + ibit8 )                    &
     1706                       - w(k-1,j,i) * rho_air_zw(k-1) *                                            &
     1707                                      (                                                            &
     1708                       REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp )                            &
     1709                     + REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp )                            &
     1710                     + REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp )                            &
     1711                                      )                                                            &
     1712                       ) * drho_air(k) * ddzw(k)
     1713
     1714       tend(k,j,i) = tend(k,j,i) - (                                                               &
     1715                     ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn)                           &
     1716                     - swap_diss_x_local(k,j,tn) ) * ddx                                           &
     1717                     + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)                           &
     1718                     - swap_diss_y_local(k,tn)   ) * ddy                                           &
     1719                     + ( ( flux_t(k) + diss_t(k) ) - ( flux_d + diss_d )                           &
     1720                       )             * drho_air(k) * ddzw(k)                                       &
     1721                                   ) + sk(k,j,i) * div
     1722
     1723
     1724       swap_flux_y_local(k,tn)   = flux_n(k)
     1725       swap_diss_y_local(k,tn)   = diss_n(k)
     1726       swap_flux_x_local(k,j,tn) = flux_r(k)
     1727       swap_diss_x_local(k,j,tn) = diss_r(k)
     1728
     1729    ENDDO
     1730
     1731    DO  k = nzb_max_l+1, nzt
     1732
     1733       flux_d    = flux_t(k-1)
     1734       diss_d    = diss_t(k-1)
     1735!
     1736!--    Calculate the divergence of the velocity field. A respective correction is needed to overcome
     1737!--    numerical instabilities introduced by an insufficient reduction of divergences near
     1738!--    topography.
     1739       div         =   ( u(k,j,i+1) - u(k,j,i) ) * ddx                                             &
     1740                     + ( v(k,j+1,i) - v(k,j,i) ) * ddy                                             &
     1741                     + ( w(k,j,i)   * rho_air_zw(k)                                                &
     1742                       - w(k-1,j,i) * rho_air_zw(k-1)                                              &
     1743                                               ) * drho_air(k) * ddzw(k)
     1744
     1745       tend(k,j,i) = tend(k,j,i) - (                                                               &
     1746                     ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn)                           &
     1747                     - swap_diss_x_local(k,j,tn) ) * ddx                                           &
     1748                     + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)                           &
     1749                     - swap_diss_y_local(k,tn)   ) * ddy                                           &
     1750                     + ( ( flux_t(k) + diss_t(k) ) - ( flux_d + diss_d )                           &
     1751                                                 ) * drho_air(k) * ddzw(k)                         &
     1752                                   ) + sk(k,j,i) * div
     1753
     1754
     1755       swap_flux_y_local(k,tn)   = flux_n(k)
     1756       swap_diss_y_local(k,tn)   = diss_n(k)
     1757       swap_flux_x_local(k,j,tn) = flux_r(k)
     1758       swap_diss_x_local(k,j,tn) = diss_r(k)
     1759
     1760    ENDDO
     1761
     1762!
     1763!-- Evaluation of statistics.
     1764    SELECT CASE ( sk_char )
     1765
     1766       CASE ( 'pt' )
     1767
     1768          DO  k = nzb, nzt
     1769             sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +                                       &
     1770                                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )     &
     1771                                                 * ( w(k,j,i) - hom(k,1,3,0)                 )     &
     1772                                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )     &
     1773                                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )     &
     1774                                     ) * weight_substep(intermediate_timestep_count)
    13731775          ENDDO
    1374 !
    1375 !--       Above to the top of the highest topography. No degradation necessary.
    1376           DO  k = nzb_max_l+1, nzt
    1377 
    1378              v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
    1379              swap_flux_y_local(k,tn) = v_comp * (                              &
    1380                                     37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )    &
    1381                                   -  8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) )    &
    1382                                   +           ( sk(k,j+2,i) + sk(k,j-3,i) )    &
    1383                                                 ) * adv_sca_5
    1384              swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                      &
    1385                                     10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )    &
    1386                                   -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )    &
    1387                                   +             sk(k,j+2,i) - sk(k,j-3,i)      &
    1388                                                         ) * adv_sca_5
    1389 
     1776
     1777       CASE ( 'sa' )
     1778
     1779          DO  k = nzb, nzt
     1780             sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +                                       &
     1781                                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )     &
     1782                                                 * ( w(k,j,i)    - hom(k,1,3,0)              )     &
     1783                                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )     &
     1784                                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )     &
     1785                                     ) * weight_substep(intermediate_timestep_count)
    13901786          ENDDO
    13911787
    1392        ENDIF
    1393 !
    1394 !--    Compute leftside fluxes of the respective PE bounds.
    1395        IF ( i == i_omp )  THEN
    1396 
    1397           DO  k = nzb+1, nzb_max_l
    1398 
    1399              ibit2 = REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp )
    1400              ibit1 = REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp )
    1401              ibit0 = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )
    1402 
    1403              u_comp                    = u(k,j,i) - u_gtrans + u_stokes_zu(k)
    1404              swap_flux_x_local(k,j,tn) = u_comp * (                            &
    1405                                                ( 37.0_wp * ibit2 * adv_sca_5   &
    1406                                             +     7.0_wp * ibit1 * adv_sca_3   &
    1407                                             +              ibit0 * adv_sca_1   &
    1408                                                ) *                             &
    1409                                             ( sk(k,j,i)   + sk(k,j,i-1)    )   &
    1410                                         -      (  8.0_wp * ibit2 * adv_sca_5   &
    1411                                             +              ibit1 * adv_sca_3   &
    1412                                                ) *                             &
    1413                                             ( sk(k,j,i+1) + sk(k,j,i-2)    )   &
    1414                                         +      (           ibit2 * adv_sca_5   &
    1415                                                ) *                             &
    1416                                             ( sk(k,j,i+2) + sk(k,j,i-3)    )   &
    1417                                                   )
    1418 
    1419              swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                    &
    1420                                                ( 10.0_wp * ibit2 * adv_sca_5   &
    1421                                             +     3.0_wp * ibit1 * adv_sca_3   &
    1422                                             +              ibit0 * adv_sca_1   &
    1423                                                ) *                             &
    1424                                             ( sk(k,j,i)   - sk(k,j,i-1)    )   &
    1425                                         -      (  5.0_wp * ibit2 * adv_sca_5   &
    1426                                             +              ibit1 * adv_sca_3   &
    1427                                                ) *                             &
    1428                                             ( sk(k,j,i+1) - sk(k,j,i-2)    )   &
    1429                                         +      (           ibit2 * adv_sca_5   &
    1430                                                ) *                             &
    1431                                             ( sk(k,j,i+2) - sk(k,j,i-3)    )   &
    1432                                                           )
    1433 
     1788       CASE ( 'q' )
     1789
     1790          DO  k = nzb, nzt
     1791             sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                                        &
     1792                                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )     &
     1793                                                 * ( w(k,j,i) - hom(k,1,3,0)                 )     &
     1794                                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )     &
     1795                                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )     &
     1796                                     ) * weight_substep(intermediate_timestep_count)
    14341797          ENDDO
    14351798
    1436           DO  k = nzb_max_l+1, nzt
    1437 
    1438              u_comp                    = u(k,j,i) - u_gtrans + u_stokes_zu(k)
    1439              swap_flux_x_local(k,j,tn) = u_comp * (                            &
    1440                                       37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) )  &
    1441                                     -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) )  &
    1442                                     +           ( sk(k,j,i+2) + sk(k,j,i-3) )  &
    1443                                                   ) * adv_sca_5
    1444 
    1445              swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                    &
    1446                                       10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) )  &
    1447                                     -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) )  &
    1448                                     +           ( sk(k,j,i+2) - sk(k,j,i-3) )  &
    1449                                                           ) * adv_sca_5
    1450 
     1799       CASE ( 'qc' )
     1800
     1801          DO  k = nzb, nzt
     1802             sums_wsqcs_ws_l(k,tn)  = sums_wsqcs_ws_l(k,tn) +                                      &
     1803                                      ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
     1804                                                  * ( w(k,j,i) - hom(k,1,3,0)                 )    &
     1805                                      + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )    &
     1806                                                  *   ABS( w(k,j,i) - hom(k,1,3,0)            )    &
     1807                                      ) * weight_substep(intermediate_timestep_count)
    14511808          ENDDO
    14521809
    1453        ENDIF
    1454 !
    1455 !--    Now compute the fluxes for the horizontal termns up to the highest
    1456 !--    topography.
    1457        DO  k = nzb+1, nzb_max_l
    1458 
    1459           ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
    1460           ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp )
    1461           ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp )
    1462 
    1463           u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
    1464           flux_r(k) = u_comp * (                                               &
    1465                      ( 37.0_wp * ibit2 * adv_sca_5                             &
    1466                   +     7.0_wp * ibit1 * adv_sca_3                             &
    1467                   +              ibit0 * adv_sca_1                             &
    1468                      ) *                                                       &
    1469                              ( sk(k,j,i+1) + sk(k,j,i)   )                     &
    1470               -      (  8.0_wp * ibit2 * adv_sca_5                             &
    1471                   +              ibit1 * adv_sca_3                             &
    1472                      ) *                                                       &
    1473                              ( sk(k,j,i+2) + sk(k,j,i-1) )                     &
    1474               +      (           ibit2 * adv_sca_5                             &
    1475                      ) *                                                       &
    1476                              ( sk(k,j,i+3) + sk(k,j,i-2) )                     &
    1477                                )
    1478 
    1479           diss_r(k) = -ABS( u_comp ) * (                                       &
    1480                      ( 10.0_wp * ibit2 * adv_sca_5                             &
    1481                   +     3.0_wp * ibit1 * adv_sca_3                             &
    1482                   +              ibit0 * adv_sca_1                             &
    1483                      ) *                                                       &
    1484                              ( sk(k,j,i+1) - sk(k,j,i)  )                      &
    1485               -      (  5.0_wp * ibit2 * adv_sca_5                             &
    1486                   +              ibit1 * adv_sca_3                             &
    1487                      ) *                                                       &
    1488                              ( sk(k,j,i+2) - sk(k,j,i-1) )                     &
    1489               +      (           ibit2 * adv_sca_5                             &
    1490                      ) *                                                       &
    1491                              ( sk(k,j,i+3) - sk(k,j,i-2) )                     &
    1492                                        )
    1493 
    1494           ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp )
    1495           ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp )
    1496           ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp )
    1497 
    1498           v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
    1499           flux_n(k) = v_comp * (                                              &
    1500                      ( 37.0_wp * ibit5 * adv_sca_5                            &
    1501                   +     7.0_wp * ibit4 * adv_sca_3                            &
    1502                   +              ibit3 * adv_sca_1                            &
    1503                      ) *                                                      &
    1504                              ( sk(k,j+1,i) + sk(k,j,i)   )                    &
    1505               -      (  8.0_wp * ibit5 * adv_sca_5                            &
    1506                   +              ibit4 * adv_sca_3                            &
    1507                      ) *                                                      &
    1508                              ( sk(k,j+2,i) + sk(k,j-1,i) )                    &
    1509               +      (           ibit5 * adv_sca_5                            &
    1510                      ) *                                                      &
    1511                              ( sk(k,j+3,i) + sk(k,j-2,i) )                    &
    1512                                )
    1513 
    1514           diss_n(k) = -ABS( v_comp ) * (                                      &
    1515                      ( 10.0_wp * ibit5 * adv_sca_5                            &
    1516                   +     3.0_wp * ibit4 * adv_sca_3                            &
    1517                   +              ibit3 * adv_sca_1                            &
    1518                      ) *                                                      &
    1519                              ( sk(k,j+1,i) - sk(k,j,i)   )                    &
    1520               -      (  5.0_wp * ibit5 * adv_sca_5                            &
    1521                   +              ibit4 * adv_sca_3                            &
    1522                      ) *                                                      &
    1523                              ( sk(k,j+2,i) - sk(k,j-1,i) )                    &
    1524               +      (           ibit5 * adv_sca_5                            &
    1525                      ) *                                                      &
    1526                              ( sk(k,j+3,i) - sk(k,j-2,i) )                    &
    1527                                        )
    1528        ENDDO
    1529 !
    1530 !--    Now compute the fluxes for the horizontal terms above the topography
    1531 !--    where no degradation along the horizontal parts is necessary (except
    1532 !--    for the non-cyclic lateral boundaries treated by nzb_max_l).
    1533        DO  k = nzb_max_l+1, nzt
    1534 
    1535           u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
    1536           flux_r(k) = u_comp * (                                              &
    1537                       37.0_wp * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
    1538                     -  8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
    1539                     +           ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
    1540           diss_r(k) = -ABS( u_comp ) * (                                      &
    1541                       10.0_wp * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
    1542                     -  5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
    1543                     +           ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
    1544 
    1545           v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
    1546           flux_n(k) = v_comp * (                                              &
    1547                       37.0_wp * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
    1548                     -  8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
    1549                     +           ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
    1550           diss_n(k) = -ABS( v_comp ) * (                                      &
    1551                       10.0_wp * ( sk(k,j+1,i) - sk(k,j,i)   )                 &
    1552                     -  5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
    1553                     +           ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
    1554 
    1555        ENDDO
    1556 !
    1557 !--    Now, compute vertical fluxes. Split loop into a part treating the
    1558 !--    lowest grid points with indirect indexing, a main loop without
    1559 !--    indirect indexing, and a loop for the uppermost grip points with
    1560 !--    indirect indexing. This allows better vectorization for the main loop.
    1561 !--    First, compute the flux at model surface, which need has to be
    1562 !--    calculated explicetely for the tendency at
    1563 !--    the first w-level. For topography wall this is done implicitely by
    1564 !--    advc_flag.
    1565        flux_t(nzb) = 0.0_wp
    1566        diss_t(nzb) = 0.0_wp
    1567 
    1568        DO  k = nzb+1, nzb+1
    1569           ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
    1570           ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
    1571           ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
    1572 !
    1573 !--       k index has to be modified near bottom and top, else array
    1574 !--       subscripts will be exceeded.
    1575           k_ppp = k + 3 * ibit8
    1576           k_pp  = k + 2 * ( 1 - ibit6  )
    1577           k_mm  = k - 2 * ibit8
    1578 
    1579           flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
    1580                      ( 37.0_wp * ibit8 * adv_sca_5                            &
    1581                   +     7.0_wp * ibit7 * adv_sca_3                            &
    1582                   +              ibit6 * adv_sca_1                            &
    1583                      ) *                                                      &
    1584                              ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
    1585               -      (  8.0_wp * ibit8 * adv_sca_5                            &
    1586                   +              ibit7 * adv_sca_3                            &
    1587                      ) *                                                      &
    1588                              ( sk(k_pp,j,i) + sk(k-1,j,i)  )                  &
    1589               +      (           ibit8 * adv_sca_5                            &
    1590                      ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                  &
    1591                                                  )
    1592 
    1593           diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
    1594                      ( 10.0_wp * ibit8 * adv_sca_5                            &
    1595                   +     3.0_wp * ibit7 * adv_sca_3                            &
    1596                   +              ibit6 * adv_sca_1                            &
    1597                      ) *                                                      &
    1598                              ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
    1599               -      (  5.0_wp * ibit8 * adv_sca_5                            &
    1600                   +              ibit7 * adv_sca_3                            &
    1601                      ) *                                                      &
    1602                              ( sk(k_pp,j,i)  - sk(k-1,j,i)  )                 &
    1603               +      (           ibit8 * adv_sca_5                            &
    1604                      ) *                                                      &
    1605                              ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    1606                                                          )
    1607        ENDDO
    1608 
    1609        DO  k = nzb+2, nzt-2
    1610           ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
    1611           ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
    1612           ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
    1613 
    1614           flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
    1615                      ( 37.0_wp * ibit8 * adv_sca_5                            &
    1616                   +     7.0_wp * ibit7 * adv_sca_3                            &
    1617                   +              ibit6 * adv_sca_1                            &
    1618                      ) *                                                      &
    1619                              ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
    1620               -      (  8.0_wp * ibit8 * adv_sca_5                            &
    1621                   +              ibit7 * adv_sca_3                            &
    1622                      ) *                                                      &
    1623                              ( sk(k+2,j,i) + sk(k-1,j,i)  )                   &
    1624               +      (           ibit8 * adv_sca_5                            &
    1625                      ) *     ( sk(k+3,j,i)+ sk(k-2,j,i) )                     &
    1626                                                  )
    1627 
    1628           diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
    1629                      ( 10.0_wp * ibit8 * adv_sca_5                            &
    1630                   +     3.0_wp * ibit7 * adv_sca_3                            &
    1631                   +              ibit6 * adv_sca_1                            &
    1632                      ) *                                                      &
    1633                              ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
    1634               -      (  5.0_wp * ibit8 * adv_sca_5                            &
    1635                   +              ibit7 * adv_sca_3                            &
    1636                      ) *                                                      &
    1637                              ( sk(k+2,j,i)  - sk(k-1,j,i)  )                  &
    1638               +      (           ibit8 * adv_sca_5                            &
    1639                      ) *                                                      &
    1640                              ( sk(k+3,j,i) - sk(k-2,j,i) )                    &
    1641                                                          )
    1642        ENDDO
    1643 
    1644        DO  k = nzt-1, nzt-symmetry_flag
    1645           ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
    1646           ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
    1647           ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
    1648 !
    1649 !--       k index has to be modified near bottom and top, else array
    1650 !--       subscripts will be exceeded.
    1651           k_ppp = k + 3 * ibit8
    1652           k_pp  = k + 2 * ( 1 - ibit6  )
    1653           k_mm  = k - 2 * ibit8
    1654 
    1655 
    1656           flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
    1657                      ( 37.0_wp * ibit8 * adv_sca_5                            &
    1658                   +     7.0_wp * ibit7 * adv_sca_3                            &
    1659                   +              ibit6 * adv_sca_1                            &
    1660                      ) *                                                      &
    1661                              ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
    1662               -      (  8.0_wp * ibit8 * adv_sca_5                            &
    1663                   +              ibit7 * adv_sca_3                            &
    1664                      ) *                                                      &
    1665                              ( sk(k_pp,j,i) + sk(k-1,j,i)  )                  &
    1666               +      (           ibit8 * adv_sca_5                            &
    1667                      ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                  &
    1668                                                  )
    1669 
    1670           diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
    1671                      ( 10.0_wp * ibit8 * adv_sca_5                            &
    1672                   +     3.0_wp * ibit7 * adv_sca_3                            &
    1673                   +              ibit6 * adv_sca_1                            &
    1674                      ) *                                                      &
    1675                              ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
    1676               -      (  5.0_wp * ibit8 * adv_sca_5                            &
    1677                   +              ibit7 * adv_sca_3                            &
    1678                      ) *                                                      &
    1679                              ( sk(k_pp,j,i)  - sk(k-1,j,i)  )                 &
    1680               +      (           ibit8 * adv_sca_5                            &
    1681                      ) *                                                      &
    1682                              ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    1683                                                          )
    1684        ENDDO
    1685 
    1686 !
    1687 !--    Set resolved/turbulent flux at model top to zero (w-level). In case that
    1688 !--    a symmetric behavior between bottom and top shall be guaranteed (closed
    1689 !--    channel flow), the flux at nzt is also set to zero.
    1690        IF ( symmetry_flag == 1 ) THEN
    1691           flux_t(nzt) = 0.0_wp
    1692           diss_t(nzt) = 0.0_wp
    1693        ENDIF
    1694        flux_t(nzt+1) = 0.0_wp
    1695        diss_t(nzt+1) = 0.0_wp
    1696 
    1697 
    1698        IF ( limiter )  THEN
    1699 !
    1700 !--       Compute monotone first-order fluxes which are required for mononte
    1701 !--       flux limitation.
    1702           flux_t_1st(nzb) = 0.0_wp
    1703           DO  k = nzb+1, nzb_max_l
    1704              flux_t_1st(k) = ( w(k,j,i)   * ( sk(k+1,j,i)  + sk(k,j,i) )       &
    1705                        -  ABS( w(k,j,i) ) * ( sk(k+1,j,i)  - sk(k,j,i) ) )     &
    1706                            * rho_air_zw(k) * adv_sca_1
    1707 !
    1708 !--          In flux limitation the total flux will be corrected. For the sake
    1709 !--          of cleariness the higher-order advective and disspative fluxes
    1710 !--          will be merged onto flux_t.
    1711              flux_t(k) = flux_t(k) + diss_t(k)
    1712              diss_t(k) = 0.0_wp
     1810
     1811       CASE ( 'qi' )
     1812
     1813          DO  k = nzb, nzt
     1814             sums_wsqis_ws_l(k,tn)  = sums_wsqis_ws_l(k,tn) +                                      &
     1815                                      ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
     1816                                                  * ( w(k,j,i) - hom(k,1,3,0)                 )    &
     1817                                      + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )    &
     1818                                                  *   ABS( w(k,j,i) - hom(k,1,3,0)            )    &
     1819                                      ) * weight_substep(intermediate_timestep_count)
    17131820          ENDDO
    1714 !
    1715 !--       Flux limitation of vertical fluxes according to Skamarock (2006).
    1716 !--       Please note, as flux limitation implies linear dependencies of fluxes,
    1717 !--       flux limitation is only made for the vertical advection term. Limitation
    1718 !--       of the horizontal terms cannot be parallelized.
    1719 !--       Due to the linear dependency, the following loop will not be vectorized.
    1720 !--       Further, note that the flux limiter is only applied within the urban
    1721 !--       layer, i.e up to the topography top.
    1722           DO  k = nzb+1, nzb_max_l
    1723 !
    1724 !--          Compute one-dimensional divergence along the vertical direction,
    1725 !--          which is used to correct the advection discretization. This is
    1726 !--          necessary as in one-dimensional space the advection velocity
    1727 !--          should actually be constant.
    1728              div = ( w(k,j,i)   * rho_air_zw(k)                                &
    1729                    - w(k-1,j,i) * rho_air_zw(k-1)                              &
    1730                    ) * drho_air(k) * ddzw(k)
    1731 !
    1732 !--          Compute monotone solution of the advection equation from
    1733 !--          1st-order fluxes. Please note, the advection equation is corrected
    1734 !--          by the divergence term (in 1D the advective flow should be divergence
    1735 !--          free). Moreover, please note, as time-increment the full timestep
    1736 !--          is used, even though a Runge-Kutta scheme will be used. However,
    1737 !--          the length of the actual time increment is not important at all
    1738 !--          since it cancels out later when the fluxes are limited.
    1739              mon = sk(k,j,i) + ( - ( flux_t_1st(k) - flux_t_1st(k-1) )         &
    1740                              * drho_air(k) * ddzw(k)                           &
    1741                              + div * sk(k,j,i)                                 &
    1742                                ) * dt_3d
    1743 !
    1744 !--          Determine minimum and maximum values along the numerical stencil.
    1745              k_mmm = MAX( k - 3, nzb + 1 )
    1746              k_ppp = MIN( k + 3, nzt + 1 )
    1747 
    1748              min_val = MINVAL( sk(k_mmm:k_ppp,j,i) )
    1749              max_val = MAXVAL( sk(k_mmm:k_ppp,j,i) )
    1750 !
    1751 !--          Compute difference between high- and low-order fluxes, which may
    1752 !--          act as correction fluxes
    1753              f_corr_t = flux_t(k)   - flux_t_1st(k)
    1754              f_corr_d = flux_t(k-1) - flux_t_1st(k-1)
    1755 !
    1756 !--          Determine outgoing fluxes, i.e. the part of the fluxes which can
    1757 !--          decrease the value within the grid box
    1758              f_corr_t_out = MAX( 0.0_wp, f_corr_t )
    1759              f_corr_d_out = MIN( 0.0_wp, f_corr_d )
    1760 !
    1761 !--          Determine ingoing fluxes, i.e. the part of the fluxes which can
    1762 !--          increase the value within the grid box
    1763              f_corr_t_in = MIN( 0.0_wp, f_corr_t)
    1764              f_corr_d_in = MAX( 0.0_wp, f_corr_d)
    1765 !
    1766 !--          Compute divergence of outgoing correction fluxes
    1767              div_out = - ( f_corr_t_out - f_corr_d_out ) * drho_air(k)         &
    1768                                                          * ddzw(k) * dt_3d
    1769 !
    1770 !--          Compute divergence of ingoing correction fluxes
    1771              div_in = - ( f_corr_t_in - f_corr_d_in )    * drho_air(k)         &
    1772                                                          * ddzw(k) * dt_3d
    1773 !
    1774 !--          Check if outgoing fluxes can lead to undershoots, i.e. values smaller
    1775 !--          than the minimum value within the numerical stencil. If so, limit
    1776 !--          them.
    1777              IF ( mon - min_val < - div_out  .AND.  ABS( div_out ) > 0.0_wp )  &
    1778              THEN
    1779                 fac_correction = ( mon - min_val ) / ( - div_out )
    1780                 f_corr_t_out = f_corr_t_out * fac_correction
    1781                 f_corr_d_out = f_corr_d_out * fac_correction
    1782              ENDIF
    1783 !
    1784 !--          Check if ingoing fluxes can lead to overshoots, i.e. values larger
    1785 !--          than the maximum value within the numerical stencil. If so, limit
    1786 !--          them.
    1787              IF ( mon - max_val > - div_in  .AND.  ABS( div_in ) > 0.0_wp )    &
    1788              THEN
    1789                 fac_correction = ( mon - max_val ) / ( - div_in )
    1790                 f_corr_t_in = f_corr_t_in * fac_correction
    1791                 f_corr_d_in = f_corr_d_in * fac_correction
    1792              ENDIF
    1793 !
    1794 !--          Finally add the limited fluxes to the original ones. If no
    1795 !--          flux limitation was done, the fluxes equal the original ones.
    1796              flux_t(k)   = flux_t_1st(k)   + f_corr_t_out + f_corr_t_in
    1797              flux_t(k-1) = flux_t_1st(k-1) + f_corr_d_out + f_corr_d_in
     1821
     1822       CASE ( 'qr' )
     1823
     1824          DO  k = nzb, nzt
     1825             sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn) +                                      &
     1826                                      ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
     1827                                                  * ( w(k,j,i) - hom(k,1,3,0)                 )    &
     1828                                      + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )    &
     1829                                                  *   ABS( w(k,j,i) - hom(k,1,3,0)            )    &
     1830                                      ) * weight_substep(intermediate_timestep_count)
    17981831          ENDDO
    1799        ENDIF
    1800 !
    1801 !--    Now compute the tendency term including divergence correction.
    1802        DO  k = nzb+1, nzb_max_l
    1803 
    1804           flux_d    = flux_t(k-1)
    1805           diss_d    = diss_t(k-1)
    1806 
    1807           ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
    1808           ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp )
    1809           ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp )
    1810 
    1811           ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp )
    1812           ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp )
    1813           ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp )
    1814 
    1815           ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
    1816           ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
    1817           ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
    1818 !
    1819 !--       Calculate the divergence of the velocity field. A respective
    1820 !--       correction is needed to overcome numerical instabilities introduced
    1821 !--       by a not sufficient reduction of divergences near topography.
    1822           div         =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
    1823                           - u(k,j,i)   * (                                     &
    1824                         REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )       &
    1825                       + REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp )       &
    1826                       + REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp )       &
    1827                                          )                                     &
    1828                           ) * ddx                                              &
    1829                         + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
    1830                           - v(k,j,i)   * (                                     &
    1831                         REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp )       &
    1832                       + REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp )       &
    1833                       + REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp )       &
    1834                                          )                                     &
    1835                           ) * ddy                                              &
    1836                         + ( w(k,j,i)   * rho_air_zw(k) *                       &
    1837                                          ( ibit6 + ibit7 + ibit8 )             &
    1838                           - w(k-1,j,i) * rho_air_zw(k-1) *                     &
    1839                                          (                                     &
    1840                         REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp )       &
    1841                       + REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp )       &
    1842                       + REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp )       &
    1843                                          )                                     &
    1844                           ) * drho_air(k) * ddzw(k)
    1845 
    1846           tend(k,j,i) = tend(k,j,i) - (                                        &
    1847                         ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) -  &
    1848                           swap_diss_x_local(k,j,tn)            ) * ddx         &
    1849                       + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   -  &
    1850                           swap_diss_y_local(k,tn)              ) * ddy         &
    1851                       + ( ( flux_t(k) + diss_t(k) ) -                          &
    1852                           ( flux_d    + diss_d    )                            &
    1853                                                     ) * drho_air(k) * ddzw(k)  &
    1854                                       ) + sk(k,j,i) * div
    1855 
    1856 
    1857           swap_flux_y_local(k,tn)   = flux_n(k)
    1858           swap_diss_y_local(k,tn)   = diss_n(k)
    1859           swap_flux_x_local(k,j,tn) = flux_r(k)
    1860           swap_diss_x_local(k,j,tn) = diss_r(k)
    1861 
    1862        ENDDO
    1863 
    1864        DO  k = nzb_max_l+1, nzt
    1865 
    1866           flux_d    = flux_t(k-1)
    1867           diss_d    = diss_t(k-1)
    1868 !
    1869 !--       Calculate the divergence of the velocity field. A respective
    1870 !--       correction is needed to overcome numerical instabilities introduced
    1871 !--       by a not sufficient reduction of divergences near topography.
    1872           div         =   ( u(k,j,i+1) - u(k,j,i) ) * ddx                      &
    1873                         + ( v(k,j+1,i) - v(k,j,i) ) * ddy                      &
    1874                         + ( w(k,j,i)   * rho_air_zw(k)                         &
    1875                           - w(k-1,j,i) * rho_air_zw(k-1)                       &
    1876                                                   ) * drho_air(k) * ddzw(k)
    1877 
    1878           tend(k,j,i) = tend(k,j,i) - (                                        &
    1879                         ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) -  &
    1880                           swap_diss_x_local(k,j,tn)            ) * ddx         &
    1881                       + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   -  &
    1882                           swap_diss_y_local(k,tn)              ) * ddy         &
    1883                       + ( ( flux_t(k) + diss_t(k) ) -                          &
    1884                           ( flux_d    + diss_d    )                            &
    1885                                                     ) * drho_air(k) * ddzw(k)  &
    1886                                       ) + sk(k,j,i) * div
    1887 
    1888 
    1889           swap_flux_y_local(k,tn)   = flux_n(k)
    1890           swap_diss_y_local(k,tn)   = diss_n(k)
    1891           swap_flux_x_local(k,j,tn) = flux_r(k)
    1892           swap_diss_x_local(k,j,tn) = diss_r(k)
    1893 
    1894        ENDDO
    1895 
    1896 !
    1897 !--    Evaluation of statistics.
    1898        SELECT CASE ( sk_char )
    1899 
    1900           CASE ( 'pt' )
    1901 
    1902              DO  k = nzb, nzt
    1903                 sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +                &
    1904                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    1905                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
    1906                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
    1907                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
    1908                     ) * weight_substep(intermediate_timestep_count)
    1909              ENDDO
    1910 
    1911           CASE ( 'sa' )
    1912 
    1913              DO  k = nzb, nzt
    1914                 sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +                &
    1915                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    1916                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
    1917                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
    1918                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
    1919                     ) * weight_substep(intermediate_timestep_count)
    1920              ENDDO
    1921 
    1922           CASE ( 'q' )
    1923 
    1924              DO  k = nzb, nzt
    1925                 sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                 &
    1926                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    1927                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
    1928                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
    1929                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
    1930                     ) * weight_substep(intermediate_timestep_count)
    1931              ENDDO
    1932 
    1933           CASE ( 'qc' )
    1934 
    1935              DO  k = nzb, nzt
    1936                 sums_wsqcs_ws_l(k,tn)  = sums_wsqcs_ws_l(k,tn) +               &
    1937                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    1938                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
    1939                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
    1940                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
    1941                     ) * weight_substep(intermediate_timestep_count)
    1942              ENDDO
    1943 
    1944           CASE ( 'qi' )
    1945 
    1946              DO  k = nzb, nzt
    1947                 sums_wsqis_ws_l(k,tn)  = sums_wsqis_ws_l(k,tn) +               &
    1948                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    1949                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
    1950                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
    1951                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
    1952                     ) * weight_substep(intermediate_timestep_count)
    1953              ENDDO
    1954 
    1955           CASE ( 'qr' )
    1956 
    1957              DO  k = nzb, nzt
    1958                 sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn) +               &
    1959                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    1960                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
    1961                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
    1962                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
    1963                     ) * weight_substep(intermediate_timestep_count)
    1964              ENDDO
    1965 
    1966           CASE ( 'nc' )
    1967 
    1968              DO  k = nzb, nzt
    1969                 sums_wsncs_ws_l(k,tn)  = sums_wsncs_ws_l(k,tn) +               &
    1970                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    1971                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
    1972                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
    1973                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
    1974                     ) * weight_substep(intermediate_timestep_count)
    1975              ENDDO
    1976 
    1977           CASE ( 'ni' )
    1978 
    1979              DO  k = nzb, nzt
    1980                 sums_wsnis_ws_l(k,tn)  = sums_wsnis_ws_l(k,tn) +               &
    1981                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    1982                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
    1983                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
    1984                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
    1985                     ) * weight_substep(intermediate_timestep_count)
    1986              ENDDO
    1987 
    1988           CASE ( 'nr' )
    1989 
    1990              DO  k = nzb, nzt
    1991                 sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +               &
    1992                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    1993                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
    1994                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
    1995                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
    1996                     ) * weight_substep(intermediate_timestep_count)
    1997              ENDDO
    1998 
    1999           CASE ( 's' )
    2000 
    2001              DO  k = nzb, nzt
    2002                 sums_wsss_ws_l(k,tn)  = sums_wsss_ws_l(k,tn) +                 &
    2003                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    2004                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
    2005                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
    2006                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
    2007                     ) * weight_substep(intermediate_timestep_count)
    2008              ENDDO
    2009 
    2010          CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' )
    2011 
    2012              DO  k = nzb, nzt
    2013                 sums_salsa_ws_l(k,tn)  = sums_salsa_ws_l(k,tn) +               &
    2014                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    2015                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
    2016                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
    2017                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
    2018                     ) * weight_substep(intermediate_timestep_count)
    2019              ENDDO
    2020 
    2021 !          CASE ( 'kc' )
    2022           !kk Has to be implemented for kpp chemistry
    2023 
    2024          END SELECT
    2025 
    2026     END SUBROUTINE advec_s_ws_ij
    2027 
    2028 
    2029 
    2030 
    2031 !------------------------------------------------------------------------------!
     1832
     1833       CASE ( 'nc' )
     1834
     1835          DO  k = nzb, nzt
     1836             sums_wsncs_ws_l(k,tn)  = sums_wsncs_ws_l(k,tn) +                                      &
     1837                                      ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
     1838                                                  * ( w(k,j,i) - hom(k,1,3,0)                 )    &
     1839                                      + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )    &
     1840                                                  *   ABS( w(k,j,i) - hom(k,1,3,0)            )    &
     1841                                      ) * weight_substep(intermediate_timestep_count)
     1842          ENDDO
     1843
     1844       CASE ( 'ni' )
     1845
     1846          DO  k = nzb, nzt
     1847             sums_wsnis_ws_l(k,tn)  = sums_wsnis_ws_l(k,tn) +                                      &
     1848                                      ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
     1849                                                  * ( w(k,j,i) - hom(k,1,3,0)                 )    &
     1850                                      + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )    &
     1851                                                  *   ABS( w(k,j,i) - hom(k,1,3,0)            )    &
     1852                                      ) * weight_substep(intermediate_timestep_count)
     1853          ENDDO
     1854
     1855       CASE ( 'nr' )
     1856
     1857          DO  k = nzb, nzt
     1858             sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +                                      &
     1859                                      ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
     1860                                                  * ( w(k,j,i) - hom(k,1,3,0)                 )    &
     1861                                      + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )    &
     1862                                                  *   ABS( w(k,j,i) - hom(k,1,3,0)            )    &
     1863                                      ) * weight_substep(intermediate_timestep_count)
     1864          ENDDO
     1865
     1866       CASE ( 's' )
     1867
     1868          DO  k = nzb, nzt
     1869             sums_wsss_ws_l(k,tn)  = sums_wsss_ws_l(k,tn) +                                        &
     1870                                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )     &
     1871                                                 * ( w(k,j,i) - hom(k,1,3,0)                 )     &
     1872                                     + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )     &
     1873                                                 *   ABS( w(k,j,i) - hom(k,1,3,0)            )     &
     1874                                     ) * weight_substep(intermediate_timestep_count)
     1875          ENDDO
     1876
     1877      CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' )
     1878
     1879          DO  k = nzb, nzt
     1880             sums_salsa_ws_l(k,tn)  = sums_salsa_ws_l(k,tn) +                                      &
     1881                                      ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
     1882                                                  * ( w(k,j,i) - hom(k,1,3,0)                 )    &
     1883                                      + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )    &
     1884                                                  *   ABS( w(k,j,i) - hom(k,1,3,0)            )    &
     1885                                      ) * weight_substep(intermediate_timestep_count)
     1886          ENDDO
     1887
     1888!       CASE ( 'kc' )
     1889        !kk Has to be implemented for kpp chemistry
     1890
     1891    END SELECT
     1892
     1893 END SUBROUTINE advec_s_ws_ij
     1894
     1895
     1896
     1897
     1898!--------------------------------------------------------------------------------------------------!
    20321899! Description:
    20331900! ------------
    20341901!> Advection of u-component - Call for grid point i,j
    2035 !------------------------------------------------------------------------------!
    2036     SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
    2037 
    2038 
    2039        INTEGER(iwp) ::  i         !< grid index along x-direction
    2040        INTEGER(iwp) ::  i_omp     !< leftmost index on subdomain, or in case of OpenMP, on thread
    2041        INTEGER(iwp) ::  j         !< grid index along y-direction
    2042        INTEGER(iwp) ::  k         !< grid index along z-direction
    2043        INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
    2044        INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
    2045        INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
    2046        INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
    2047        INTEGER(iwp) ::  tn        !< number of OpenMP thread
    2048 
    2049        REAL(wp)    ::  ibit0   !< flag indicating 1st-order scheme along x-direction
    2050        REAL(wp)    ::  ibit1   !< flag indicating 3rd-order scheme along x-direction
    2051        REAL(wp)    ::  ibit2   !< flag indicating 5th-order scheme along x-direction
    2052        REAL(wp)    ::  ibit3   !< flag indicating 1st-order scheme along y-direction
    2053        REAL(wp)    ::  ibit4   !< flag indicating 3rd-order scheme along y-direction
    2054        REAL(wp)    ::  ibit5   !< flag indicating 5th-order scheme along y-direction
    2055        REAL(wp)    ::  ibit6   !< flag indicating 1st-order scheme along z-direction
    2056        REAL(wp)    ::  ibit7   !< flag indicating 3rd-order scheme along z-direction
    2057        REAL(wp)    ::  ibit8   !< flag indicating 5th-order scheme along z-direction
    2058        REAL(wp)    ::  diss_d   !< artificial dissipation term at grid box bottom
    2059        REAL(wp)    ::  div      !< diverence on u-grid
    2060        REAL(wp)    ::  flux_d   !< 6th-order flux at grid box bottom
    2061        REAL(wp)    ::  gu       !< Galilei-transformation velocity along x
    2062        REAL(wp)    ::  gv       !< Galilei-transformation velocity along y
    2063        REAL(wp)    ::  u_comp_l !< advection velocity along x at leftmost grid point on subdomain
    2064 
    2065        REAL(wp), DIMENSION(nzb:nzt+1) ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
    2066        REAL(wp), DIMENSION(nzb:nzt+1) ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
    2067        REAL(wp), DIMENSION(nzb:nzt+1) ::  diss_t !< discretized artificial dissipation at top of the grid box
    2068        REAL(wp), DIMENSION(nzb:nzt+1) ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
    2069        REAL(wp), DIMENSION(nzb:nzt+1) ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
    2070        REAL(wp), DIMENSION(nzb:nzt+1) ::  flux_t !< discretized 6th-order flux at top of the grid box
    2071        REAL(wp), DIMENSION(nzb:nzt+1) ::  u_comp !< advection velocity along x
    2072        REAL(wp), DIMENSION(nzb:nzt+1) ::  v_comp !< advection velocity along y
    2073        REAL(wp), DIMENSION(nzb:nzt+1) ::  w_comp !< advection velocity along z
    2074 !
    2075 !--    Used local modified copy of nzb_max (used to degrade order of
    2076 !--    discretization) at non-cyclic boundaries. Modify only at relevant points
    2077 !--    instead of the entire subdomain. This should lead to better
    2078 !--    load balance between boundary and non-boundary PEs.
    2079        IF( ( bc_dirichlet_l  .OR.  bc_radiation_l )  .AND.  i <= nxl + 2  .OR. &
    2080            ( bc_dirichlet_r  .OR.  bc_radiation_r )  .AND.  i >= nxr - 2  .OR. &
    2081            ( bc_dirichlet_s  .OR.  bc_radiation_s )  .AND.  j <= nys + 2  .OR. &
    2082            ( bc_dirichlet_n  .OR.  bc_radiation_n )  .AND.  j >= nyn - 2 )  THEN
    2083           nzb_max_l = nzt
    2084        ELSE
    2085           nzb_max_l = nzb_max
    2086        END IF
    2087 
    2088        gu = 2.0_wp * u_gtrans
    2089        gv = 2.0_wp * v_gtrans
    2090 !
    2091 !--    Compute southside fluxes for the respective boundary of PE
    2092        IF ( j == nys  )  THEN
    2093 
     1902!--------------------------------------------------------------------------------------------------!
     1903 SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
     1904
     1905
     1906    INTEGER(iwp) ::  i         !< grid index along x-direction
     1907    INTEGER(iwp) ::  i_omp     !< leftmost index on subdomain, or in case of OpenMP, on thread
     1908    INTEGER(iwp) ::  j         !< grid index along y-direction
     1909    INTEGER(iwp) ::  k         !< grid index along z-direction
     1910    INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
     1911    INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
     1912    INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
     1913    INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
     1914    INTEGER(iwp) ::  tn        !< number of OpenMP thread
     1915
     1916    REAL(wp)    ::  diss_d   !< artificial dissipation term at grid box bottom
     1917    REAL(wp)    ::  div      !< diverence on u-grid
     1918    REAL(wp)    ::  flux_d   !< 6th-order flux at grid box bottom
     1919    REAL(wp)    ::  gu       !< Galilei-transformation velocity along x
     1920    REAL(wp)    ::  gv       !< Galilei-transformation velocity along y
     1921    REAL(wp)    ::  ibit0   !< flag indicating 1st-order scheme along x-direction
     1922    REAL(wp)    ::  ibit1   !< flag indicating 3rd-order scheme along x-direction
     1923    REAL(wp)    ::  ibit2   !< flag indicating 5th-order scheme along x-direction
     1924    REAL(wp)    ::  ibit3   !< flag indicating 1st-order scheme along y-direction
     1925    REAL(wp)    ::  ibit4   !< flag indicating 3rd-order scheme along y-direction
     1926    REAL(wp)    ::  ibit5   !< flag indicating 5th-order scheme along y-direction
     1927    REAL(wp)    ::  ibit6   !< flag indicating 1st-order scheme along z-direction
     1928    REAL(wp)    ::  ibit7   !< flag indicating 3rd-order scheme along z-direction
     1929    REAL(wp)    ::  ibit8   !< flag indicating 5th-order scheme along z-direction
     1930    REAL(wp)    ::  u_comp_l !< advection velocity along x at leftmost grid point on subdomain
     1931
     1932    REAL(wp), DIMENSION(nzb:nzt+1) ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
     1933    REAL(wp), DIMENSION(nzb:nzt+1) ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
     1934    REAL(wp), DIMENSION(nzb:nzt+1) ::  diss_t !< discretized artificial dissipation at top of the grid box
     1935    REAL(wp), DIMENSION(nzb:nzt+1) ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
     1936    REAL(wp), DIMENSION(nzb:nzt+1) ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
     1937    REAL(wp), DIMENSION(nzb:nzt+1) ::  flux_t !< discretized 6th-order flux at top of the grid box
     1938    REAL(wp), DIMENSION(nzb:nzt+1) ::  u_comp !< advection velocity along x
     1939    REAL(wp), DIMENSION(nzb:nzt+1) ::  v_comp !< advection velocity along y
     1940    REAL(wp), DIMENSION(nzb:nzt+1) ::  w_comp !< advection velocity along z
     1941!
     1942!-- Used local modified copy of nzb_max (used to degrade order of discretization) at non-cyclic
     1943!-- boundaries. Modify only at relevant points instead of the entire subdomain. This should lead to
     1944!-- better load balance between boundary and non-boundary PEs.
     1945    IF( ( bc_dirichlet_l  .OR.  bc_radiation_l )  .AND.  i <= nxl + 2  .OR.                        &
     1946        ( bc_dirichlet_r  .OR.  bc_radiation_r )  .AND.  i >= nxr - 2  .OR.                        &
     1947        ( bc_dirichlet_s  .OR.  bc_radiation_s )  .AND.  j <= nys + 2  .OR.                        &
     1948        ( bc_dirichlet_n  .OR.  bc_radiation_n )  .AND.  j >= nyn - 2 )  THEN
     1949       nzb_max_l = nzt
     1950    ELSE
     1951       nzb_max_l = nzb_max
     1952    END IF
     1953
     1954    gu = 2.0_wp * u_gtrans
     1955    gv = 2.0_wp * v_gtrans
     1956!
     1957!-- Compute southside fluxes for the respective boundary of PE
     1958    IF ( j == nys  )  THEN
     1959
     1960       DO  k = nzb+1, nzb_max_l
     1961
     1962          ibit5 = REAL( IBITS(advc_flags_m(k,j-1,i),5,1), KIND = wp )
     1963          ibit4 = REAL( IBITS(advc_flags_m(k,j-1,i),4,1), KIND = wp )
     1964          ibit3 = REAL( IBITS(advc_flags_m(k,j-1,i),3,1), KIND = wp )
     1965
     1966          v_comp(k)      = v(k,j,i) + v(k,j,i-1) - gv
     1967          flux_s_u(k,tn) = v_comp(k) * (                                                           &
     1968                           (  37.0_wp * ibit5 * adv_mom_5                                          &
     1969                           +   7.0_wp * ibit4 * adv_mom_3                                          &
     1970                           +            ibit3 * adv_mom_1 ) * ( u(k,j,i)    + u(k,j-1,i) )         &
     1971                           - ( 8.0_wp * ibit5 * adv_mom_5                                          &
     1972                           +            ibit4 * adv_mom_3 ) * ( u(k,j+1,i)  + u(k,j-2,i) )         &
     1973                           + (          ibit5 * adv_mom_5 ) *  ( u(k,j+2,i) + u(k,j-3,i) )         &
     1974                                       )
     1975
     1976          diss_s_u(k,tn) = - ABS ( v_comp(k) ) * (                                                 &
     1977                           (  10.0_wp * ibit5 * adv_mom_5                                          &
     1978                           +   3.0_wp * ibit4 * adv_mom_3                                          &
     1979                           +            ibit3 * adv_mom_1 ) * ( u(k,j,i)   - u(k,j-1,i) )          &
     1980                           - ( 5.0_wp * ibit5 * adv_mom_5                                          &
     1981                           +            ibit4 * adv_mom_3 ) * ( u(k,j+1,i) - u(k,j-2,i) )          &
     1982                           + (          ibit5 * adv_mom_5 ) * ( u(k,j+2,i) - u(k,j-3,i) )          &
     1983                                                )
     1984
     1985       ENDDO
     1986
     1987       DO  k = nzb_max_l+1, nzt
     1988
     1989          v_comp(k)      = v(k,j,i) + v(k,j,i-1) - gv
     1990          flux_s_u(k,tn) = v_comp(k) * (                                                           &
     1991                             37.0_wp * ( u(k,j,i)   + u(k,j-1,i) )                                 &
     1992                           -  8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) )                                 &
     1993                           +           ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
     1994          diss_s_u(k,tn) = - ABS(v_comp(k)) * (                                                    &
     1995                                    10.0_wp * ( u(k,j,i)   - u(k,j-1,i) )                          &
     1996                                  -  5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) )                          &
     1997                                  +           ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
     1998
     1999       ENDDO
     2000
     2001    ENDIF
     2002!
     2003!-- Compute leftside fluxes for the respective boundary of PE
     2004    IF ( i == i_omp  .OR.  i == nxlu )  THEN
     2005
     2006       DO  k = nzb+1, nzb_max_l
     2007
     2008          ibit2 = REAL( IBITS(advc_flags_m(k,j,i-1),2,1), KIND = wp )
     2009          ibit1 = REAL( IBITS(advc_flags_m(k,j,i-1),1,1), KIND = wp )
     2010          ibit0 = REAL( IBITS(advc_flags_m(k,j,i-1),0,1), KIND = wp )
     2011
     2012          u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
     2013          flux_l_u(k,j,tn) = u_comp_l * (                                                          &
     2014                             (  37.0_wp * ibit2 * adv_mom_5                                        &
     2015                             +   7.0_wp * ibit1 * adv_mom_3                                        &
     2016                             +            ibit0 * adv_mom_1 ) * ( u(k,j,i)   + u(k,j,i-1) )        &
     2017                             - ( 8.0_wp * ibit2 * adv_mom_5                                        &
     2018                             +            ibit1 * adv_mom_3 ) * ( u(k,j,i+1) + u(k,j,i-2) )        &
     2019                             + (          ibit2 * adv_mom_5 ) * ( u(k,j,i+2) + u(k,j,i-3) )        &
     2020                                        )
     2021
     2022          diss_l_u(k,j,tn) = - ABS( u_comp_l ) * (                                                 &
     2023                             (  10.0_wp * ibit2 * adv_mom_5                                        &
     2024                             +   3.0_wp * ibit1 * adv_mom_3                                        &
     2025                             +           ibit0  * adv_mom_1 ) * ( u(k,j,i)   - u(k,j,i-1) )        &
     2026                             - ( 5.0_wp * ibit2 * adv_mom_5                                        &
     2027                             +            ibit1 * adv_mom_3 ) * ( u(k,j,i+1) - u(k,j,i-2) )        &
     2028                             + (          ibit2 * adv_mom_5 ) * ( u(k,j,i+2) - u(k,j,i-3) )        &
     2029                                                 )
     2030
     2031       ENDDO
     2032
     2033       DO  k = nzb_max_l+1, nzt
     2034
     2035          u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
     2036          flux_l_u(k,j,tn) = u_comp_l * (                                                          &
     2037                               37.0_wp * ( u(k,j,i)   + u(k,j,i-1) )                               &
     2038                             -  8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) )                               &
     2039                             +           ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
     2040          diss_l_u(k,j,tn) = - ABS(u_comp_l) * (                                                   &
     2041                               10.0_wp * ( u(k,j,i)   - u(k,j,i-1) )                               &
     2042                             -  5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) )                               &
     2043                             +           ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
     2044
     2045       ENDDO
     2046
     2047    ENDIF
     2048!
     2049!-- Now compute the fluxes tendency terms for the horizontal and vertical parts.
     2050    DO  k = nzb+1, nzb_max_l
     2051
     2052       ibit2 = REAL( IBITS(advc_flags_m(k,j,i),2,1), KIND = wp )
     2053       ibit1 = REAL( IBITS(advc_flags_m(k,j,i),1,1), KIND = wp )
     2054       ibit0 = REAL( IBITS(advc_flags_m(k,j,i),0,1), KIND = wp )
     2055
     2056       u_comp(k) = u(k,j,i+1) + u(k,j,i)
     2057       flux_r(k) = ( u_comp(k) - gu ) * (                                                          &
     2058                   (  37.0_wp * ibit2 * adv_mom_5                                                  &
     2059                   +   7.0_wp * ibit1 * adv_mom_3                                                  &
     2060                   +            ibit0 * adv_mom_1 ) * ( u(k,j,i+1) + u(k,j,i)   )                  &
     2061                   - ( 8.0_wp * ibit2 * adv_mom_5                                                  &
     2062                   +            ibit1 * adv_mom_3 ) * ( u(k,j,i+2) + u(k,j,i-1) )                  &
     2063                   + (          ibit2 * adv_mom_5 ) * ( u(k,j,i+3) + u(k,j,i-2) )                  &
     2064                                        )
     2065
     2066       diss_r(k) = - ABS( u_comp(k) - gu ) * (                                                     &
     2067                   (  10.0_wp * ibit2 * adv_mom_5                                                  &
     2068                   +   3.0_wp * ibit1 * adv_mom_3                                                  &
     2069                   +            ibit0 * adv_mom_1 ) * ( u(k,j,i+1) - u(k,j,i)   )                  &
     2070                   - ( 5.0_wp * ibit2 * adv_mom_5                                                  &
     2071                   +            ibit1 * adv_mom_3 ) * ( u(k,j,i+2) - u(k,j,i-1) )                  &
     2072                   + (          ibit2 * adv_mom_5 ) * ( u(k,j,i+3) - u(k,j,i-2) )                  &
     2073                                             )
     2074
     2075       ibit5 = REAL( IBITS(advc_flags_m(k,j,i),5,1), KIND = wp )
     2076       ibit4 = REAL( IBITS(advc_flags_m(k,j,i),4,1), KIND = wp )
     2077       ibit3 = REAL( IBITS(advc_flags_m(k,j,i),3,1), KIND = wp )
     2078
     2079       v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
     2080       flux_n(k) = v_comp(k) * (                                                                   &
     2081                   (  37.0_wp * ibit5 * adv_mom_5                                                  &
     2082                   +   7.0_wp * ibit4 * adv_mom_3                                                  &
     2083                   +            ibit3 * adv_mom_1 ) * ( u(k,j+1,i) + u(k,j,i)   )                  &
     2084                   - ( 8.0_wp * ibit5 * adv_mom_5                                                  &
     2085                   +            ibit4 * adv_mom_3 ) * ( u(k,j+2,i) + u(k,j-1,i) )                  &
     2086                   + (          ibit5 * adv_mom_5 ) * ( u(k,j+3,i) + u(k,j-2,i) )                  &
     2087                               )
     2088
     2089       diss_n(k) = - ABS ( v_comp(k) ) * (                                                         &
     2090                   (  10.0_wp * ibit5 * adv_mom_5                                                  &
     2091                   +   3.0_wp * ibit4 * adv_mom_3                                                  &
     2092                   +    ibit3 * adv_mom_1 ) * ( u(k,j+1,i) - u(k,j,i)   )                          &
     2093                   - ( 5.0_wp * ibit5 * adv_mom_5                                                  &
     2094                   +    ibit4 * adv_mom_3 ) * ( u(k,j+2,i) - u(k,j-1,i) )                          &
     2095                   + (  ibit5 * adv_mom_5 ) * ( u(k,j+3,i) - u(k,j-2,i) )                          &
     2096                                         )
     2097    ENDDO
     2098
     2099    DO  k = nzb_max_l+1, nzt
     2100
     2101       u_comp(k) = u(k,j,i+1) + u(k,j,i)
     2102       flux_r(k) = ( u_comp(k) - gu ) * (                                                          &
     2103                      37.0_wp * ( u(k,j,i+1) + u(k,j,i)   )                                        &
     2104                    -  8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) )                                        &
     2105                    +           ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
     2106       diss_r(k) = - ABS( u_comp(k) - gu ) * (                                                     &
     2107                      10.0_wp * ( u(k,j,i+1) - u(k,j,i)   )                                        &
     2108                    -  5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) )                                        &
     2109                    +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
     2110
     2111       v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
     2112       flux_n(k) = v_comp(k) * (                                                                   &
     2113                      37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                                        &
     2114                    -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                                        &
     2115                    +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
     2116       diss_n(k) = - ABS( v_comp(k) ) * (                                                          &
     2117                      10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                                        &
     2118                    -  5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) )                                        &
     2119                    +           ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
     2120
     2121    ENDDO
     2122!
     2123!-- Now, compute vertical fluxes. Split loop into a part treating the lowest grid points with
     2124!-- indirect indexing, a main loop without indirect indexing, and a loop for the uppermost grid
     2125!-- points with indirect indexing. This allows better vectorization for the main loop.
     2126!-- First, compute the flux at model surface, which need has to be calculated explicitly for the
     2127!-- tendency at the first w-level. For topography wall this is done implicitely by advc_flags_m.
     2128    flux_t(nzb) = 0.0_wp
     2129    diss_t(nzb) = 0.0_wp
     2130    w_comp(nzb) = 0.0_wp
     2131
     2132    DO  k = nzb+1, nzb+1
     2133!
     2134!--    k index has to be modified near bottom and top, else array subscripts will be exceeded.
     2135       ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
     2136       ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
     2137       ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp )
     2138
     2139       k_ppp = k + 3 * ibit8
     2140       k_pp  = k + 2 * ( 1 - ibit6 )
     2141       k_mm  = k - 2 * ibit8
     2142
     2143       w_comp(k) = w(k,j,i) + w(k,j,i-1)
     2144       flux_t(k) = w_comp(k) * rho_air_zw(k) * (                                                   &
     2145                   (  37.0_wp * ibit8 * adv_mom_5                                                  &
     2146                   +   7.0_wp * ibit7 * adv_mom_3                                                  &
     2147                   +            ibit6 * adv_mom_1 ) * ( u(k+1,j,i)   + u(k,j,i)    )               &
     2148                   - ( 8.0_wp * ibit8 * adv_mom_5                                                  &
     2149                   +            ibit7 * adv_mom_3 ) * ( u(k_pp,j,i)  + u(k-1,j,i)  )               &
     2150                   + (          ibit8 * adv_mom_5 ) * ( u(k_ppp,j,i) + u(k_mm,j,i) )               &
     2151                                               )
     2152
     2153       diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                                          &
     2154                   (  10.0_wp * ibit8 * adv_mom_5                                                  &
     2155                   +   3.0_wp * ibit7 * adv_mom_3                                                  &
     2156                   +            ibit6 * adv_mom_1 ) * ( u(k+1,j,i)   - u(k,j,i)    )               &
     2157                   - ( 5.0_wp * ibit8 * adv_mom_5                                                  &
     2158                   +            ibit7 * adv_mom_3 ) * ( u(k_pp,j,i)  - u(k-1,j,i)  )               &
     2159                   + (          ibit8 * adv_mom_5 ) * ( u(k_ppp,j,i) - u(k_mm,j,i) )               &
     2160                                                        )
     2161    ENDDO
     2162
     2163    DO  k = nzb+2, nzt-2
     2164
     2165       ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
     2166       ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
     2167       ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp )
     2168
     2169       w_comp(k) = w(k,j,i) + w(k,j,i-1)
     2170       flux_t(k) = w_comp(k) * rho_air_zw(k) * (                                                   &
     2171                   (  37.0_wp * ibit8 * adv_mom_5                                                  &
     2172                   +   7.0_wp * ibit7 * adv_mom_3                                                  &
     2173                   +            ibit6 * adv_mom_1 ) * ( u(k+1,j,i)  + u(k,j,i)  )                  &
     2174                   - ( 8.0_wp * ibit8 * adv_mom_5                                                  &
     2175                   +            ibit7 * adv_mom_3 ) * ( u(k+2,j,i) + u(k-1,j,i) )                  &
     2176                   + (          ibit8 * adv_mom_5 ) * ( u(k+3,j,i) + u(k-2,j,i) )                  &
     2177                                               )
     2178
     2179       diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                                          &
     2180                   (  10.0_wp * ibit8 * adv_mom_5                                                  &
     2181                   +   3.0_wp * ibit7 * adv_mom_3                                                  &
     2182                   +            ibit6 * adv_mom_1 ) * ( u(k+1,j,i)  - u(k,j,i)   )                 &
     2183                   - ( 5.0_wp * ibit8 * adv_mom_5                                                  &
     2184                   +            ibit7 * adv_mom_3 ) * ( u(k+2,j,i)  - u(k-1,j,i) )                 &
     2185                   + (          ibit8 * adv_mom_5 ) * ( u(k+3,j,i) - u(k-2,j,i)  )                 &
     2186                                                        )
     2187    ENDDO
     2188
     2189    DO  k = nzt-1, nzt-symmetry_flag
     2190!
     2191!--    k index has to be modified near bottom and top, else array subscripts will be exceeded.
     2192       ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
     2193       ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
     2194       ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp )
     2195
     2196       k_ppp = k + 3 * ibit8
     2197       k_pp  = k + 2 * ( 1 - ibit6 )
     2198       k_mm  = k - 2 * ibit8
     2199
     2200       w_comp(k) = w(k,j,i) + w(k,j,i-1)
     2201       flux_t(k) = w_comp(k) * rho_air_zw(k) * (                                                   &
     2202                   (  37.0_wp * ibit8 * adv_mom_5                                                  &
     2203                   +   7.0_wp * ibit7 * adv_mom_3                                                  &
     2204                   +            ibit6 * adv_mom_1 ) * ( u(k+1,j,i)   + u(k,j,i)    )               &
     2205                   - ( 8.0_wp * ibit8 * adv_mom_5                                                  &
     2206                   +            ibit7 * adv_mom_3 ) * ( u(k_pp,j,i)  + u(k-1,j,i)  )               &
     2207                   + (          ibit8 * adv_mom_5 ) * ( u(k_ppp,j,i) + u(k_mm,j,i) )               &
     2208                                               )
     2209
     2210       diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                                          &
     2211                   (  10.0_wp * ibit8 * adv_mom_5                                                  &
     2212                   +   3.0_wp * ibit7 * adv_mom_3                                                  &
     2213                   +            ibit6 * adv_mom_1 ) * ( u(k+1,j,i)   - u(k,j,i)    )               &
     2214                   - ( 5.0_wp * ibit8 * adv_mom_5                                                  &
     2215                   +            ibit7 * adv_mom_3 ) * ( u(k_pp,j,i)  - u(k-1,j,i)  )               &
     2216                   + (          ibit8 * adv_mom_5 ) * ( u(k_ppp,j,i) - u(k_mm,j,i) )               &
     2217                                                        )
     2218    ENDDO
     2219
     2220!
     2221!-- Set resolved/turbulent flux at model top to zero (w-level). In case that a symmetric behavior
     2222!-- between bottom and top shall be guaranteed (closed channel flow), the flux at nzt is also set to
     2223!-- zero.
     2224    IF ( symmetry_flag == 1 )  THEN
     2225       flux_t(nzt) = 0.0_wp
     2226       diss_t(nzt) = 0.0_wp
     2227       w_comp(nzt) = 0.0_wp
     2228    ENDIF
     2229    flux_t(nzt+1) = 0.0_wp
     2230    diss_t(nzt+1) = 0.0_wp
     2231    w_comp(nzt+1) = 0.0_wp
     2232
     2233    DO  k = nzb+1, nzb_max_l
     2234
     2235       flux_d    = flux_t(k-1)
     2236       diss_d    = diss_t(k-1)
     2237
     2238       ibit2 = REAL( IBITS(advc_flags_m(k,j,i),2,1), KIND = wp )
     2239       ibit1 = REAL( IBITS(advc_flags_m(k,j,i),1,1), KIND = wp )
     2240       ibit0 = REAL( IBITS(advc_flags_m(k,j,i),0,1), KIND = wp )
     2241
     2242       ibit5 = REAL( IBITS(advc_flags_m(k,j,i),5,1), KIND = wp )
     2243       ibit4 = REAL( IBITS(advc_flags_m(k,j,i),4,1), KIND = wp )
     2244       ibit3 = REAL( IBITS(advc_flags_m(k,j,i),3,1), KIND = wp )
     2245
     2246       ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
     2247       ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
     2248       ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp )
     2249!
     2250!--    Calculate the divergence of the velocity field. A respective correction is needed to overcome
     2251!--    numerical instabilities introduced by an insufficient reduction of divergences near
     2252!--    topography.
     2253       div = ( ( u_comp(k)       * ( ibit0 + ibit1 + ibit2 )                                       &
     2254             - ( u(k,j,i)   + u(k,j,i-1)   )                                                       &
     2255                                 * (                                                               &
     2256                  REAL( IBITS(advc_flags_m(k,j,i-1),0,1), KIND = wp )                              &
     2257                + REAL( IBITS(advc_flags_m(k,j,i-1),1,1), KIND = wp )                              &
     2258                + REAL( IBITS(advc_flags_m(k,j,i-1),2,1), KIND = wp )                              &
     2259                                   )                                                               &
     2260               ) * ddx                                                                             &
     2261             + ( ( v_comp(k) + gv ) * ( ibit3 + ibit4 + ibit5 )                                    &
     2262               - ( v(k,j,i)   + v(k,j,i-1 )  )                                                     &
     2263                                 * (                                                               &
     2264                  REAL( IBITS(advc_flags_m(k,j-1,i),3,1), KIND = wp )                              &
     2265                + REAL( IBITS(advc_flags_m(k,j-1,i),4,1), KIND = wp )                              &
     2266                + REAL( IBITS(advc_flags_m(k,j-1,i),5,1), KIND = wp )                              &
     2267                                   )                                                               &
     2268               ) * ddy                                                                             &
     2269             + ( w_comp(k)   * rho_air_zw(k) * ( ibit6 + ibit7 + ibit8 )                           &
     2270             -   w_comp(k-1) * rho_air_zw(k-1)                                                     &
     2271                                 * (                                                               &
     2272                  REAL( IBITS(advc_flags_m(k-1,j,i),6,1), KIND = wp )                              &
     2273                + REAL( IBITS(advc_flags_m(k-1,j,i),7,1), KIND = wp )                              &
     2274                + REAL( IBITS(advc_flags_m(k-1,j,i),8,1), KIND = wp )                              &
     2275                                   )                                                               &
     2276               ) * drho_air(k) * ddzw(k)                                                           &
     2277             ) * 0.5_wp
     2278
     2279       tend(k,j,i) = tend(k,j,i) - (                                                               &
     2280                         ( flux_r(k) + diss_r(k)                                                   &
     2281                       -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx                             &
     2282                       + ( flux_n(k) + diss_n(k)                                                   &
     2283                       -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy                             &
     2284                       + ( ( flux_t(k) + diss_t(k) )                                               &
     2285                       -   ( flux_d    + diss_d    )                                               &
     2286                         )                         * drho_air(k) * ddzw(k)                         &
     2287                                   ) + div * u(k,j,i)
     2288
     2289       flux_l_u(k,j,tn) = flux_r(k)
     2290       diss_l_u(k,j,tn) = diss_r(k)
     2291       flux_s_u(k,tn)   = flux_n(k)
     2292       diss_s_u(k,tn)   = diss_n(k)
     2293!
     2294!--    Statistical Evaluation of u'u'. The factor has to be applied for right evaluation when
     2295!--    gallilei_trans = .T. .
     2296       sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) +                                                 &
     2297                             ( flux_r(k)                                                           &
     2298                               * ( u_comp(k) - 2.0_wp * hom(k,1,1,0)                   )           &
     2299                               / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )           &
     2300                             + diss_r(k)                                                           &
     2301                               *   ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0)              )           &
     2302                               / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp                  )           &
     2303                             ) *   weight_substep(intermediate_timestep_count)
     2304!
     2305!--    Statistical Evaluation of w'u'.
     2306       sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) +                                               &
     2307                              ( flux_t(k)                                                          &
     2308                                * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )          &
     2309                                / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )          &
     2310                              + diss_t(k)                                                          &
     2311                                *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)              )          &
     2312                                / ( ABS( w_comp(k) ) + 1.0E-20_wp                       )          &
     2313                              ) *   weight_substep(intermediate_timestep_count)
     2314    ENDDO
     2315
     2316    DO  k = nzb_max_l+1, nzt
     2317
     2318       flux_d    = flux_t(k-1)
     2319       diss_d    = diss_t(k-1)
     2320!
     2321!--    Calculate the divergence of the velocity field. A respective correction is needed to overcome
     2322!--    numerical instabilities introduced by an insufficient reduction of divergences near
     2323!--    topography.
     2324       div = ( ( u_comp(k)       - ( u(k,j,i)   + u(k,j,i-1) ) ) * ddx                             &
     2325             + ( v_comp(k) + gv  - ( v(k,j,i)   + v(k,j,i-1) ) ) * ddy                             &
     2326             + ( w_comp(k)   * rho_air_zw(k)                                                       &
     2327             -   w_comp(k-1) * rho_air_zw(k-1)                                                     &
     2328               )                                   * drho_air(k) * ddzw(k)                         &
     2329             ) * 0.5_wp
     2330
     2331       tend(k,j,i) = tend(k,j,i) - (                                                               &
     2332                         ( flux_r(k) + diss_r(k)                                                   &
     2333                       -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx                             &
     2334                       + ( flux_n(k) + diss_n(k)                                                   &
     2335                       -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy                             &
     2336                       + ( ( flux_t(k) + diss_t(k) )                                               &
     2337                       -   ( flux_d    + diss_d    )                                               &
     2338                                                 ) * drho_air(k) * ddzw(k)                         &
     2339                                   ) + div * u(k,j,i)
     2340
     2341       flux_l_u(k,j,tn) = flux_r(k)
     2342       diss_l_u(k,j,tn) = diss_r(k)
     2343       flux_s_u(k,tn)   = flux_n(k)
     2344       diss_s_u(k,tn)   = diss_n(k)
     2345!
     2346!--    Statistical Evaluation of u'u'. The factor has to be applied for right evaluation when
     2347!--    gallilei_trans = .T. .
     2348       sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) +                                                 &
     2349                             ( flux_r(k)                                                           &
     2350                               * ( u_comp(k) - 2.0_wp * hom(k,1,1,0)                   )           &
     2351                               / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )           &
     2352                             + diss_r(k)                                                           &
     2353                               *   ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0)              )           &
     2354                               / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp                  )           &
     2355                             ) *   weight_substep(intermediate_timestep_count)
     2356!
     2357!--    Statistical Evaluation of w'u'.
     2358       sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) +                                               &
     2359                              ( flux_t(k)                                                          &
     2360                                * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )          &
     2361                                / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )          &
     2362                              + diss_t(k)                                                          &
     2363                                *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)              )          &
     2364                                / ( ABS( w_comp(k) ) + 1.0E-20_wp                       )          &
     2365                              ) *   weight_substep(intermediate_timestep_count)
     2366    ENDDO
     2367
     2368
     2369
     2370 END SUBROUTINE advec_u_ws_ij
     2371
     2372
     2373
     2374!--------------------------------------------------------------------------------------------------!
     2375! Description:
     2376! ------------
     2377!> Advection of v-component - Call for grid point i,j
     2378!--------------------------------------------------------------------------------------------------!
     2379 SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn )
     2380
     2381
     2382    INTEGER(iwp)  ::  i         !< grid index along x-direction
     2383    INTEGER(iwp)  ::  i_omp     !< leftmost index on subdomain, or in case of OpenMP, on thread
     2384    INTEGER(iwp)  ::  j         !< grid index along y-direction
     2385    INTEGER(iwp)  ::  k         !< grid index along z-direction
     2386    INTEGER(iwp)  ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
     2387    INTEGER(iwp)  ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
     2388    INTEGER(iwp)  ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
     2389    INTEGER(iwp)  ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
     2390    INTEGER(iwp)  ::  tn        !< number of OpenMP thread
     2391
     2392    REAL(wp)      ::  diss_d   !< artificial dissipation term at grid box bottom
     2393    REAL(wp)      ::  div      !< divergence on v-grid
     2394    REAL(wp)      ::  flux_d   !< 6th-order flux at grid box bottom
     2395    REAL(wp)      ::  gu       !< Galilei-transformation velocity along x
     2396    REAL(wp)      ::  gv       !< Galilei-transformation velocity along y
     2397    REAL(wp)      ::  ibit9    !< flag indicating 1st-order scheme along x-direction
     2398    REAL(wp)      ::  ibit10   !< flag indicating 3rd-order scheme along x-direction
     2399    REAL(wp)      ::  ibit11   !< flag indicating 5th-order scheme along x-direction
     2400    REAL(wp)      ::  ibit12   !< flag indicating 1st-order scheme along y-direction
     2401    REAL(wp)      ::  ibit13   !< flag indicating 3rd-order scheme along y-direction
     2402    REAL(wp)      ::  ibit14   !< flag indicating 3rd-order scheme along y-direction
     2403    REAL(wp)      ::  ibit15   !< flag indicating 1st-order scheme along z-direction
     2404    REAL(wp)      ::  ibit16   !< flag indicating 3rd-order scheme along z-direction
     2405    REAL(wp)      ::  ibit17   !< flag indicating 3rd-order scheme along z-direction
     2406    REAL(wp)      ::  v_comp_l !< advection velocity along y on leftmost grid point on subdomain
     2407
     2408    REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
     2409    REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
     2410    REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_t !< discretized artificial dissipation at top of the grid box
     2411    REAL(wp), DIMENSION(nzb:nzt+1)  ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
     2412    REAL(wp), DIMENSION(nzb:nzt+1)  ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
     2413    REAL(wp), DIMENSION(nzb:nzt+1)  ::  flux_t !< discretized 6th-order flux at top of the grid box
     2414    REAL(wp), DIMENSION(nzb:nzt+1)  ::  u_comp !< advection velocity along x
     2415    REAL(wp), DIMENSION(nzb:nzt+1)  ::  v_comp !< advection velocity along y
     2416    REAL(wp), DIMENSION(nzb:nzt+1)  ::  w_comp !< advection velocity along z
     2417!
     2418!-- Used local modified copy of nzb_max (used to degrade order of discretization) at non-cyclic
     2419!-- boundaries. Modify only at relevant points instead of the entire subdomain. This should lead to
     2420!-- better load balance between boundary and non-boundary PEs.
     2421    IF( ( bc_dirichlet_l  .OR.  bc_radiation_l )  .AND.  i <= nxl + 2  .OR.                        &
     2422        ( bc_dirichlet_r  .OR.  bc_radiation_r )  .AND.  i >= nxr - 2  .OR.                        &
     2423        ( bc_dirichlet_s  .OR.  bc_radiation_s )  .AND.  j <= nys + 2  .OR.                        &
     2424        ( bc_dirichlet_n  .OR.  bc_radiation_n )  .AND.  j >= nyn - 2 )  THEN
     2425       nzb_max_l = nzt
     2426    ELSE
     2427       nzb_max_l = nzb_max
     2428    END IF
     2429
     2430    gu = 2.0_wp * u_gtrans
     2431    gv = 2.0_wp * v_gtrans
     2432
     2433!
     2434!-- Compute leftside fluxes for the respective boundary.
     2435    IF ( i == i_omp )  THEN
     2436
     2437       DO  k = nzb+1, nzb_max_l
     2438
     2439          ibit11 = REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp )
     2440          ibit10 = REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp )
     2441          ibit9  = REAL( IBITS(advc_flags_m(k,j,i-1),9,1),  KIND = wp )
     2442
     2443          u_comp(k)        = u(k,j-1,i) + u(k,j,i) - gu
     2444          flux_l_v(k,j,tn) = u_comp(k) * (                                                         &
     2445                             (  37.0_wp * ibit11 * adv_mom_5                                       &
     2446                             +   7.0_wp * ibit10 * adv_mom_3                                       &
     2447                             +            ibit9  * adv_mom_1 ) * ( v(k,j,i)   + v(k,j,i-1) )       &
     2448                             - ( 8.0_wp * ibit11 * adv_mom_5                                       &
     2449                             +            ibit10 * adv_mom_3 ) * ( v(k,j,i+1) + v(k,j,i-2) )       &
     2450                             + (          ibit11 * adv_mom_5 ) * ( v(k,j,i+2) + v(k,j,i-3) )       &
     2451                                         )
     2452
     2453          diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * (                                                &
     2454                             (  10.0_wp * ibit11 * adv_mom_5                                       &
     2455                             +   3.0_wp * ibit10 * adv_mom_3                                       &
     2456                             +            ibit9  * adv_mom_1 ) * ( v(k,j,i)   - v(k,j,i-1) )       &
     2457                             - ( 5.0_wp * ibit11 * adv_mom_5                                       &
     2458                             +            ibit10 * adv_mom_3 ) * ( v(k,j,i+1) - v(k,j,i-2) )       &
     2459                             + (          ibit11 * adv_mom_5 ) * ( v(k,j,i+2) - v(k,j,i-3) )       &
     2460                                                  )
     2461
     2462       ENDDO
     2463
     2464       DO  k = nzb_max_l+1, nzt
     2465
     2466          u_comp(k)        = u(k,j-1,i) + u(k,j,i) - gu
     2467          flux_l_v(k,j,tn) = u_comp(k) * (                                                         &
     2468                               37.0_wp * ( v(k,j,i)   + v(k,j,i-1) )                               &
     2469                             -  8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) )                               &
     2470                             +           ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
     2471          diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * (                                                &
     2472                               10.0_wp * ( v(k,j,i)   - v(k,j,i-1) )                               &
     2473                             -  5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) )                               &
     2474                             +           ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
     2475
     2476       ENDDO
     2477
     2478    ENDIF
     2479!
     2480!-- Compute southside fluxes for the respective boundary.
     2481    IF ( j == nysv )  THEN
     2482
     2483       DO  k = nzb+1, nzb_max_l
     2484
     2485          ibit14 = REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp )
     2486          ibit13 = REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp )
     2487          ibit12 = REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp )
     2488
     2489          v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
     2490          flux_s_v(k,tn) = v_comp_l * (                                                            &
     2491                         (  37.0_wp * ibit14 * adv_mom_5                                           &
     2492                         +   7.0_wp * ibit13 * adv_mom_3                                           &
     2493                         +            ibit12 * adv_mom_1 ) * ( v(k,j,i)   + v(k,j-1,i) )           &
     2494                         - ( 8.0_wp * ibit14 * adv_mom_5                                           &
     2495                         +            ibit13 * adv_mom_3 ) * ( v(k,j+1,i) + v(k,j-2,i) )           &
     2496                         + (          ibit14 * adv_mom_5 ) * ( v(k,j+2,i) + v(k,j-3,i) )           &
     2497                                      )
     2498
     2499          diss_s_v(k,tn) = - ABS( v_comp_l ) * (                                                   &
     2500                           (  10.0_wp * ibit14 * adv_mom_5                                         &
     2501                           +   3.0_wp * ibit13 * adv_mom_3                                         &
     2502                           +            ibit12 * adv_mom_1 ) * ( v(k,j,i)   - v(k,j-1,i) )         &
     2503                           - ( 5.0_wp * ibit14 * adv_mom_5                                         &
     2504                           +            ibit13 * adv_mom_3 ) * ( v(k,j+1,i) - v(k,j-2,i) )         &
     2505                           + (          ibit14 * adv_mom_5 ) * ( v(k,j+2,i) - v(k,j-3,i) )         &
     2506                                               )
     2507
     2508       ENDDO
     2509
     2510       DO  k = nzb_max_l+1, nzt
     2511
     2512          v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
     2513          flux_s_v(k,tn) = v_comp_l * (                                     &
     2514                        37.0_wp * ( v(k,j,i)   + v(k,j-1,i)   )             &
     2515                      -  8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) )               &
     2516                      +           ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
     2517          diss_s_v(k,tn) = - ABS( v_comp_l ) * (                            &
     2518                        10.0_wp * ( v(k,j,i)   - v(k,j-1,i)   )             &
     2519                      -  5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) )               &
     2520                      +           ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
     2521
     2522       ENDDO
     2523
     2524    ENDIF
     2525!
     2526!--    Now compute the fluxes and tendency terms for the horizontal and
     2527!--    verical parts.
     2528    DO  k = nzb+1, nzb_max_l
     2529
     2530       ibit11 = REAL( IBITS(advc_flags_m(k,j,i),11,1), KIND = wp )
     2531       ibit10 = REAL( IBITS(advc_flags_m(k,j,i),10,1), KIND = wp )
     2532       ibit9  = REAL( IBITS(advc_flags_m(k,j,i),9,1),  KIND = wp )
     2533
     2534       u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
     2535       flux_r(k) = u_comp(k) * (                                                                   &
     2536                  (  37.0_wp * ibit11 * adv_mom_5                                                  &
     2537                  +   7.0_wp * ibit10 * adv_mom_3                                                  &
     2538                  +            ibit9  * adv_mom_1 ) * ( v(k,j,i+1) + v(k,j,i)   )                  &
     2539                  - ( 8.0_wp * ibit11 * adv_mom_5                                                  &
     2540                  +            ibit10 * adv_mom_3 ) * ( v(k,j,i+2) + v(k,j,i-1) )                  &
     2541                  + (          ibit11 * adv_mom_5 ) * ( v(k,j,i+3) + v(k,j,i-2) )                  &
     2542                               )