Changeset 4466


Ignore:
Timestamp:
Mar 20, 2020 4:14:41 PM (13 months ago)
Author:
suehring
Message:

Vector branch in advec_ws optimized, symmetric boundary conditions implemented in vector branch

Location:
palm/trunk/SOURCE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r4458 r4466  
    2525# -----------------
    2626# $Id$
     27# cpu measures in advec_ws added
     28#
     29# 4458 2020-03-11 15:37:31Z raasch
    2730# bugfix for r4457: missing dependency added
    28 # 
     31#
    2932# 4457 2020-03-11 14:20:43Z raasch
    3033# exchange horiz has been modularized and exchange horiz 2d has been merged, dependencies updated
    3134# accordingly
    32 # 
     35#
    3336# 4453 2020-03-11 08:10:13Z raasch
    3437# dependencies for exchange horiz modified
    35 # 
     38#
    3639# 4434 2020-03-03 10:02:18Z oliver.maas
    3740# Added output control for wind turbine model
    38 # 
     41#
    3942# 4414 2020-02-19 20:16:04Z suehring
    4043# Move dependencies for init grid from advection scheme and multigrid solver
    4144# to module_interface
    42 # 
     45#
    4346# 4411 2020-02-18 14:28:02Z maronga
    4447# Added output routines for WTM
    45 # 
     48#
    4649# 4400 2020-02-10 20:32:41Z suehring
    4750# Add dependency for data-output module
    48 # 
     51#
    4952# 4392 2020-01-31 16:14:57Z pavelkrc
    5053# add dependency on fft for transpose
    51 # 
     54#
    5255# 4347 2019-12-18 13:18:33Z suehring
    5356# add dependency to basic_constants_and_equations_mod for dynamics_mod
    54 # 
     57#
    5558# 4331 2019-12-10 18:25:02Z suehring
    5659# Move diagnostic surface output to diagnostic_output_quantities
    57 # 
     60#
    5861# 4309 2019-11-26 18:49:59Z suehring
    5962# Add dependency to parallel random generator for synthetic turbulence generator
    60 # 
     63#
    6164# 4286 2019-10-30 16:01:14Z resler
    6265# delete boundary_conds, added missing dependencies
    63 # 
     66#
    6467# 4270 2019-10-23 10:46:20Z monakurppa
    6568# - Implement offline nesting for salsa and add dependency of nesting_offl_mod
     
    6871#   for salsa
    6972# - Add dependency on basic_constants_and_equations_mod for salsa_mod
    70 # 
     73#
    7174# 4258 2019-10-07 13:29:08Z suehring
    7275# Add dependency of land-surface model on pmc_handle_communicator and cpu_log
    73 # 
     76#
    7477# 4245 2019-09-30 08:40:37Z pavelkrc
    7578# Remove no longer needed dependencies on surface_mod
    76 # 
     79#
    7780# 4227 2019-09-10 18:04:34Z gronemeier
    7881# Add palm_date_time_mod, remove date_and_time_mod
     
    8083# 4223 2019-09-10 09:20:47Z gronemeier
    8184# Corrected "Former revisions" section
    82 # 
     85#
    8386# 4174 2019-08-20 12:41:13Z gronemeier
    8487# bugfix: add missing dependencies for vdi_internal_controls
    85 # 
     88#
    8689# 4173 2019-08-20 12:04:06Z gronemeier
    8790# add vdi_internal_controls
    88 # 
     91#
    8992# 4168 2019-08-16 13:50:17Z suehring
    9093# Remove some dependencies on surface_mod that are no longer required without
    9194# get_topography_top_index functions
    92 # 
     95#
    9396# 4167 2019-08-16 11:01:48Z suehring
    9497# Remove no longer needed dependencies on surface_mod
    95 # 
     98#
    9699# 4127 2019-07-30 14:47:10Z suehring
    97 # Add dependency of data_output_3d on plant_canopy_model_mod 
     100# Add dependency of data_output_3d on plant_canopy_model_mod
    98101# (merge from branch resler)
    99102#
    100103# 4106 2019-07-19 08:54:42Z gronemeier
    101104# Remove dependency on pmc_interface for boundary_conds
    102 # 
     105#
    103106# 4070 2019-07-03 13:51:40Z gronemeier
    104107# Add new data output modules
     
    358361        modules.o
    359362advec_ws.o: \
     363        cpulog_mod.o \
    360364        exchange_horiz_mod.o \
    361365        mod_kinds.o \
     
    825829        modules.o \
    826830        netcdf_data_input_mod.o \
    827         salsa_mod.o 
     831        salsa_mod.o
    828832netcdf_data_input_mod.o: \
    829833        chem_modules.o \
  • palm/trunk/SOURCE/advec_ws.f90

    r4457 r4466  
    2020! Current revisions:
    2121! ------------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! - vector branch further optimized (linear dependencies along z removed and
     28!   loops are splitted)
     29! - topography closed channel flow with symmetric boundaries also implemented
     30!   in vector branch
     31! - some formatting adjustments made and comments added
     32! - cpu measures for vector branch added
     33!
     34! 4457 2020-03-11 14:20:43Z raasch
    2735! use statement for exchange horiz added
    28 ! 
     36!
    2937! 4414 2020-02-19 20:16:04Z suehring
    3038! Move call for initialization of control flags to ws_init
    31 ! 
     39!
    3240! 4360 2020-01-07 11:25:50Z suehring
    3341! Introduction of wall_flags_total_0, which currently sets bits based on static
    3442! topography information used in wall_flags_static_0
    35 ! 
     43!
    3644! 4340 2019-12-16 08:17:03Z Giersch
    37 ! Topography closed channel flow with symmetric boundaries implemented 
    38 ! 
     45! Topography closed channel flow with symmetric boundaries implemented
     46!
    3947! 4330 2019-12-10 16:16:33Z knoop
    4048! Bugix: removed syntax error introduced by last commit
    41 ! 
     49!
    4250! 4329 2019-12-10 15:46:36Z motisi
    4351! Renamed wall_flags_0 to wall_flags_static_0
    44 ! 
     52!
    4553! 4328 2019-12-09 18:53:04Z suehring
    4654! Minor formatting adjustments
    47 ! 
     55!
    4856! 4327 2019-12-06 14:48:31Z Giersch
    4957! Setting of advection flags for vertical fluxes of w revised, air density for
    5058! vertical flux calculation of w at k=1 is considered now
    51 ! 
     59!
    5260! 4325 2019-12-06 07:14:04Z Giersch
    53 ! Vertical fluxes of w are now set to zero at nzt and nzt+1, setting of 
     61! Vertical fluxes of w are now set to zero at nzt and nzt+1, setting of
    5462! advection flags for fluxes in z-direction revised, comments extended
    55 ! 
     63!
    5664! 4324 2019-12-06 07:11:33Z Giersch
    5765! Indirect indexing for calculating vertical fluxes close to boundaries is only
    5866! used for loop indizes where it is really necessary
    59 ! 
     67!
    6068! 4317 2019-12-03 12:43:22Z Giersch
    61 ! Comments revised/added, formatting improved, fluxes for u,v, and scalars are 
     69! Comments revised/added, formatting improved, fluxes for u,v, and scalars are
    6270! explicitly set to zero at nzt+1, fluxes of w-component are now calculated only
    6371! until nzt-1 (Prognostic equation for w-velocity component ends at nzt-1)
    64 ! 
     72!
    6573! 4204 2019-08-30 12:30:17Z knoop
    6674! Bugfix: Changed sk_num initialization default to avoid implicit SAVE-Attribut
    67 ! 
     75!
    6876! 4182 2019-08-22 15:20:23Z scharf
    6977! Corrected "Former revisions" section
    70 ! 
     78!
    7179! 4110 2019-07-22 17:05:21Z suehring
    7280! - Separate initialization of advection flags for momentum and scalars. In this
    73 !   context, resort the bits and do some minor formatting. 
    74 ! - Make flag initialization for scalars more flexible, introduce an 
    75 !   arguemnt list to indicate non-cyclic boundaries (required for decycled 
     81!   context, resort the bits and do some minor formatting.
     82! - Make flag initialization for scalars more flexible, introduce an
     83!   arguemnt list to indicate non-cyclic boundaries (required for decycled
    7684!   scalars such as chemical species or aerosols)
    77 ! - Introduce extended 'degradation zones', where horizontal advection of 
    78 !   passive scalars is discretized by first-order scheme at all grid points 
    79 !   that in the vicinity of buildings (<= 3 grid points). Even though no 
    80 !   building is within the numerical stencil, first-order scheme is used. 
     85! - Introduce extended 'degradation zones', where horizontal advection of
     86!   passive scalars is discretized by first-order scheme at all grid points
     87!   that in the vicinity of buildings (<= 3 grid points). Even though no
     88!   building is within the numerical stencil, first-order scheme is used.
    8189!   At fourth and fifth grid point the order of the horizontal advection scheme
    8290!   is successively upgraded.
    83 !   These extended degradation zones are used to avoid stationary numerical 
     91!   These extended degradation zones are used to avoid stationary numerical
    8492!   oscillations, which are responsible for high concentration maxima that may
    85 !   appear under shear-free stable conditions. 
    86 ! - Change interface for scalar advection routine. 
     93!   appear under shear-free stable conditions.
     94! - Change interface for scalar advection routine.
    8795! - Bugfix, avoid uninitialized value sk_num in vector version of scalar
    8896!   advection
    89 ! 
     97!
    9098! 4109 2019-07-22 17:00:34Z suehring
    91 ! Implementation of a flux limiter according to Skamarock (2006) for the 
    92 ! vertical scalar advection. Please note, this is only implemented for the 
     99! Implementation of a flux limiter according to Skamarock (2006) for the
     100! vertical scalar advection. Please note, this is only implemented for the
    93101! cache-optimized version at the moment. Implementation for the vector-
    94 ! optimized version will follow after critical issues concerning 
     102! optimized version will follow after critical issues concerning
    95103! vectorization are fixed.
    96 ! 
     104!
    97105! 3873 2019-04-08 15:44:30Z knoop
    98106! Moved ocean_mode specific code to ocean_mod
    99 ! 
     107!
    100108! 3872 2019-04-08 15:03:06Z knoop
    101109! Moved all USE statements to module level + removed salsa dependency
    102 ! 
     110!
    103111! 3871 2019-04-08 14:38:39Z knoop
    104112! Moving initialization of bcm specific flux arrays into bulk_cloud_model_mod
    105 ! 
     113!
    106114! 3864 2019-04-05 09:01:56Z monakurppa
    107115! Remove tailing white spaces
    108 ! 
     116!
    109117! 3696 2019-01-24 16:37:35Z suehring
    110118! Bugfix in degradation height
    111 ! 
     119!
    112120! 3661 2019-01-08 18:22:50Z suehring
    113 ! - Minor bugfix in divergence correction (only has implications at 
     121! - Minor bugfix in divergence correction (only has implications at
    114122!   downward-facing wall surfaces)
    115123! - Remove setting of Neumann condition for horizontal velocity variances
    116124! - Split loops for tendency calculation and divergence correction in order to
    117125!   reduce bit queries
    118 ! - Introduce new parameter nzb_max_l to better control order degradation at 
     126! - Introduce new parameter nzb_max_l to better control order degradation at
    119127!   non-cyclic boundaries
    120 ! 
     128!
    121129! 3655 2019-01-07 16:51:22Z knoop
    122130! OpenACC port for SPEC
     
    133141! ------------
    134142!> Advection scheme for scalars and momentum using the flux formulation of
    135 !> Wicker and Skamarock 5th order. Additionally the module contains of a 
    136 !> routine using for initialisation and steering of the statical evaluation. 
     143!> Wicker and Skamarock 5th order. Additionally the module contains of a
     144!> routine using for initialisation and steering of the statical evaluation.
    137145!> The computation of turbulent fluxes takes place inside the advection
    138146!> routines.
     
    140148!> degraded.
    141149!> A divergence correction is applied. It is necessary for topography, since
    142 !> the divergence is not sufficiently reduced, resulting in erroneous fluxes 
    143 !> and could lead to numerical instabilities. 
     150!> the divergence is not sufficiently reduced, resulting in erroneous fluxes
     151!> and could lead to numerical instabilities.
    144152!>
    145153!> @todo Implement monotonic flux limiter also for vector version.
     
    154162               u_stokes_zu, v_stokes_zu,                                       &
    155163               diss_l_diss, diss_l_e, diss_l_pt, diss_l_q,                     &
    156                diss_l_s, diss_l_sa, diss_l_u, diss_l_v, diss_l_w,              &
     164               diss_l_s, diss_l_u, diss_l_v, diss_l_w,                         &
    157165               flux_l_diss, flux_l_e, flux_l_pt, flux_l_q, flux_l_s,           &
    158                flux_l_sa, flux_l_u, flux_l_v, flux_l_w,                        &
     166               flux_l_u, flux_l_v, flux_l_w,                                   &
    159167               diss_s_diss, diss_s_e, diss_s_pt, diss_s_q, diss_s_s,           &
    160                diss_s_sa, diss_s_u, diss_s_v, diss_s_w,                        &
     168               diss_s_u, diss_s_v, diss_s_w,                                   &
    161169               flux_s_diss, flux_s_e, flux_s_pt, flux_s_q, flux_s_s,           &
    162                flux_s_sa, flux_s_u, flux_s_v, flux_s_w
     170               flux_s_u, flux_s_v, flux_s_w
    163171
    164172    USE control_parameters,                                                    &
    165         ONLY:  air_chemistry,                                                  &
    166                bc_dirichlet_l,                                                 &
     173        ONLY:  bc_dirichlet_l,                                                 &
    167174               bc_dirichlet_n,                                                 &
    168175               bc_dirichlet_r,                                                 &
     
    176183               passive_scalar,                                                 &
    177184               rans_tke_e,                                                     &
    178                momentum_advec,                                                 &
    179                salsa,                                                          &
    180                scalar_advec,                                                   &
    181185               symmetry_flag,                                                  &
    182186               intermediate_timestep_count,                                    &
     
    186190               ws_scheme_sca,                                                  &
    187191               dt_3d
     192
     193    USE cpulog,                                                                &
     194        ONLY:  cpu_log,                                                        &
     195               log_point_s
    188196
    189197    USE exchange_horiz_mod,                                                    &
     
    198206               nxlg,                                                           &
    199207               nxlu,                                                           &
    200                nxr,                                                            & 
    201                nxrg,                                                           & 
     208               nxr,                                                            &
     209               nxrg,                                                           &
    202210               ny,                                                             &
    203                nyn,                                                            & 
    204                nyng,                                                           & 
     211               nyn,                                                            &
     212               nyng,                                                           &
    205213               nys,                                                            &
    206214               nysg,                                                           &
     
    242250    INTERFACE ws_init
    243251       MODULE PROCEDURE ws_init
    244     END INTERFACE ws_init         
    245              
     252    END INTERFACE ws_init
     253
    246254    INTERFACE ws_init_flags_momentum
    247255       MODULE PROCEDURE ws_init_flags_momentum
    248256    END INTERFACE ws_init_flags_momentum
    249    
     257
    250258    INTERFACE ws_init_flags_scalar
    251259       MODULE PROCEDURE ws_init_flags_scalar
     
    287295
    288296!
    289 !--    Set the appropriate factors for scalar and momentum advection.
     297!--    Set factors for scalar and momentum advection.
    290298       adv_sca_5 = 1.0_wp /  60.0_wp
    291299       adv_sca_3 = 1.0_wp /  12.0_wp
     
    294302       adv_mom_3 = 1.0_wp /  24.0_wp
    295303       adv_mom_1 = 1.0_wp /   4.0_wp
    296 !         
     304!
    297305!--    Arrays needed for statical evaluation of fluxes.
    298306       IF ( ws_scheme_mom )  THEN
     
    330338
    331339!
    332 !--    Arrays needed for reasons of speed optimization for cache version.
    333 !--    For the vector version the buffer arrays are not necessary,
    334 !--    because the the fluxes can swapped directly inside the loops of the
    335 !--    advection routines.
     340!--    Arrays needed for reasons of speed optimization
     341       IF ( ws_scheme_mom )  THEN
     342
     343          ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1),               &
     344                    flux_s_v(nzb+1:nzt,0:threads_per_task-1),               &
     345                    flux_s_w(nzb+1:nzt,0:threads_per_task-1),               &
     346                    diss_s_u(nzb+1:nzt,0:threads_per_task-1),               &
     347                    diss_s_v(nzb+1:nzt,0:threads_per_task-1),               &
     348                    diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
     349          ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
     350                    flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
     351                    flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
     352                    diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
     353                    diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
     354                    diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
     355
     356       ENDIF
     357!
     358!--    For the vector version the buffer arrays for scalars are not necessary,
     359!--    since internal arrays are used in the vector version.
    336360       IF ( loop_optimization /= 'vector' )  THEN
    337 
    338           IF ( ws_scheme_mom )  THEN
    339 
    340              ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1),               &
    341                        flux_s_v(nzb+1:nzt,0:threads_per_task-1),               &
    342                        flux_s_w(nzb+1:nzt,0:threads_per_task-1),               &
    343                        diss_s_u(nzb+1:nzt,0:threads_per_task-1),               &
    344                        diss_s_v(nzb+1:nzt,0:threads_per_task-1),               &
    345                        diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
    346              ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
    347                        flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
    348                        flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
    349                        diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
    350                        diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
    351                        diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    352 
    353           ENDIF
    354 
    355361          IF ( ws_scheme_sca )  THEN
    356362
     
    358364                       flux_s_e(nzb+1:nzt,0:threads_per_task-1),               &
    359365                       diss_s_pt(nzb+1:nzt,0:threads_per_task-1),              &
    360                        diss_s_e(nzb+1:nzt,0:threads_per_task-1) ) 
     366                       diss_s_e(nzb+1:nzt,0:threads_per_task-1) )
    361367             ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),      &
    362368                       flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1),       &
     
    386392
    387393          ENDIF
    388 
    389394       ENDIF
    390395!
    391 !--    Initialize the flag arrays controlling degradation near walls, i.e. 
     396!--    Initialize the flag arrays controlling degradation near walls, i.e.
    392397!--    to decrease the numerical stencil appropriately. The order of the scheme
    393398!--    is degraded near solid walls as well as near non-cyclic inflow and outflow
     
    412417!> Initialization of flags to control the order of the advection scheme near
    413418!> solid walls and non-cyclic inflow boundaries, where the order is sucessively
    414 !> degraded. 
     419!> degraded.
    415420!------------------------------------------------------------------------------!
    416421    SUBROUTINE ws_init_flags_momentum
     
    437442          DO  j = nys, nyn
    438443             DO  k = nzb+1, nzt
    439 ! 
    440 !--             At first, set flags to WS1. 
    441 !--             Since fluxes are swapped in advec_ws.f90, this is necessary to 
     444!
     445!--             At first, set flags to WS1.
     446!--             Since fluxes are swapped in advec_ws.f90, this is necessary to
    442447!--             in order to handle the left/south flux.
    443 !--             near vertical walls. 
     448!--             near vertical walls.
    444449                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 )
    445450                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 )
     
    452457                         ( ( bc_dirichlet_r .OR. bc_radiation_r )              &
    453458                           .AND. i == nxr   ) )                                &
    454                 THEN                                                           
    455                     advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 )     
     459                THEN
     460                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 )
    456461                ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),1)  .AND.   &
    457462                                 BTEST(wall_flags_total_0(k,j,i+1),1)  .OR.    &
     
    462467                         ( ( bc_dirichlet_l .OR. bc_radiation_l )              &
    463468                           .AND. i == nxlu+1) )                                &
    464                 THEN                                                           
    465                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 1 )       
    466 !                                                                             
    467 !--                Clear flag for WS1                                         
    468                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 )       
    469                 ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),1)  .AND.          &
    470                          BTEST(wall_flags_total_0(k,j,i+2),1)  .AND.          &
    471                          BTEST(wall_flags_total_0(k,j,i-1),1) )               &
    472                 THEN                                                           
    473                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 2 )       
    474 !                                                                             
    475 !--                Clear flag for WS1                                         
    476                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 )       
    477                 ENDIF                                                         
    478 !                                                                             
    479 !--             u component - y-direction                                     
    480 !--             WS1 (3), WS3 (4), WS5 (5)                                       
     469                THEN
     470                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 1 )
     471!
     472!--                Clear flag for WS1
     473                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 )
     474                ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),1)  .AND.           &
     475                         BTEST(wall_flags_total_0(k,j,i+2),1)  .AND.           &
     476                         BTEST(wall_flags_total_0(k,j,i-1),1) )                &
     477                THEN
     478                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 2 )
     479!
     480!--                Clear flag for WS1
     481                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 )
     482                ENDIF
     483!
     484!--             u component - y-direction
     485!--             WS1 (3), WS3 (4), WS5 (5)
    481486                IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),1)   .OR.         &
    482487                         ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
     
    484489                         ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
    485490                           .AND. j == nyn   ) )                                &
    486                 THEN                                                           
    487                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 )       
     491                THEN
     492                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 )
    488493                ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),1)  .AND.   &
    489494                                 BTEST(wall_flags_total_0(k,j+1,i),1)  .OR.    &
     
    494499                         ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
    495500                           .AND. j == nyn-1 ) )                                &
    496                 THEN                                                           
    497                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 4 )       
    498 !                                                                             
    499 !--                Clear flag for WS1                                         
    500                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 )       
     501                THEN
     502                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 4 )
     503!
     504!--                Clear flag for WS1
     505                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 )
    501506                ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),1)  .AND.          &
    502507                         BTEST(wall_flags_total_0(k,j+2,i),1)  .AND.          &
    503508                         BTEST(wall_flags_total_0(k,j-1,i),1) )               &
    504                 THEN                                                           
    505                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 5 )       
    506 !                                                                             
    507 !--                Clear flag for WS1                                         
    508                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 )       
    509                 ENDIF                                                         
    510 !                                                                             
    511 !--             u component - z-direction. Fluxes are calculated on w-grid                                     
    512 !--             level. Boundary u-values at/within walls aren't used. 
    513 !--             WS1 (6), WS3 (7), WS5 (8)                                     
    514                 IF ( k == nzb+1 )  THEN                                       
    515                    k_mm = nzb                                                 
    516                 ELSE                                                           
    517                    k_mm = k - 2                                               
    518                 ENDIF                                                         
    519                 IF ( k > nzt-1 )  THEN                                         
    520                    k_pp = nzt+1                                               
    521                 ELSE                                                           
    522                    k_pp = k + 2                                               
    523                 ENDIF                                                         
    524                 IF ( k > nzt-2 )  THEN                                         
    525                    k_ppp = nzt+1                                               
    526                 ELSE                                                           
    527                    k_ppp = k + 3                                               
    528                 ENDIF                                                         
    529                                                                                
    530                 flag_set = .FALSE.                                             
    531                 IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),1)       .AND.     &
    532                              BTEST(wall_flags_total_0(k,j,i),1)         .AND.     &
    533                              BTEST(wall_flags_total_0(k+1,j,i),1) )     .OR.      &
    534                      ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),1)      .AND.     &                             
    535                              BTEST(wall_flags_total_0(k+1,j,i),1)       .AND.     &
    536                              BTEST(wall_flags_total_0(k,j,i),1) )       .OR.      &
    537                      ( k == nzt .AND. symmetry_flag == 0 ) )                      &
    538                 THEN                                                           
    539                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 6 )       
    540                    flag_set = .TRUE.                                           
    541                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),1)    .OR.    &
    542                            .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),1) ) .AND.   &
    543                                  BTEST(wall_flags_total_0(k-1,j,i),1)     .AND.   &
    544                                  BTEST(wall_flags_total_0(k,j,i),1)       .AND.   &
    545                                  BTEST(wall_flags_total_0(k+1,j,i),1)     .AND.   &
    546                                  BTEST(wall_flags_total_0(k_pp,j,i),1)    .AND.   &
    547                            .NOT. flag_set                                  .OR.   &
    548                          ( k == nzt - 1 .AND. symmetry_flag == 0 ) )              &
    549                 THEN                                                           
    550                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 )       
    551                    flag_set = .TRUE.                                           
    552                 ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),1)            .AND.   &
    553                          BTEST(wall_flags_total_0(k-1,j,i),1)             .AND.   &
    554                          BTEST(wall_flags_total_0(k,j,i),1)               .AND.   &
    555                          BTEST(wall_flags_total_0(k+1,j,i),1)             .AND.   &
    556                          BTEST(wall_flags_total_0(k_pp,j,i),1)            .AND.   &
    557                          BTEST(wall_flags_total_0(k_ppp,j,i),1)           .AND.   &
    558                          .NOT. flag_set )                                         &
     509                THEN
     510                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 5 )
     511!
     512!--                Clear flag for WS1
     513                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 )
     514                ENDIF
     515!
     516!--             u component - z-direction. Fluxes are calculated on w-grid
     517!--             level. Boundary u-values at/within walls aren't used.
     518!--             WS1 (6), WS3 (7), WS5 (8)
     519                IF ( k == nzb+1 )  THEN
     520                   k_mm = nzb
     521                ELSE
     522                   k_mm = k - 2
     523                ENDIF
     524                IF ( k > nzt-1 )  THEN
     525                   k_pp = nzt+1
     526                ELSE
     527                   k_pp = k + 2
     528                ENDIF
     529                IF ( k > nzt-2 )  THEN
     530                   k_ppp = nzt+1
     531                ELSE
     532                   k_ppp = k + 3
     533                ENDIF
     534
     535                flag_set = .FALSE.
     536                IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),1)       .AND.  &
     537                             BTEST(wall_flags_total_0(k,j,i),1)         .AND.  &
     538                             BTEST(wall_flags_total_0(k+1,j,i),1) )     .OR.   &
     539                     ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),1)      .AND.  &
     540                             BTEST(wall_flags_total_0(k+1,j,i),1)       .AND.  &
     541                             BTEST(wall_flags_total_0(k,j,i),1) )       .OR.   &
     542                     ( k == nzt .AND. symmetry_flag == 0 ) )                   &
     543                THEN
     544                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 6 )
     545                   flag_set = .TRUE.
     546                ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),1)    .OR. &
     547                           .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),1) ) .AND.&
     548                                 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)     .AND.&
     551                                 BTEST(wall_flags_total_0(k_pp,j,i),1)    .AND.&
     552                           .NOT. flag_set                                  .OR.&
     553                         ( k == nzt - 1 .AND. symmetry_flag == 0 ) )           &
     554                THEN
     555                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 )
     556                   flag_set = .TRUE.
     557                ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),1)            .AND.&
     558                         BTEST(wall_flags_total_0(k-1,j,i),1)             .AND.&
     559                         BTEST(wall_flags_total_0(k,j,i),1)               .AND.&
     560                         BTEST(wall_flags_total_0(k+1,j,i),1)             .AND.&
     561                         BTEST(wall_flags_total_0(k_pp,j,i),1)            .AND.&
     562                         BTEST(wall_flags_total_0(k_ppp,j,i),1)           .AND.&
     563                         .NOT. flag_set )                                      &
    559564                THEN
    560565                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 8 )
     
    570575             DO  k = nzb+1, nzt
    571576!
    572 !--             At first, set flags to WS1. 
    573 !--             Since fluxes are swapped in advec_ws.f90, this is necessary to 
     577!--             At first, set flags to WS1.
     578!--             Since fluxes are swapped in advec_ws.f90, this is necessary to
    574579!--             in order to handle the left/south flux.
    575580                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9  )
     
    583588                         ( ( bc_dirichlet_r .OR. bc_radiation_r )              &
    584589                           .AND. i == nxr   ) )                                &
    585                THEN                                                           
    586                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 )       
    587 !                                                                             
    588 !--             WS3                                                           
     590               THEN
     591                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 )
     592!
     593!--             WS3
    589594                ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),2)   .AND.  &
    590595                                 BTEST(wall_flags_total_0(k,j,i+1),2) ) .OR.   &
     
    595600                         ( ( bc_dirichlet_l .OR. bc_radiation_l )              &
    596601                           .AND. i == nxlu  ) )                                &
    597                 THEN                                                           
    598                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 10 )     
    599 !                                                                             
    600 !--                Clear flag for WS1                                         
    601                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 )     
     602                THEN
     603                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 10 )
     604!
     605!--                Clear flag for WS1
     606                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 )
    602607                ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),2)  .AND.           &
    603608                         BTEST(wall_flags_total_0(k,j,i+2),2)  .AND.           &
    604609                         BTEST(wall_flags_total_0(k,j,i-1),2) )                &
    605                 THEN                                                           
    606                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 11 )     
    607 !                                                                             
    608 !--                Clear flag for WS1                                         
    609                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 )       
    610                 ENDIF                                                         
    611 !                                                                             
    612 !--             v component - y-direction                                     
    613 !--             WS1 (12), WS3 (13), WS5 (14)                                   
     610                THEN
     611                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 11 )
     612!
     613!--                Clear flag for WS1
     614                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 )
     615                ENDIF
     616!
     617!--             v component - y-direction
     618!--             WS1 (12), WS3 (13), WS5 (14)
    614619                IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),2) .OR.           &
    615620                         ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
     
    617622                         ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
    618623                           .AND. j == nyn   ) )                                &
    619                 THEN                                                           
    620                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 )     
     624                THEN
     625                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 )
    621626                ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),2)  .AND.   &
    622627                                 BTEST(wall_flags_total_0(k,j+1,i),2)  .OR.    &
     
    627632                         ( (  bc_dirichlet_n .OR. bc_radiation_n )             &
    628633                           .AND. j == nyn-1 ) )                                &
    629                 THEN                                                           
    630                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 13 )     
    631 !                                                                             
    632 !--                Clear flag for WS1                                         
    633                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 )     
     634                THEN
     635                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 13 )
     636!
     637!--                Clear flag for WS1
     638                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 )
    634639                ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),2)  .AND.           &
    635640                         BTEST(wall_flags_total_0(k,j+2,i),2)  .AND.           &
    636                          BTEST(wall_flags_total_0(k,j-1,i),2) )                & 
    637                 THEN                                                         
     641                         BTEST(wall_flags_total_0(k,j-1,i),2) )                &
     642                THEN
    638643                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 14 )
    639644!
     
    641646                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 )
    642647                ENDIF
    643 !                                                             
    644 !--             v component - z-direction. Fluxes are calculated on w-grid                                     
    645 !--             level. Boundary v-values at/within walls aren't used. 
     648!
     649!--             v component - z-direction. Fluxes are calculated on w-grid
     650!--             level. Boundary v-values at/within walls aren't used.
    646651!--             WS1 (15), WS3 (16), WS5 (17)
    647652                IF ( k == nzb+1 )  THEN
    648653                   k_mm = nzb
    649                 ELSE 
     654                ELSE
    650655                   k_mm = k - 2
    651656                ENDIF
    652657                IF ( k > nzt-1 )  THEN
    653658                   k_pp = nzt+1
    654                 ELSE 
     659                ELSE
    655660                   k_pp = k + 2
    656661                ENDIF
    657662                IF ( k > nzt-2 )  THEN
    658663                   k_ppp = nzt+1
    659                 ELSE 
     664                ELSE
    660665                   k_ppp = k + 3
    661                 ENDIF 
    662                
     666                ENDIF
     667
    663668                flag_set = .FALSE.
    664                 IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),2)          .AND.   &
    665                              BTEST(wall_flags_total_0(k,j,i),2)            .AND.   &
    666                              BTEST(wall_flags_total_0(k+1,j,i),2) )        .OR.    &
    667                      ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),2)         .AND.   &                             
    668                              BTEST(wall_flags_total_0(k+1,j,i),2)          .AND.   &
    669                              BTEST(wall_flags_total_0(k,j,i),2) )          .OR.    &
    670                      ( k == nzt .AND. symmetry_flag == 0 ) )                       &
    671                 THEN                                                           
    672                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 15 )     
    673                    flag_set = .TRUE.                                           
    674                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),2)    .OR.     &
    675                            .NOT. BTEST(wall_flags_total_0(k_ppp,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                            .NOT. flag_set                                  .OR.    &
    681                          ( k == nzt - 1 .AND. symmetry_flag == 0 ) )               &
    682                 THEN                                                           
    683                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 )     
    684                    flag_set = .TRUE.                                           
    685                 ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),2)             .AND.   &
    686                          BTEST(wall_flags_total_0(k-1,j,i),2)              .AND.   &
    687                          BTEST(wall_flags_total_0(k,j,i),2)                .AND.   &
    688                          BTEST(wall_flags_total_0(k+1,j,i),2)              .AND.   &
    689                          BTEST(wall_flags_total_0(k_pp,j,i),2)             .AND.   &
    690                          BTEST(wall_flags_total_0(k_ppp,j,i),2)            .AND.   &
    691                          .NOT. flag_set )                                          &
     669                IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),2)         .AND.&
     670                             BTEST(wall_flags_total_0(k,j,i),2)           .AND.&
     671                             BTEST(wall_flags_total_0(k+1,j,i),2) )       .OR. &
     672                     ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),2)        .AND.&
     673                             BTEST(wall_flags_total_0(k+1,j,i),2)         .AND.&
     674                             BTEST(wall_flags_total_0(k,j,i),2) )         .OR. &
     675                     ( k == nzt .AND. symmetry_flag == 0 ) )                   &
     676                THEN
     677                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 15 )
     678                   flag_set = .TRUE.
     679                ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),2)    .OR. &
     680                           .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),2) ) .AND.&
     681                                 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)     .AND.&
     684                                 BTEST(wall_flags_total_0(k_pp,j,i),2)    .AND.&
     685                           .NOT. flag_set                                  .OR.&
     686                         ( k == nzt - 1 .AND. symmetry_flag == 0 ) )           &
     687                THEN
     688                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 )
     689                   flag_set = .TRUE.
     690                ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),2)            .AND.&
     691                         BTEST(wall_flags_total_0(k-1,j,i),2)             .AND.&
     692                         BTEST(wall_flags_total_0(k,j,i),2)               .AND.&
     693                         BTEST(wall_flags_total_0(k+1,j,i),2)             .AND.&
     694                         BTEST(wall_flags_total_0(k_pp,j,i),2)            .AND.&
     695                         BTEST(wall_flags_total_0(k_ppp,j,i),2)           .AND.&
     696                         .NOT. flag_set )                                      &
    692697                THEN
    693698                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 17 )
     
    703708             DO  k = nzb+1, nzt
    704709!
    705 !--             At first, set flags to WS1. 
    706 !--             Since fluxes are swapped in advec_ws.f90, this is necessary to 
     710!--             At first, set flags to WS1.
     711!--             Since fluxes are swapped in advec_ws.f90, this is necessary to
    707712!--             in order to handle the left/south flux.
    708713                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 )
     
    716721                         ( (  bc_dirichlet_r .OR. bc_radiation_r )             &
    717722                           .AND. i == nxr   ) )                                &
    718                 THEN                                                           
    719                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 )     
     723                THEN
     724                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 )
    720725                ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),3)  .AND.   &
    721726                                 BTEST(wall_flags_total_0(k,j,i+1),3)  .OR.    &
     
    726731                         ( ( bc_dirichlet_l .OR.  bc_radiation_l )             &
    727732                           .AND. i == nxlu  ) )                                &
    728                 THEN                                                           
    729                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 19 )     
    730 !                                                                             
    731 !--                Clear flag for WS1                                         
    732                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 )     
     733                THEN
     734                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 19 )
     735!
     736!--                Clear flag for WS1
     737                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 )
    733738                ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),3)  .AND.           &
    734739                         BTEST(wall_flags_total_0(k,j,i+2),3)  .AND.           &
    735740                         BTEST(wall_flags_total_0(k,j,i-1),3) )                &
    736                 THEN                                                           
    737                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i),20 )       
    738 !                                                                             
    739 !--                Clear flag for WS1                                         
    740                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 )     
    741                 ENDIF                                                         
    742 !                                                                             
    743 !--             w component - y-direction                                     
    744 !--             WS1 (21), WS3 (22), WS5 (23)                                   
     741                THEN
     742                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i),20 )
     743!
     744!--                Clear flag for WS1
     745                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 )
     746                ENDIF
     747!
     748!--             w component - y-direction
     749!--             WS1 (21), WS3 (22), WS5 (23)
    745750                IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),3) .OR.           &
    746751                         ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
     
    748753                         ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
    749754                           .AND. j == nyn   ) )                                &
    750                 THEN                                                           
    751                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 )     
     755                THEN
     756                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 )
    752757                ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),3)  .AND.   &
    753758                                 BTEST(wall_flags_total_0(k,j+1,i),3)  .OR.    &
     
    758763                         ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
    759764                           .AND. j == nyn-1 ) )                                &
    760                 THEN                                                           
    761                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 22 )     
    762 !                                                                             
    763 !--                Clear flag for WS1                                         
    764                    advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 )       
     765                THEN
     766                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 22 )
     767!
     768!--                Clear flag for WS1
     769                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 )
    765770                ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),3)  .AND.           &
    766771                         BTEST(wall_flags_total_0(k,j+2,i),3)  .AND.           &
    767772                         BTEST(wall_flags_total_0(k,j-1,i),3) )                &
    768                 THEN                                                           
     773                THEN
    769774                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 23 )
    770775!
     
    772777                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 )
    773778                ENDIF
    774 !                                                                             
     779!
    775780!--             w component - z-direction. Fluxes are calculated on scalar grid
    776 !--             level. Boundary w-values at walls are used. Flux at k=i is 
    777 !--             defined at scalar position k=i+1 with i being an integer. 
     781!--             level. Boundary w-values at walls are used. Flux at k=i is
     782!--             defined at scalar position k=i+1 with i being an integer.
    778783!--             WS1 (24), WS3 (25), WS5 (26)
    779784                IF ( k == nzb+1 )  THEN
     
    784789                IF ( k > nzt-1 )  THEN
    785790                   k_pp = nzt+1
    786                 ELSE 
     791                ELSE
    787792                   k_pp = k + 2
    788793                ENDIF
    789794                IF ( k > nzt-2 )  THEN
    790795                   k_ppp = nzt+1
    791                 ELSE 
     796                ELSE
    792797                   k_ppp = k + 3
    793                 ENDIF 
    794                
    795                 flag_set = .FALSE.                                       
    796                 IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i),3)          .AND.   &
    797                              BTEST(wall_flags_total_0(k+1,j,i),3) )      .OR.    &
    798                      ( .NOT. BTEST(wall_flags_total_0(k+1,j,i),3)        .AND.   &
    799                              BTEST(wall_flags_total_0(k,j,i),3) )        .OR.    &       
    800                      k == nzt -1 )                                               &
     798                ENDIF
     799
     800                flag_set = .FALSE.
     801                IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i),3)          .AND. &
     802                             BTEST(wall_flags_total_0(k+1,j,i),3) )      .OR.  &
     803                     ( .NOT. BTEST(wall_flags_total_0(k+1,j,i),3)        .AND. &
     804                             BTEST(wall_flags_total_0(k,j,i),3) )        .OR.  &
     805                     k == nzt -1 )                                             &
    801806                THEN
    802807!
    803808!--                Please note, at k == nzb_w_inner(j,i) a flag is explicitly
    804 !--                set, although this is not a prognostic level. However, 
     809!--                set, although this is not a prognostic level. However,
    805810!--                contrary to the advection of u,v and s this is necessary
    806811!--                because flux_t(nzb_w_inner(j,i)) is used for the tendency
     
    808813                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 24 )
    809814                   flag_set = .TRUE.
    810                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),3)    .AND.   &
    811                                  BTEST(wall_flags_total_0(k,j,i),3)      .AND.   &
    812                                  BTEST(wall_flags_total_0(k+1,j,i),3)    .AND.   &
    813                                  BTEST(wall_flags_total_0(k_pp,j,i),3) ) .OR.    &
    814                          ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),3)   .AND.   &
    815                                  BTEST(wall_flags_total_0(k+1,j,i),3)    .AND.   &
    816                                  BTEST(wall_flags_total_0(k,j,i),3)      .AND.   &
    817                                  BTEST(wall_flags_total_0(k-1,j,i),3) )  .AND.   &
    818                            .NOT. flag_set                          .OR.          &
    819                          k == nzt - 2 )                                          &
    820                 THEN                                                           
    821                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 25 )     
    822                    flag_set = .TRUE.                                           
    823                 ELSEIF ( BTEST(wall_flags_total_0(k-1,j,i),3)            .AND.   &
    824                          BTEST(wall_flags_total_0(k,j,i),3)              .AND.   &
    825                          BTEST(wall_flags_total_0(k+1,j,i),3)            .AND.   &
    826                          BTEST(wall_flags_total_0(k_pp,j,i),3)           .AND.   &
    827                          .NOT. flag_set )                                        &
    828                 THEN                                                               
     815                ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),3)    .AND. &
     816                                 BTEST(wall_flags_total_0(k,j,i),3)      .AND. &
     817                                 BTEST(wall_flags_total_0(k+1,j,i),3)    .AND. &
     818                                 BTEST(wall_flags_total_0(k_pp,j,i),3) ) .OR.  &
     819                         ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),3)   .AND. &
     820                                 BTEST(wall_flags_total_0(k+1,j,i),3)    .AND. &
     821                                 BTEST(wall_flags_total_0(k,j,i),3)      .AND. &
     822                                 BTEST(wall_flags_total_0(k-1,j,i),3) )  .AND. &
     823                           .NOT. flag_set                          .OR.        &
     824                         k == nzt - 2 )                                        &
     825                THEN
     826                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 25 )
     827                   flag_set = .TRUE.
     828                ELSEIF ( BTEST(wall_flags_total_0(k-1,j,i),3)            .AND. &
     829                         BTEST(wall_flags_total_0(k,j,i),3)              .AND. &
     830                         BTEST(wall_flags_total_0(k+1,j,i),3)            .AND. &
     831                         BTEST(wall_flags_total_0(k_pp,j,i),3)           .AND. &
     832                         .NOT. flag_set )                                      &
     833                THEN
    829834                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 26 )
    830835                ENDIF
     
    837842       CALL exchange_horiz_int( advc_flags_m, nys, nyn, nxl, nxr, nzt, nbgp )
    838843!
    839 !--    Set boundary flags at inflow and outflow boundary in case of 
     844!--    Set boundary flags at inflow and outflow boundary in case of
    840845!--    non-cyclic boundary conditions.
    841846       IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
     
    844849
    845850       IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
    846          advc_flags_m(:,:,nxr+1) = advc_flags_m(:,:,nxr)
     851          advc_flags_m(:,:,nxr+1) = advc_flags_m(:,:,nxr)
    847852       ENDIF
    848853
     
    863868!> Initialization of flags to control the order of the advection scheme near
    864869!> solid walls and non-cyclic inflow boundaries, where the order is sucessively
    865 !> degraded. 
     870!> degraded.
    866871!------------------------------------------------------------------------------!
    867872    SUBROUTINE ws_init_flags_scalar( non_cyclic_l, non_cyclic_n, non_cyclic_r, &
     
    885890       LOGICAL ::  non_cyclic_s !< flag that indicates non-cyclic boundary on the south
    886891
    887        LOGICAL, OPTIONAL ::  extensive_degrad !< flag indicating that extensive degradation is required, e.g. for 
     892       LOGICAL, OPTIONAL ::  extensive_degrad !< flag indicating that extensive degradation is required, e.g. for
    888893                                              !< passive scalars nearby topography along the horizontal directions,
    889                                               !< as no monotonic limiter can be applied there 
     894                                              !< as no monotonic limiter can be applied there
    890895!
    891896!--    Set flags to steer the degradation of the advection scheme in advec_ws
     
    900905!--             scalar - x-direction
    901906!--             WS1 (0), WS3 (1), WS5 (2)
    902                 IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),0)      .OR.      &
    903                        .NOT. BTEST(wall_flags_total_0(k,j,i+2),0)      .OR.      &   
    904                        .NOT. BTEST(wall_flags_total_0(k,j,i-1),0) )    .OR.      &
    905                        ( non_cyclic_l  .AND.  i == 0  )                .OR.      &
    906                        ( non_cyclic_r  .AND.  i == nx ) )  THEN           
    907                  advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )             
    908                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+3),0)  .AND.     &
    909                                  BTEST(wall_flags_total_0(k,j,i+1),0)  .AND.     &
    910                                  BTEST(wall_flags_total_0(k,j,i+2),0)  .AND.     &
    911                                  BTEST(wall_flags_total_0(k,j,i-1),0)            &
    912                          )                       .OR.                            &
    913                          ( .NOT. BTEST(wall_flags_total_0(k,j,i-2),0)  .AND.     &
    914                                  BTEST(wall_flags_total_0(k,j,i+1),0)  .AND.     &
    915                                  BTEST(wall_flags_total_0(k,j,i+2),0)  .AND.     &
    916                                  BTEST(wall_flags_total_0(k,j,i-1),0)            &
    917                          )                                                       &
    918                                                  .OR.                            &
    919                          ( non_cyclic_r  .AND.  i == nx-1 )  .OR.                &
    920                          ( non_cyclic_l  .AND.  i == 1    ) )  THEN           
    921                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )             
    922                 ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),0)           .AND.    &
    923                          BTEST(wall_flags_total_0(k,j,i+2),0)           .AND.    &
    924                          BTEST(wall_flags_total_0(k,j,i+3),0)           .AND.    &
    925                          BTEST(wall_flags_total_0(k,j,i-1),0)           .AND.    &
    926                          BTEST(wall_flags_total_0(k,j,i-2),0) )                  &
    927                 THEN                                                           
    928                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 )             
    929                 ENDIF                                                         
    930 !                                                                             
    931 !--             scalar - y-direction                                           
    932 !--             WS1 (3), WS3 (4), WS5 (5)                                     
    933                 IF ( ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),0)        .OR.    &
    934                        .NOT. BTEST(wall_flags_total_0(k,j+2,i),0)        .OR.    &   
    935                        .NOT. BTEST(wall_flags_total_0(k,j-1,i),0))       .OR.    &
    936                      ( non_cyclic_s  .AND.  j == 0  )                    .OR.    &
    937                      ( non_cyclic_n  .AND.  j == ny ) )  THEN                                                           
    938                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )             
    939 !                                                                             
    940 !--             WS3                                                           
    941                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+3,i),0)    .AND.   &
    942                                  BTEST(wall_flags_total_0(k,j+1,i),0)    .AND.   &
    943                                  BTEST(wall_flags_total_0(k,j+2,i),0)    .AND.   &
    944                                  BTEST(wall_flags_total_0(k,j-1,i),0)            &
    945                          )                       .OR.                            &
    946                          ( .NOT. BTEST(wall_flags_total_0(k,j-2,i),0)    .AND.   &
    947                                  BTEST(wall_flags_total_0(k,j+1,i),0)    .AND.   &
    948                                  BTEST(wall_flags_total_0(k,j+2,i),0)    .AND.   &
    949                                  BTEST(wall_flags_total_0(k,j-1,i),0)            &
    950                          )                                                       &
    951                                                  .OR.                            &
    952                          ( non_cyclic_s  .AND.  j == 1    )  .OR.                &
    953                          ( non_cyclic_n  .AND.  j == ny-1 ) )  THEN           
    954                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )             
    955 !                                                                             
    956 !--             WS5                                                           
    957                 ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),0)           .AND.    &
    958                          BTEST(wall_flags_total_0(k,j+2,i),0)           .AND.    &
    959                          BTEST(wall_flags_total_0(k,j+3,i),0)           .AND.    &
    960                          BTEST(wall_flags_total_0(k,j-1,i),0)           .AND.    &
    961                          BTEST(wall_flags_total_0(k,j-2,i),0) )                  &
    962                 THEN                                                           
    963                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 )             
     907                IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),0)      .OR.    &
     908                       .NOT. BTEST(wall_flags_total_0(k,j,i+2),0)      .OR.    &
     909                       .NOT. BTEST(wall_flags_total_0(k,j,i-1),0) )    .OR.    &
     910                       ( non_cyclic_l  .AND.  i == 0  )                .OR.    &
     911                       ( non_cyclic_r  .AND.  i == nx ) )  THEN
     912                 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
     913                ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+3),0)  .AND.   &
     914                                 BTEST(wall_flags_total_0(k,j,i+1),0)  .AND.   &
     915                                 BTEST(wall_flags_total_0(k,j,i+2),0)  .AND.   &
     916                                 BTEST(wall_flags_total_0(k,j,i-1),0)          &
     917                         )                       .OR.                          &
     918                         ( .NOT. BTEST(wall_flags_total_0(k,j,i-2),0)  .AND.   &
     919                                 BTEST(wall_flags_total_0(k,j,i+1),0)  .AND.   &
     920                                 BTEST(wall_flags_total_0(k,j,i+2),0)  .AND.   &
     921                                 BTEST(wall_flags_total_0(k,j,i-1),0)          &
     922                         )                                                     &
     923                                                 .OR.                          &
     924                         ( non_cyclic_r  .AND.  i == nx-1 )  .OR.              &
     925                         ( non_cyclic_l  .AND.  i == 1    ) )  THEN
     926                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
     927                ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),0)           .AND.  &
     928                         BTEST(wall_flags_total_0(k,j,i+2),0)           .AND.  &
     929                         BTEST(wall_flags_total_0(k,j,i+3),0)           .AND.  &
     930                         BTEST(wall_flags_total_0(k,j,i-1),0)           .AND.  &
     931                         BTEST(wall_flags_total_0(k,j,i-2),0) )                &
     932                THEN
     933                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 )
    964934                ENDIF
    965935!
    966 !--             Near topography, set horizontal advection scheme to 1st order
    967 !--             for passive scalars, even if only one direction may be
     936!--             scalar - y-direction
     937!--             WS1 (3), WS3 (4), WS5 (5)
     938                IF ( ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),0)        .OR.  &
     939                       .NOT. BTEST(wall_flags_total_0(k,j+2,i),0)        .OR.  &
     940                       .NOT. BTEST(wall_flags_total_0(k,j-1,i),0))       .OR.  &
     941                     ( non_cyclic_s  .AND.  j == 0  )                    .OR.  &
     942                     ( non_cyclic_n  .AND.  j == ny ) )  THEN
     943                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
     944!
     945!--             WS3
     946                ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+3,i),0)    .AND. &
     947                                 BTEST(wall_flags_total_0(k,j+1,i),0)    .AND. &
     948                                 BTEST(wall_flags_total_0(k,j+2,i),0)    .AND. &
     949                                 BTEST(wall_flags_total_0(k,j-1,i),0)          &
     950                         )                       .OR.                          &
     951                         ( .NOT. BTEST(wall_flags_total_0(k,j-2,i),0)    .AND. &
     952                                 BTEST(wall_flags_total_0(k,j+1,i),0)    .AND. &
     953                                 BTEST(wall_flags_total_0(k,j+2,i),0)    .AND. &
     954                                 BTEST(wall_flags_total_0(k,j-1,i),0)          &
     955                         )                                                     &
     956                                                 .OR.                          &
     957                         ( non_cyclic_s  .AND.  j == 1    )  .OR.              &
     958                         ( non_cyclic_n  .AND.  j == ny-1 ) )  THEN
     959                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
     960!
     961!--             WS5
     962                ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),0)           .AND.  &
     963                         BTEST(wall_flags_total_0(k,j+2,i),0)           .AND.  &
     964                         BTEST(wall_flags_total_0(k,j+3,i),0)           .AND.  &
     965                         BTEST(wall_flags_total_0(k,j-1,i),0)           .AND.  &
     966                         BTEST(wall_flags_total_0(k,j-2,i),0) )                &
     967                THEN
     968                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 )
     969                ENDIF
     970!
     971!--             Near topography, set horizontal advection scheme to 1st order
     972!--             for passive scalars, even if only one direction may be
    968973!--             blocked by topography. These locations will be identified
    969 !--             by wall_flags_total_0 bit 31. Note, since several modules define 
    970 !--             advection flags but may apply different scalar boundary 
    971 !--             conditions, bit 31 is temporarily stored on advc_flags. 
    972 !--             Moreover, note that this extended degradtion for passive 
    973 !--             scalars is not required for the vertical direction as there 
    974 !--             the monotonic limiter can be applied. 
     974!--             by wall_flags_total_0 bit 31. Note, since several modules define
     975!--             advection flags but may apply different scalar boundary
     976!--             conditions, bit 31 is temporarily stored on advc_flags.
     977!--             Moreover, note that this extended degradtion for passive
     978!--             scalars is not required for the vertical direction as there
     979!--             the monotonic limiter can be applied.
    975980                IF ( PRESENT( extensive_degrad ) )  THEN
    976981                   IF ( extensive_degrad )  THEN
     
    980985                      IF( BTEST( advc_flag(k,j,i), 31 ) )  THEN
    981986!
    982 !--                      Clear flags that might indicate higher-order 
     987!--                      Clear flags that might indicate higher-order
    983988!--                      advection along x- and y-direction.
    984989                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
     
    990995!--                      x- and y-direction.
    991996                         advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
    992                          advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) 
     997                         advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
    993998                      ENDIF
    994999!
    9951000!--                   Adjacent to this extended degradation zone, successively
    996 !--                   upgrade the order of the scheme if this grid point isn't 
    997 !--                   flagged with bit 31 (indicating extended degradation 
    998 !--                   zone). 
     1001!--                   upgrade the order of the scheme if this grid point isn't
     1002!--                   flagged with bit 31 (indicating extended degradation
     1003!--                   zone).
    9991004                      IF ( .NOT. BTEST( advc_flag(k,j,i), 31 ) )  THEN
    10001005!
    10011006!--                      x-direction. First, clear all previous settings, than
    10021007!--                      set flag for 3rd-order scheme.
    1003                          IF ( BTEST( advc_flag(k,j,i-1), 31 )  .AND.        &
     1008                         IF ( BTEST( advc_flag(k,j,i-1), 31 )  .AND.           &
    10041009                              BTEST( advc_flag(k,j,i+1), 31 ) )  THEN
    10051010                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
    10061011                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
    10071012                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
    1008                            
     1013
    10091014                            advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
    10101015                         ENDIF
     
    10121017!--                      x-direction. First, clear all previous settings, than
    10131018!--                      set flag for 5rd-order scheme.
    1014                          IF ( .NOT. BTEST( advc_flag(k,j,i-1), 31 )  .AND.  &
    1015                                     BTEST( advc_flag(k,j,i-2), 31 )  .AND.  &
    1016                               .NOT. BTEST( advc_flag(k,j,i+1), 31 )  .AND.  &
     1019                         IF ( .NOT. BTEST( advc_flag(k,j,i-1), 31 )  .AND.     &
     1020                                    BTEST( advc_flag(k,j,i-2), 31 )  .AND.     &
     1021                              .NOT. BTEST( advc_flag(k,j,i+1), 31 )  .AND.     &
    10171022                                    BTEST( advc_flag(k,j,i+2), 31 ) )  THEN
    10181023                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
    10191024                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
    10201025                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
    1021                            
     1026
    10221027                            advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 )
    10231028                         ENDIF
     
    10251030!--                      y-direction. First, clear all previous settings, than
    10261031!--                      set flag for 3rd-order scheme.
    1027                          IF ( BTEST( advc_flag(k,j-1,i), 31 )  .AND.        &
     1032                         IF ( BTEST( advc_flag(k,j-1,i), 31 )  .AND.           &
    10281033                              BTEST( advc_flag(k,j+1,i), 31 ) )  THEN
    10291034                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
    10301035                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
    10311036                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
    1032                            
     1037
    10331038                            advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
    10341039                         ENDIF
     
    10361041!--                      y-direction. First, clear all previous settings, than
    10371042!--                      set flag for 5rd-order scheme.
    1038                          IF ( .NOT. BTEST( advc_flag(k,j-1,i), 31 )  .AND.  &
    1039                                     BTEST( advc_flag(k,j-2,i), 31 )  .AND.  &
    1040                               .NOT. BTEST( advc_flag(k,j+1,i), 31 )  .AND.  &
     1043                         IF ( .NOT. BTEST( advc_flag(k,j-1,i), 31 )  .AND.     &
     1044                                    BTEST( advc_flag(k,j-2,i), 31 )  .AND.     &
     1045                              .NOT. BTEST( advc_flag(k,j+1,i), 31 )  .AND.     &
    10411046                                    BTEST( advc_flag(k,j+2,i), 31 ) )  THEN
    10421047                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
    10431048                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
    10441049                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
    1045                            
     1050
    10461051                            advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 )
    10471052                         ENDIF
    10481053                      ENDIF
    1049                    
     1054
    10501055                   ENDIF
    1051                    
     1056
    10521057!
    10531058!--                Near lateral boundary flags might be overwritten. Set
    1054 !--                them again. 
     1059!--                them again.
    10551060!--                x-direction
    1056                    IF ( ( non_cyclic_l  .AND.  i == 0  )  .OR.             &
     1061                   IF ( ( non_cyclic_l  .AND.  i == 0  )  .OR.                 &
    10571062                        ( non_cyclic_r  .AND.  i == nx ) )  THEN
    10581063                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
    10591064                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
    10601065                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
    1061                          
     1066
    10621067                      advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
    10631068                   ENDIF
    1064                    
    1065                    IF ( ( non_cyclic_l  .AND.  i == 1    )  .OR.           &
     1069
     1070                   IF ( ( non_cyclic_l  .AND.  i == 1    )  .OR.               &
    10661071                        ( non_cyclic_r  .AND.  i == nx-1 ) )  THEN
    10671072                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
    10681073                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
    10691074                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
    1070                          
     1075
    10711076                      advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
    10721077                   ENDIF
    10731078!
    10741079!--                y-direction
    1075                    IF ( ( non_cyclic_n  .AND.  j == 0  )  .OR.             &
     1080                   IF ( ( non_cyclic_n  .AND.  j == 0  )  .OR.                 &
    10761081                        ( non_cyclic_s  .AND.  j == ny ) )  THEN
    10771082                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
    10781083                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
    10791084                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
    1080                          
     1085
    10811086                      advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
    10821087                   ENDIF
    1083                    
    1084                    IF ( ( non_cyclic_n  .AND.  j == 1    )  .OR.           &
     1088
     1089                   IF ( ( non_cyclic_n  .AND.  j == 1    )  .OR.               &
    10851090                        ( non_cyclic_s  .AND.  j == ny-1 ) )  THEN
    10861091                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
    10871092                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
    10881093                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
    1089                          
     1094
    10901095                      advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
    10911096                   ENDIF
    1092                    
     1097
    10931098                ENDIF
    1094                
    1095                
    1096 !                                                                             
    1097 !--             scalar - z-direction. Fluxes are calculated on w-grid                                     
    1098 !--             level. Boundary values at/within walls aren't used.                                           
    1099 !--             WS1 (6), WS3 (7), WS5 (8)                                     
    1100                 IF ( k == nzb+1 )  THEN                                       
    1101                    k_mm = nzb                                                 
    1102                 ELSE                                                           
    1103                    k_mm = k - 2                                               
    1104                 ENDIF                                                         
    1105                 IF ( k > nzt-1 )  THEN                                         
    1106                    k_pp = nzt+1                                               
    1107                 ELSE                                                           
    1108                    k_pp = k + 2                                               
    1109                 ENDIF                                                         
    1110                 IF ( k > nzt-2 )  THEN                                         
    1111                    k_ppp = nzt+1                                               
    1112                 ELSE                                                           
    1113                    k_ppp = k + 3                                               
    1114                 ENDIF                                                         
    1115                                                                                
    1116                 flag_set = .FALSE.                                             
    1117                 IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),0)          .AND.  &
    1118                              BTEST(wall_flags_total_0(k,j,i),0)            .AND.  &
    1119                              BTEST(wall_flags_total_0(k+1,j,i),0) )        .OR.   &
    1120                      ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),0)         .AND.  &                             
    1121                              BTEST(wall_flags_total_0(k+1,j,i),0)          .AND.  &
    1122                              BTEST(wall_flags_total_0(k,j,i),0) )          .OR.   &
    1123                      ( k == nzt .AND. symmetry_flag == 0 ) )                      &
    1124                 THEN                                                           
    1125                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 6 )             
    1126                    flag_set = .TRUE.                                           
    1127                 ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),0)    .OR.    &
    1128                            .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),0) ) .AND.   &
    1129                                  BTEST(wall_flags_total_0(k-1,j,i),0)     .AND.   &
    1130                                  BTEST(wall_flags_total_0(k,j,i),0)       .AND.   &
    1131                                  BTEST(wall_flags_total_0(k+1,j,i),0)     .AND.   &
    1132                                  BTEST(wall_flags_total_0(k_pp,j,i),0)    .AND.   &
    1133                            .NOT. flag_set                                  .OR.   &
    1134                          ( k == nzt - 1 .AND. symmetry_flag == 0 ) )              &
    1135                 THEN                                                           
    1136                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 7 )             
    1137                    flag_set = .TRUE.                                           
    1138                 ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),0)            .AND.   &
    1139                          BTEST(wall_flags_total_0(k-1,j,i),0)             .AND.   &
    1140                          BTEST(wall_flags_total_0(k,j,i),0)               .AND.   &
    1141                          BTEST(wall_flags_total_0(k+1,j,i),0)             .AND.   &
    1142                          BTEST(wall_flags_total_0(k_pp,j,i),0)            .AND.   &
    1143                          BTEST(wall_flags_total_0(k_ppp,j,i),0)           .AND.   &
    1144                         .NOT. flag_set )                                          &
     1099
     1100
     1101!
     1102!--             scalar - z-direction. Fluxes are calculated on w-grid
     1103!--             level. Boundary values at/within walls aren't used.
     1104!--             WS1 (6), WS3 (7), WS5 (8)
     1105                IF ( k == nzb+1 )  THEN
     1106                   k_mm = nzb
     1107                ELSE
     1108                   k_mm = k - 2
     1109                ENDIF
     1110                IF ( k > nzt-1 )  THEN
     1111                   k_pp = nzt+1
     1112                ELSE
     1113                   k_pp = k + 2
     1114                ENDIF
     1115                IF ( k > nzt-2 )  THEN
     1116                   k_ppp = nzt+1
     1117                ELSE
     1118                   k_ppp = k + 3
     1119                ENDIF
     1120
     1121                flag_set = .FALSE.
     1122                IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),0)       .AND.  &
     1123                             BTEST(wall_flags_total_0(k,j,i),0)         .AND.  &
     1124                             BTEST(wall_flags_total_0(k+1,j,i),0) )     .OR.   &
     1125                     ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),0)      .AND.  &
     1126                             BTEST(wall_flags_total_0(k+1,j,i),0)       .AND.  &
     1127                             BTEST(wall_flags_total_0(k,j,i),0) )       .OR.   &
     1128                     ( k == nzt .AND. symmetry_flag == 0 ) )                   &
     1129                THEN
     1130                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 6 )
     1131                   flag_set = .TRUE.
     1132                ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),0)    .OR. &
     1133                           .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),0) ) .AND.&
     1134                                 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)     .AND.&
     1137                                 BTEST(wall_flags_total_0(k_pp,j,i),0)    .AND.&
     1138                           .NOT. flag_set                                  .OR.&
     1139                         ( k == nzt - 1 .AND. symmetry_flag == 0 ) )           &
     1140                THEN
     1141                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 7 )
     1142                   flag_set = .TRUE.
     1143                ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),0)         .AND.   &
     1144                         BTEST(wall_flags_total_0(k-1,j,i),0)          .AND.   &
     1145                         BTEST(wall_flags_total_0(k,j,i),0)            .AND.   &
     1146                         BTEST(wall_flags_total_0(k+1,j,i),0)          .AND.   &
     1147                         BTEST(wall_flags_total_0(k_pp,j,i),0)         .AND.   &
     1148                         BTEST(wall_flags_total_0(k_ppp,j,i),0)        .AND.   &
     1149                        .NOT. flag_set )                                       &
    11451150                THEN
    11461151                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 8 )
     
    11561161       CALL exchange_horiz_int( advc_flag, nys, nyn, nxl, nxr, nzt, nbgp )
    11571162!
    1158 !--    Set boundary flags at inflow and outflow boundary in case of 
     1163!--    Set boundary flags at inflow and outflow boundary in case of
    11591164!--    non-cyclic boundary conditions.
    11601165       IF ( non_cyclic_l )  THEN
     
    11731178          advc_flag(:,nys-1,:) = advc_flag(:,nys,:)
    11741179       ENDIF
    1175  
    1176 
    1177 
    1178     END SUBROUTINE ws_init_flags_scalar   
    1179    
     1180
     1181
     1182
     1183    END SUBROUTINE ws_init_flags_scalar
     1184
    11801185!------------------------------------------------------------------------------!
    11811186! Description:
     
    11871192
    11881193!
    1189 !--    The arrays needed for statistical evaluation are set to to 0 at the 
     1194!--    The arrays needed for statistical evaluation are set to to 0 at the
    11901195!--    beginning of prognostic_equations.
    11911196       IF ( ws_scheme_mom )  THEN
     
    12251230
    12261231
    1227        CHARACTER (LEN = *), INTENT(IN) ::  sk_char !< string identifier, used for assign fluxes to the correct dimension in the analysis array
    1228        
     1232       CHARACTER (LEN = *), INTENT(IN) ::  sk_char !< string identifier, used for assign fluxes to the
     1233                                                   !<correct dimension in the analysis array
     1234
    12291235       INTEGER(iwp) ::  i         !< grid index along x-direction
    12301236       INTEGER(iwp) ::  i_omp     !< leftmost index on subdomain, or in case of OpenMP, on thread
     
    12351241       INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
    12361242       INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
    1237        INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 
     1243       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
    12381244       INTEGER(iwp) ::  tn        !< number of OpenMP thread
    1239        
     1245
    12401246       INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::   &
    12411247                                                  advc_flag !< flag array to control order of scalar advection
     
    12441250       LOGICAL           ::  non_cyclic_n    !< flag that indicates non-cyclic boundary on the north
    12451251       LOGICAL           ::  non_cyclic_r    !< flag that indicates non-cyclic boundary on the right
    1246        LOGICAL           ::  non_cyclic_s    !< flag that indicates non-cyclic boundary on the south                                           
     1252       LOGICAL           ::  non_cyclic_s    !< flag that indicates non-cyclic boundary on the south
    12471253       LOGICAL, OPTIONAL ::  flux_limitation !< flag indicating flux limitation of the vertical advection
    12481254       LOGICAL           ::  limiter         !< control flag indicating the application of flux limitation
     
    12511257       REAL(wp) ::  div           !< velocity diverence on scalar grid
    12521258       REAL(wp) ::  div_in        !< vertical flux divergence of ingoing fluxes
    1253        REAL(wp) ::  div_out       !< vertical flux divergence of outgoing fluxes     
     1259       REAL(wp) ::  div_out       !< vertical flux divergence of outgoing fluxes
    12541260       REAL(wp) ::  f_corr_t      !< correction flux at grid-cell top, i.e. the difference between high and low-order flux
    12551261       REAL(wp) ::  f_corr_d      !< correction flux at grid-cell bottom, i.e. the difference between high and low-order flux
     
    12721278       REAL(wp) ::  min_val       !< maximum value of the quanitity along the numerical stencil (in vertical direction)
    12731279       REAL(wp) ::  mon           !< monotone solution of the advection equation using 1st-order fluxes
    1274        REAL(wp) ::  u_comp        !< advection velocity along x-direction 
    1275        REAL(wp) ::  v_comp        !< advection velocity along y-direction 
    1276 !       
    1277 !--    sk is an array from parameter list. It should not be a pointer, because 
    1278 !--    in that case the compiler can not assume a stride 1 and cannot perform 
    1279 !--    a strided one vector load. Adding the CONTIGUOUS keyword makes things 
    1280 !--    even worse, because the compiler cannot assume strided one in the 
     1280       REAL(wp) ::  u_comp        !< advection velocity along x-direction
     1281       REAL(wp) ::  v_comp        !< advection velocity along y-direction
     1282!
     1283!--    sk is an array from parameter list. It should not be a pointer, because
     1284!--    in that case the compiler can not assume a stride 1 and cannot perform
     1285!--    a strided one vector load. Adding the CONTIGUOUS keyword makes things
     1286!--    even worse, because the compiler cannot assume strided one in the
    12811287!--    caller side.
    12821288       REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !<  advected scalar
    12831289
    1284        REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_n     !< discretized artificial dissipation at northward-side of the grid box
    1285        REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_r     !< discretized artificial dissipation at rightward-side of the grid box
    1286        REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_t     !< discretized artificial dissipation at top of the grid box
    1287        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_n     !< discretized 6th-order flux at northward-side of the grid box
    1288        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_r     !< discretized 6th-order flux at rightward-side of the grid box
    1289        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t     !< discretized 6th-order flux at top of the grid box
    1290        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t_1st !< discretized 1st-order flux at top of the grid box
    1291        
    1292        REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_diss_y_local !< discretized artificial dissipation at southward-side of the grid box
    1293        REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_flux_y_local !< discretized 6th-order flux at northward-side of the grid box
    1294        
    1295        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_diss_x_local !< discretized artificial dissipation at leftward-side of the grid box
    1296        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_flux_x_local !< discretized 6th-order flux at leftward-side of the grid box
    1297 !
    1298 !--    Used local modified copy of nzb_max (used to degrade order of
     1290       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_n     !< discretized artificial dissipation at northward-side
     1291       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_r     !< discretized artificial dissipation at rightward-side
     1292       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_t     !< discretized artificial dissipation at top
     1293       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_n     !< discretized 6th-order flux at northward-side
     1294       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_r     !< discretized 6th-order flux at rightward-side
     1295       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t     !< discretized 6th-order flux at top
     1296       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t_1st !< discretized 1st-order flux at top
     1297
     1298       REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_diss_y_local !< discretized artificial dissipation at southward-side
     1299       REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_flux_y_local !< discretized 6th-order flux at northward-side
     1300       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_diss_x_local !< discretized artificial dissipation at leftward-side
     1301       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::  swap_flux_x_local !< discretized 6th-order flux at leftward-side
     1302!
     1303!--    Used local modified copy of nzb_max (used to degrade order of
    12991304!--    discretization) at non-cyclic boundaries. Modify only at relevant points
    1300 !--    instead of the entire subdomain. This should lead to better 
     1305!--    instead of the entire subdomain. This should lead to better
    13011306!--    load balance between boundary and non-boundary PEs.
    13021307       IF( non_cyclic_l  .AND.  i <= nxl + 2  .OR.                             &
     
    13091314       END IF
    13101315!
    1311 !--    Set control flag for flux limiter 
     1316!--    Set control flag for flux limiter
    13121317       limiter = .FALSE.
    13131318       IF ( PRESENT( flux_limitation) )  limiter = flux_limitation
     
    13241329
    13251330             v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
    1326              swap_flux_y_local(k,tn) = v_comp *         (                     &
    1327                                                ( 37.0_wp * ibit5 * adv_sca_5  &
    1328                                             +     7.0_wp * ibit4 * adv_sca_3  &
    1329                                             +              ibit3 * adv_sca_1  &
    1330                                                ) *                            &
    1331                                            ( sk(k,j,i)  + sk(k,j-1,i)     )   &
    1332                                          -     (  8.0_wp * ibit5 * adv_sca_5  &
    1333                                             +              ibit4 * adv_sca_3  &
    1334                                                 ) *                           &
    1335                                            ( sk(k,j+1,i) + sk(k,j-2,i)    )   &
    1336                                          +     (           ibit5 * adv_sca_5  &
    1337                                                ) *                            &
    1338                                            ( sk(k,j+2,i) + sk(k,j-3,i)    )   &
     1331             swap_flux_y_local(k,tn) = v_comp *         (                      &
     1332                                               ( 37.0_wp * ibit5 * adv_sca_5   &
     1333                                            +     7.0_wp * ibit4 * adv_sca_3   &
     1334                                            +              ibit3 * adv_sca_1   &
     1335                                               ) *                             &
     1336                                           ( sk(k,j,i)  + sk(k,j-1,i)     )    &
     1337                                         -     (  8.0_wp * ibit5 * adv_sca_5   &
     1338                                            +              ibit4 * adv_sca_3   &
     1339                                                ) *                            &
     1340                                           ( sk(k,j+1,i) + sk(k,j-2,i)    )    &
     1341                                         +     (           ibit5 * adv_sca_5   &
     1342                                               ) *                             &
     1343                                           ( sk(k,j+2,i) + sk(k,j-3,i)    )    &
    13391344                                                        )
    13401345
    1341              swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                     &
    1342                                                ( 10.0_wp * ibit5 * adv_sca_5  &
    1343                                             +     3.0_wp * ibit4 * adv_sca_3  &
    1344                                             +              ibit3 * adv_sca_1  &
    1345                                                ) *                            &
    1346                                             ( sk(k,j,i)   - sk(k,j-1,i)  )    &
    1347                                         -      (  5.0_wp * ibit5 * adv_sca_5  &
    1348                                             +              ibit4 * adv_sca_3  &
    1349                                             ) *                               &
    1350                                             ( sk(k,j+1,i) - sk(k,j-2,i)  )    &
    1351                                         +      (           ibit5 * adv_sca_5  &
    1352                                                ) *                            &
    1353                                             ( sk(k,j+2,i) - sk(k,j-3,i)  )    &
     1346             swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                      &
     1347                                               ( 10.0_wp * ibit5 * adv_sca_5   &
     1348                                            +     3.0_wp * ibit4 * adv_sca_3   &
     1349                                            +              ibit3 * adv_sca_1   &
     1350                                               ) *                             &
     1351                                            ( sk(k,j,i)   - sk(k,j-1,i)  )     &
     1352                                        -      (  5.0_wp * ibit5 * adv_sca_5   &
     1353                                            +              ibit4 * adv_sca_3   &
     1354                                            ) *                                &
     1355                                            ( sk(k,j+1,i) - sk(k,j-2,i)  )     &
     1356                                        +      (           ibit5 * adv_sca_5   &
     1357                                               ) *                             &
     1358                                            ( sk(k,j+2,i) - sk(k,j-3,i)  )     &
    13541359                                                        )
    13551360
     
    13601365
    13611366             v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
    1362              swap_flux_y_local(k,tn) = v_comp * (                             &
    1363                                     37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )   &
    1364                                   -  8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) )   &
    1365                                   +           ( sk(k,j+2,i) + sk(k,j-3,i) )   &
     1367             swap_flux_y_local(k,tn) = v_comp * (                              &
     1368                                    37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )    &
     1369                                  -  8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) )    &
     1370                                  +           ( sk(k,j+2,i) + sk(k,j-3,i) )    &
    13661371                                                ) * adv_sca_5
    1367              swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                     &
    1368                                     10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )   &
    1369                                   -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )   &
    1370                                   +             sk(k,j+2,i) - sk(k,j-3,i)     &
     1372             swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                      &
     1373                                    10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )    &
     1374                                  -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )    &
     1375                                  +             sk(k,j+2,i) - sk(k,j-3,i)      &
    13711376                                                        ) * adv_sca_5
    13721377
     
    13771382!--    Compute leftside fluxes of the respective PE bounds.
    13781383       IF ( i == i_omp )  THEN
    1379        
     1384
    13801385          DO  k = nzb+1, nzb_max_l
    13811386
     
    13841389             ibit0 = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )
    13851390
    1386              u_comp                     = u(k,j,i) - u_gtrans + u_stokes_zu(k)
    1387              swap_flux_x_local(k,j,tn) = u_comp * (                           &
    1388                                                ( 37.0_wp * ibit2 * adv_sca_5  &
    1389                                             +     7.0_wp * ibit1 * adv_sca_3  &
    1390                                             +              ibit0 * adv_sca_1  &
    1391                                                ) *                            &
    1392                                             ( sk(k,j,i)   + sk(k,j,i-1)    )  &
    1393                                         -      (  8.0_wp * ibit2 * adv_sca_5  &
    1394                                             +              ibit1 * adv_sca_3  &
    1395                                                ) *                            &
    1396                                             ( sk(k,j,i+1) + sk(k,j,i-2)    )  &
    1397                                         +      (           ibit2 * adv_sca_5  &
    1398                                                ) *                            &
    1399                                             ( sk(k,j,i+2) + sk(k,j,i-3)    )  &
     1391             u_comp                    = u(k,j,i) - u_gtrans + u_stokes_zu(k)
     1392             swap_flux_x_local(k,j,tn) = u_comp * (                            &
     1393                                               ( 37.0_wp * ibit2 * adv_sca_5   &
     1394                                            +     7.0_wp * ibit1 * adv_sca_3   &
     1395                                            +              ibit0 * adv_sca_1   &
     1396                                               ) *                             &
     1397                                            ( sk(k,j,i)   + sk(k,j,i-1)    )   &
     1398                                        -      (  8.0_wp * ibit2 * adv_sca_5   &
     1399                                            +              ibit1 * adv_sca_3   &
     1400                                               ) *                             &
     1401                                            ( sk(k,j,i+1) + sk(k,j,i-2)    )   &
     1402                                        +      (           ibit2 * adv_sca_5   &
     1403                                               ) *                             &
     1404                                            ( sk(k,j,i+2) + sk(k,j,i-3)    )   &
    14001405                                                  )
    14011406
    1402              swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                   &
    1403                                                ( 10.0_wp * ibit2 * adv_sca_5  &
    1404                                             +     3.0_wp * ibit1 * adv_sca_3  &
    1405                                             +              ibit0 * adv_sca_1  &
    1406                                                ) *                            &
    1407                                             ( sk(k,j,i)   - sk(k,j,i-1)    )  &
    1408                                         -      (  5.0_wp * ibit2 * adv_sca_5  &
    1409                                             +              ibit1 * adv_sca_3  &
    1410                                                ) *                            &
    1411                                             ( sk(k,j,i+1) - sk(k,j,i-2)    )  &
    1412                                         +      (           ibit2 * adv_sca_5  &
    1413                                                ) *                            &
    1414                                             ( sk(k,j,i+2) - sk(k,j,i-3)    )  &
     1407             swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                    &
     1408                                               ( 10.0_wp * ibit2 * adv_sca_5   &
     1409                                            +     3.0_wp * ibit1 * adv_sca_3   &
     1410                                            +              ibit0 * adv_sca_1   &
     1411                                               ) *                             &
     1412                                            ( sk(k,j,i)   - sk(k,j,i-1)    )   &
     1413                                        -      (  5.0_wp * ibit2 * adv_sca_5   &
     1414                                            +              ibit1 * adv_sca_3   &
     1415                                               ) *                             &
     1416                                            ( sk(k,j,i+1) - sk(k,j,i-2)    )   &
     1417                                        +      (           ibit2 * adv_sca_5   &
     1418                                               ) *                             &
     1419                                            ( sk(k,j,i+2) - sk(k,j,i-3)    )   &
    14151420                                                          )
    14161421
     
    14191424          DO  k = nzb_max_l+1, nzt
    14201425
    1421              u_comp                 = u(k,j,i) - u_gtrans + u_stokes_zu(k)
    1422              swap_flux_x_local(k,j,tn) = u_comp * (                           &
    1423                                       37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) ) &
    1424                                     -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
    1425                                     +           ( sk(k,j,i+2) + sk(k,j,i-3) ) &
     1426             u_comp                    = u(k,j,i) - u_gtrans + u_stokes_zu(k)
     1427             swap_flux_x_local(k,j,tn) = u_comp * (                            &
     1428                                      37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) )  &
     1429                                    -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) )  &
     1430                                    +           ( sk(k,j,i+2) + sk(k,j,i-3) )  &
    14261431                                                  ) * adv_sca_5
    14271432
    1428              swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                   &
    1429                                       10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) ) &
    1430                                     -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
    1431                                     +           ( sk(k,j,i+2) - sk(k,j,i-3) ) &
     1433             swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                    &
     1434                                      10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) )  &
     1435                                    -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) )  &
     1436                                    +           ( sk(k,j,i+2) - sk(k,j,i-3) )  &
    14321437                                                          ) * adv_sca_5
    14331438
    14341439          ENDDO
    1435            
     1440
    14361441       ENDIF
    1437 !       
    1438 !--    Now compute the fluxes and tendency terms for the horizontal and
    1439 !--    vertical parts up to the top of the highest topography.
     1442!
     1443!--    Now compute the fluxes for the horizontal termns up to the highest
     1444!--    topography.
    14401445       DO  k = nzb+1, nzb_max_l
    1441 !
    1442 !--       Note: It is faster to conduct all multiplications explicitly, e.g.
    1443 !--       * adv_sca_5 ... than to determine a factor and multiplicate the
    1444 !--       flux at the end.
    14451446
    14461447          ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
     
    14491450
    14501451          u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
    1451           flux_r(k) = u_comp * (                                              &
    1452                      ( 37.0_wp * ibit2 * adv_sca_5                            &
    1453                   +     7.0_wp * ibit1 * adv_sca_3                            &
    1454                   +              ibit0 * adv_sca_1                            &
    1455                      ) *                                                      &
    1456                              ( sk(k,j,i+1) + sk(k,j,i)   )                    &
    1457               -      (  8.0_wp * ibit2 * adv_sca_5                            &
    1458                   +              ibit1 * adv_sca_3                            &
    1459                      ) *                                                      &
    1460                              ( sk(k,j,i+2) + sk(k,j,i-1) )                    &
    1461               +      (           ibit2 * adv_sca_5                            &
    1462                      ) *                                                      &
    1463                              ( sk(k,j,i+3) + sk(k,j,i-2) )                    &
     1452          flux_r(k) = u_comp * (                                               &
     1453                     ( 37.0_wp * ibit2 * adv_sca_5                             &
     1454                  +     7.0_wp * ibit1 * adv_sca_3                             &
     1455                  +              ibit0 * adv_sca_1                             &
     1456                     ) *                                                       &
     1457                             ( sk(k,j,i+1) + sk(k,j,i)   )                     &
     1458              -      (  8.0_wp * ibit2 * adv_sca_5                             &
     1459                  +              ibit1 * adv_sca_3                             &
     1460                     ) *                                                       &
     1461                             ( sk(k,j,i+2) + sk(k,j,i-1) )                     &
     1462              +      (           ibit2 * adv_sca_5                             &
     1463                     ) *                                                       &
     1464                             ( sk(k,j,i+3) + sk(k,j,i-2) )                     &
    14641465                               )
    14651466
    1466           diss_r(k) = -ABS( u_comp ) * (                                      &
    1467                      ( 10.0_wp * ibit2 * adv_sca_5                            &
    1468                   +     3.0_wp * ibit1 * adv_sca_3                            &
    1469                   +              ibit0 * adv_sca_1                            &
    1470                      ) *                                                      &
    1471                              ( sk(k,j,i+1) - sk(k,j,i)  )                     &
    1472               -      (  5.0_wp * ibit2 * adv_sca_5                            &
    1473                   +              ibit1 * adv_sca_3                            &
    1474                      ) *                                                      &
    1475                              ( sk(k,j,i+2) - sk(k,j,i-1) )                    &
    1476               +      (           ibit2 * adv_sca_5                            &
    1477                      ) *                                                      &
    1478                              ( sk(k,j,i+3) - sk(k,j,i-2) )                    &
     1467          diss_r(k) = -ABS( u_comp ) * (                                       &
     1468                     ( 10.0_wp * ibit2 * adv_sca_5                             &
     1469                  +     3.0_wp * ibit1 * adv_sca_3                             &
     1470                  +              ibit0 * adv_sca_1                             &
     1471                     ) *                                                       &
     1472                             ( sk(k,j,i+1) - sk(k,j,i)  )                      &
     1473              -      (  5.0_wp * ibit2 * adv_sca_5                             &
     1474                  +              ibit1 * adv_sca_3                             &
     1475                     ) *                                                       &
     1476                             ( sk(k,j,i+2) - sk(k,j,i-1) )                     &
     1477              +      (           ibit2 * adv_sca_5                             &
     1478                     ) *                                                       &
     1479                             ( sk(k,j,i+3) - sk(k,j,i-2) )                     &
    14791480                                       )
    14801481
     
    15151516       ENDDO
    15161517!
    1517 !--    Now compute the fluxes and tendency terms for the horizontal and
    1518 !--    vertical parts above the top of the highest topography. No degradation
    1519 !--    for the horizontal parts, but for the vertical it is stell needed.
     1518!--    Now compute the fluxes for the horizontal terms above the topography
     1519!--    where no degradation along the horizontal parts is necessary (except
     1520!--    for the non-cyclic lateral boundaries treated by nzb_max_l).
    15201521       DO  k = nzb_max_l+1, nzt
    15211522
     
    15421543       ENDDO
    15431544!
    1544 !--    Now, compute vertical fluxes. Split loop into a part treating the 
     1545!--    Now, compute vertical fluxes. Split loop into a part treating the
    15451546!--    lowest grid points with indirect indexing, a main loop without
    15461547!--    indirect indexing, and a loop for the uppermost grip points with
    15471548!--    indirect indexing. This allows better vectorization for the main loop.
    1548 !--    First, compute the flux at model surface, which need has to be 
     1549!--    First, compute the flux at model surface, which need has to be
    15491550!--    calculated explicetely for the tendency at
    15501551!--    the first w-level. For topography wall this is done implicitely by
     
    15521553       flux_t(nzb) = 0.0_wp
    15531554       diss_t(nzb) = 0.0_wp
    1554        
     1555
    15551556       DO  k = nzb+1, nzb+1
    15561557          ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     
    15631564          k_pp  = k + 2 * ( 1 - ibit6  )
    15641565          k_mm  = k - 2 * ibit8
    1565 
    1566 
    1567           flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
    1568                      ( 37.0_wp * ibit8 * adv_sca_5                            &
    1569                   +     7.0_wp * ibit7 * adv_sca_3                            &
    1570                   +              ibit6 * adv_sca_1                            &
    1571                      ) *                                                      &
    1572                              ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
    1573               -      (  8.0_wp * ibit8 * adv_sca_5                            &
    1574                   +              ibit7 * adv_sca_3                            &
    1575                      ) *                                                      &
    1576                              ( sk(k_pp,j,i) + sk(k-1,j,i)  )                  &
    1577               +      (           ibit8 * adv_sca_5                            &
    1578                      ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                  &
    1579                                  )
    1580 
    1581           diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
    1582                      ( 10.0_wp * ibit8 * adv_sca_5                            &
    1583                   +     3.0_wp * ibit7 * adv_sca_3                            &
    1584                   +              ibit6 * adv_sca_1                            &
    1585                      ) *                                                      &
    1586                              ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
    1587               -      (  5.0_wp * ibit8 * adv_sca_5                            &
    1588                   +              ibit7 * adv_sca_3                            &
    1589                      ) *                                                      &
    1590                              ( sk(k_pp,j,i)  - sk(k-1,j,i)  )                 &
    1591               +      (           ibit8 * adv_sca_5                            &
    1592                      ) *                                                      &
    1593                              ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    1594                                          )
    1595        ENDDO
    1596        
    1597        DO  k = nzb+2, nzt-2
    1598           ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
    1599           ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
    1600           ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
    1601 
    1602           flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
    1603                      ( 37.0_wp * ibit8 * adv_sca_5                            &
    1604                   +     7.0_wp * ibit7 * adv_sca_3                            &
    1605                   +              ibit6 * adv_sca_1                            &
    1606                      ) *                                                      &
    1607                              ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
    1608               -      (  8.0_wp * ibit8 * adv_sca_5                            &
    1609                   +              ibit7 * adv_sca_3                            &
    1610                      ) *                                                      &
    1611                              ( sk(k+2,j,i) + sk(k-1,j,i)  )                   &
    1612               +      (           ibit8 * adv_sca_5                            &
    1613                      ) *     ( sk(k+3,j,i)+ sk(k-2,j,i) )                     &
    1614                                                  )
    1615 
    1616           diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
    1617                      ( 10.0_wp * ibit8 * adv_sca_5                            &
    1618                   +     3.0_wp * ibit7 * adv_sca_3                            &
    1619                   +              ibit6 * adv_sca_1                            &
    1620                      ) *                                                      &
    1621                              ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
    1622               -      (  5.0_wp * ibit8 * adv_sca_5                            &
    1623                   +              ibit7 * adv_sca_3                            &
    1624                      ) *                                                      &
    1625                              ( sk(k+2,j,i)  - sk(k-1,j,i)  )                  &
    1626               +      (           ibit8 * adv_sca_5                            &
    1627                      ) *                                                      &
    1628                              ( sk(k+3,j,i) - sk(k-2,j,i) )                    &
    1629                                                          )
    1630        ENDDO       
    1631        
    1632        DO  k = nzt-1, nzt-symmetry_flag
    1633           ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
    1634           ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
    1635           ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
    1636 !
    1637 !--       k index has to be modified near bottom and top, else array
    1638 !--       subscripts will be exceeded.
    1639           k_ppp = k + 3 * ibit8
    1640           k_pp  = k + 2 * ( 1 - ibit6  )
    1641           k_mm  = k - 2 * ibit8
    1642 
    16431566
    16441567          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
     
    16711594                                                         )
    16721595       ENDDO
    1673        
     1596
     1597       DO  k = nzb+2, nzt-2
     1598          ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     1599          ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     1600          ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
     1601
     1602          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
     1603                     ( 37.0_wp * ibit8 * adv_sca_5                            &
     1604                  +     7.0_wp * ibit7 * adv_sca_3                            &
     1605                  +              ibit6 * adv_sca_1                            &
     1606                     ) *                                                      &
     1607                             ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
     1608              -      (  8.0_wp * ibit8 * adv_sca_5                            &
     1609                  +              ibit7 * adv_sca_3                            &
     1610                     ) *                                                      &
     1611                             ( sk(k+2,j,i) + sk(k-1,j,i)  )                   &
     1612              +      (           ibit8 * adv_sca_5                            &
     1613                     ) *     ( sk(k+3,j,i)+ sk(k-2,j,i) )                     &
     1614                                                 )
     1615
     1616          diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
     1617                     ( 10.0_wp * ibit8 * adv_sca_5                            &
     1618                  +     3.0_wp * ibit7 * adv_sca_3                            &
     1619                  +              ibit6 * adv_sca_1                            &
     1620                     ) *                                                      &
     1621                             ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
     1622              -      (  5.0_wp * ibit8 * adv_sca_5                            &
     1623                  +              ibit7 * adv_sca_3                            &
     1624                     ) *                                                      &
     1625                             ( sk(k+2,j,i)  - sk(k-1,j,i)  )                  &
     1626              +      (           ibit8 * adv_sca_5                            &
     1627                     ) *                                                      &
     1628                             ( sk(k+3,j,i) - sk(k-2,j,i) )                    &
     1629                                                         )
     1630       ENDDO
     1631
     1632       DO  k = nzt-1, nzt-symmetry_flag
     1633          ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     1634          ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     1635          ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
     1636!
     1637!--       k index has to be modified near bottom and top, else array
     1638!--       subscripts will be exceeded.
     1639          k_ppp = k + 3 * ibit8
     1640          k_pp  = k + 2 * ( 1 - ibit6  )
     1641          k_mm  = k - 2 * ibit8
     1642
     1643
     1644          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
     1645                     ( 37.0_wp * ibit8 * adv_sca_5                            &
     1646                  +     7.0_wp * ibit7 * adv_sca_3                            &
     1647                  +              ibit6 * adv_sca_1                            &
     1648                     ) *                                                      &
     1649                             ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
     1650              -      (  8.0_wp * ibit8 * adv_sca_5                            &
     1651                  +              ibit7 * adv_sca_3                            &
     1652                     ) *                                                      &
     1653                             ( sk(k_pp,j,i) + sk(k-1,j,i)  )                  &
     1654              +      (           ibit8 * adv_sca_5                            &
     1655                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                  &
     1656                                                 )
     1657
     1658          diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
     1659                     ( 10.0_wp * ibit8 * adv_sca_5                            &
     1660                  +     3.0_wp * ibit7 * adv_sca_3                            &
     1661                  +              ibit6 * adv_sca_1                            &
     1662                     ) *                                                      &
     1663                             ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
     1664              -      (  5.0_wp * ibit8 * adv_sca_5                            &
     1665                  +              ibit7 * adv_sca_3                            &
     1666                     ) *                                                      &
     1667                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )                 &
     1668              +      (           ibit8 * adv_sca_5                            &
     1669                     ) *                                                      &
     1670                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
     1671                                                         )
     1672       ENDDO
     1673
    16741674!
    16751675!--    Set resolved/turbulent flux at model top to zero (w-level). In case that
    16761676!--    a symmetric behavior between bottom and top shall be guaranteed (closed
    1677 !--    channel flow), the flux at nzt is also set to zero. 
     1677!--    channel flow), the flux at nzt is also set to zero.
    16781678       IF ( symmetry_flag == 1 ) THEN
    16791679          flux_t(nzt) = 0.0_wp
     
    16821682       flux_t(nzt+1) = 0.0_wp
    16831683       diss_t(nzt+1) = 0.0_wp
    1684        
    1685        
     1684
     1685
    16861686       IF ( limiter )  THEN
    16871687!
     
    16951695!
    16961696!--          In flux limitation the total flux will be corrected. For the sake
    1697 !--          of cleariness the higher-order advective and disspative fluxes 
    1698 !--          will be merged onto flux_t. 
     1697!--          of cleariness the higher-order advective and disspative fluxes
     1698!--          will be merged onto flux_t.
    16991699             flux_t(k) = flux_t(k) + diss_t(k)
    17001700             diss_t(k) = 0.0_wp
    17011701          ENDDO
    17021702!
    1703 !--       Flux limitation of vertical fluxes according to Skamarock (2006). 
     1703!--       Flux limitation of vertical fluxes according to Skamarock (2006).
    17041704!--       Please note, as flux limitation implies linear dependencies of fluxes,
    17051705!--       flux limitation is only made for the vertical advection term. Limitation
    1706 !--       of the horizontal terms cannot be parallelized. 
     1706!--       of the horizontal terms cannot be parallelized.
    17071707!--       Due to the linear dependency, the following loop will not be vectorized.
    17081708!--       Further, note that the flux limiter is only applied within the urban
    1709 !--       layer, i.e up to the topography top. 
     1709!--       layer, i.e up to the topography top.
    17101710          DO  k = nzb+1, nzb_max_l
    17111711!
    17121712!--          Compute one-dimensional divergence along the vertical direction,
    1713 !--          which is used to correct the advection discretization. This is 
     1713!--          which is used to correct the advection discretization. This is
    17141714!--          necessary as in one-dimensional space the advection velocity
    1715 !--          should actually be constant. 
     1715!--          should actually be constant.
    17161716             div = ( w(k,j,i)   * rho_air_zw(k)                                &
    1717                    - w(k-1,j,i) * rho_air_zw(k-1)                              &     
     1717                   - w(k-1,j,i) * rho_air_zw(k-1)                              &
    17181718                   ) * drho_air(k) * ddzw(k)
    17191719!
    1720 !--          Compute monotone solution of the advection equation from 
     1720!--          Compute monotone solution of the advection equation from
    17211721!--          1st-order fluxes. Please note, the advection equation is corrected
    17221722!--          by the divergence term (in 1D the advective flow should be divergence
    1723 !--          free). Moreover, please note, as time-increment the full timestep 
    1724 !--          is used, even though a Runge-Kutta scheme will be used. However, 
    1725 !--          the length of the actual time increment is not important at all 
    1726 !--          since it cancels out later when the fluxes are limited.   
     1723!--          free). Moreover, please note, as time-increment the full timestep
     1724!--          is used, even though a Runge-Kutta scheme will be used. However,
     1725!--          the length of the actual time increment is not important at all
     1726!--          since it cancels out later when the fluxes are limited.
    17271727             mon = sk(k,j,i) + ( - ( flux_t_1st(k) - flux_t_1st(k-1) )         &
    17281728                             * drho_air(k) * ddzw(k)                           &
    17291729                             + div * sk(k,j,i)                                 &
    1730                                ) * dt_3d 
     1730                               ) * dt_3d
    17311731!
    17321732!--          Determine minimum and maximum values along the numerical stencil.
    17331733             k_mmm = MAX( k - 3, nzb + 1 )
    1734              k_ppp = MIN( k + 3, nzt + 1 ) 
     1734             k_ppp = MIN( k + 3, nzt + 1 )
    17351735
    17361736             min_val = MINVAL( sk(k_mmm:k_ppp,j,i) )
    17371737             max_val = MAXVAL( sk(k_mmm:k_ppp,j,i) )
    17381738!
    1739 !--          Compute difference between high- and low-order fluxes, which may 
     1739!--          Compute difference between high- and low-order fluxes, which may
    17401740!--          act as correction fluxes
    17411741             f_corr_t = flux_t(k)   - flux_t_1st(k)
    17421742             f_corr_d = flux_t(k-1) - flux_t_1st(k-1)
    17431743!
    1744 !--          Determine outgoing fluxes, i.e. the part of the fluxes which can 
     1744!--          Determine outgoing fluxes, i.e. the part of the fluxes which can
    17451745!--          decrease the value within the grid box
    17461746             f_corr_t_out = MAX( 0.0_wp, f_corr_t )
    17471747             f_corr_d_out = MIN( 0.0_wp, f_corr_d )
    17481748!
    1749 !--          Determine ingoing fluxes, i.e. the part of the fluxes which can 
    1750 !--          increase the value within the grid box 
     1749!--          Determine ingoing fluxes, i.e. the part of the fluxes which can
     1750!--          increase the value within the grid box
    17511751             f_corr_t_in = MIN( 0.0_wp, f_corr_t)
    17521752             f_corr_d_in = MAX( 0.0_wp, f_corr_d)
     
    17621762!--          Check if outgoing fluxes can lead to undershoots, i.e. values smaller
    17631763!--          than the minimum value within the numerical stencil. If so, limit
    1764 !--          them. 
     1764!--          them.
    17651765             IF ( mon - min_val < - div_out  .AND.  ABS( div_out ) > 0.0_wp )  &
    17661766             THEN
     
    17721772!--          Check if ingoing fluxes can lead to overshoots, i.e. values larger
    17731773!--          than the maximum value within the numerical stencil. If so, limit
    1774 !--          them. 
     1774!--          them.
    17751775             IF ( mon - max_val > - div_in  .AND.  ABS( div_in ) > 0.0_wp )    &
    17761776             THEN
     
    17801780             ENDIF
    17811781!
    1782 !--          Finally add the limited fluxes to the original ones. If no 
    1783 !--          flux limitation was done, the fluxes equal the original ones. 
     1782!--          Finally add the limited fluxes to the original ones. If no
     1783!--          flux limitation was done, the fluxes equal the original ones.
    17841784             flux_t(k)   = flux_t_1st(k)   + f_corr_t_out + f_corr_t_in
    17851785             flux_t(k-1) = flux_t_1st(k-1) + f_corr_d_out + f_corr_d_in
    17861786          ENDDO
    17871787       ENDIF
    1788 
     1788!
     1789!--    Now compute the tendency term including divergence correction.
    17891790       DO  k = nzb+1, nzb_max_l
    17901791
    17911792          flux_d    = flux_t(k-1)
    17921793          diss_d    = diss_t(k-1)
    1793          
     1794
    17941795          ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
    17951796          ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp )
    17961797          ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp )
    1797          
     1798
    17981799          ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp )
    17991800          ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp )
    18001801          ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp )
    1801          
     1802
    18021803          ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
    18031804          ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     
    18071808!--       correction is needed to overcome numerical instabilities introduced
    18081809!--       by a not sufficient reduction of divergences near topography.
    1809           div         =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )            &
    1810                           - u(k,j,i)   * (                                    &
    1811                         REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )      &
    1812                       + REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp )      &
    1813                       + REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp )      &
    1814                                          )                                    &
    1815                           ) * ddx                                             &
    1816                         + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )            &
    1817                           - v(k,j,i)   * (                                    &
    1818                         REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp )      &
    1819                       + REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp )      &
    1820                       + REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp )      &
    1821                                          )                                    &
    1822                           ) * ddy                                             &
    1823                         + ( w(k,j,i) * rho_air_zw(k) *                        &
    1824                                          ( ibit6 + ibit7 + ibit8 )            &
    1825                           - w(k-1,j,i) * rho_air_zw(k-1) *                    &
    1826                                          (                                    &
    1827                         REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp )      &
    1828                       + REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp )      &
    1829                       + REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp )      &
    1830                                          )                                    &     
     1810          div         =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
     1811                          - u(k,j,i)   * (                                     &
     1812                        REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )       &
     1813                      + REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp )       &
     1814                      + REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp )       &
     1815                                         )                                     &
     1816                          ) * ddx                                              &
     1817                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
     1818                          - v(k,j,i)   * (                                     &
     1819                        REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp )       &
     1820                      + REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp )       &
     1821                      + REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp )       &
     1822                                         )                                     &
     1823                          ) * ddy                                              &
     1824                        + ( w(k,j,i) * rho_air_zw(k) *                         &
     1825                                         ( ibit6 + ibit7 + ibit8 )             &
     1826                          - w(k-1,j,i) * rho_air_zw(k-1) *                     &
     1827                                         (                                     &
     1828                        REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp )       &
     1829                      + REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp )       &
     1830                      + REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp )       &
     1831                                         )                                     &
    18311832                          ) * drho_air(k) * ddzw(k)
    18321833
    1833           tend(k,j,i) = tend(k,j,i) - (                                       &
    1834                         ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
    1835                           swap_diss_x_local(k,j,tn)            ) * ddx        &
    1836                       + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
    1837                           swap_diss_y_local(k,tn)              ) * ddy        &
    1838                       + ( ( flux_t(k) + diss_t(k) ) -                         &
    1839                           ( flux_d    + diss_d    )                           &
    1840                                                     ) * drho_air(k) * ddzw(k) &
     1834          tend(k,j,i) = tend(k,j,i) - (                                        &
     1835                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) -  &
     1836                          swap_diss_x_local(k,j,tn)            ) * ddx         &
     1837                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   -  &
     1838                          swap_diss_y_local(k,tn)              ) * ddy         &
     1839                      + ( ( flux_t(k) + diss_t(k) ) -                          &
     1840                          ( flux_d    + diss_d    )                            &
     1841                                                    ) * drho_air(k) * ddzw(k)  &
    18411842                                      ) + sk(k,j,i) * div
    18421843
     
    18481849
    18491850       ENDDO
    1850        
     1851
    18511852       DO  k = nzb_max_l+1, nzt
    18521853
     
    18571858!--       correction is needed to overcome numerical instabilities introduced
    18581859!--       by a not sufficient reduction of divergences near topography.
    1859           div         =   ( u(k,j,i+1) - u(k,j,i) ) * ddx                     &
    1860                         + ( v(k,j+1,i) - v(k,j,i) ) * ddy                     &
    1861                         + ( w(k,j,i) * rho_air_zw(k)                          &
    1862                           - w(k-1,j,i) * rho_air_zw(k-1)                      &
     1860          div         =   ( u(k,j,i+1) - u(k,j,i) ) * ddx                      &
     1861                        + ( v(k,j+1,i) - v(k,j,i) ) * ddy                      &
     1862                        + ( w(k,j,i) * rho_air_zw(k)                           &
     1863                          - w(k-1,j,i) * rho_air_zw(k-1)                       &
    18631864                                                  ) * drho_air(k) * ddzw(k)
    18641865
    1865           tend(k,j,i) = tend(k,j,i) - (                                       &
    1866                         ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
    1867                           swap_diss_x_local(k,j,tn)            ) * ddx        &
    1868                       + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
    1869                           swap_diss_y_local(k,tn)              ) * ddy        &
    1870                       + ( ( flux_t(k) + diss_t(k) ) -                         &
    1871                           ( flux_d    + diss_d    )                           &
    1872                                                     ) * drho_air(k) * ddzw(k) &
     1866          tend(k,j,i) = tend(k,j,i) - (                                        &
     1867                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) -  &
     1868                          swap_diss_x_local(k,j,tn)            ) * ddx         &
     1869                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   -  &
     1870                          swap_diss_y_local(k,tn)              ) * ddy         &
     1871                      + ( ( flux_t(k) + diss_t(k) ) -                          &
     1872                          ( flux_d    + diss_d    )                            &
     1873                                                    ) * drho_air(k) * ddzw(k)  &
    18731874                                      ) + sk(k,j,i) * div
    18741875
     
    18951896                    ) * weight_substep(intermediate_timestep_count)
    18961897             ENDDO
    1897            
     1898
    18981899          CASE ( 'sa' )
    18991900
     
    19061907                    ) * weight_substep(intermediate_timestep_count)
    19071908             ENDDO
    1908            
     1909
    19091910          CASE ( 'q' )
    19101911
     
    19881989          !kk Has to be implemented for kpp chemistry
    19891990
    1990 
    19911991         END SELECT
    19921992
     
    20112011       INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
    20122012       INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
    2013        INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 
     2013       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
    20142014       INTEGER(iwp) ::  tn        !< number of OpenMP thread
    2015        
     2015
    20162016       REAL(wp)    ::  ibit0   !< flag indicating 1st-order scheme along x-direction
    20172017       REAL(wp)    ::  ibit1   !< flag indicating 3rd-order scheme along x-direction
     
    20292029       REAL(wp)    ::  gv       !< Galilei-transformation velocity along y
    20302030       REAL(wp)    ::  u_comp_l !< advection velocity along x at leftmost grid point on subdomain
    2031        
     2031
    20322032       REAL(wp), DIMENSION(nzb:nzt+1) ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
    20332033       REAL(wp), DIMENSION(nzb:nzt+1) ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
     
    20402040       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_comp !< advection velocity along z
    20412041!
    2042 !--    Used local modified copy of nzb_max (used to degrade order of 
     2042!--    Used local modified copy of nzb_max (used to degrade order of
    20432043!--    discretization) at non-cyclic boundaries. Modify only at relevant points
    2044 !--    instead of the entire subdomain. This should lead to better 
     2044!--    instead of the entire subdomain. This should lead to better
    20452045!--    load balance between boundary and non-boundary PEs.
    20462046       IF( ( bc_dirichlet_l  .OR.  bc_radiation_l )  .AND.  i <= nxl + 2  .OR. &
     
    20522052          nzb_max_l = nzb_max
    20532053       END IF
    2054        
     2054
    20552055       gu = 2.0_wp * u_gtrans
    20562056       gv = 2.0_wp * v_gtrans
     
    20582058!--    Compute southside fluxes for the respective boundary of PE
    20592059       IF ( j == nys  )  THEN
    2060        
     2060
    20612061          DO  k = nzb+1, nzb_max_l
    20622062
     
    21042104                           37.0_wp * ( u(k,j,i)   + u(k,j-1,i)   )             &
    21052105                         -  8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) )               &
    2106                          +           ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 
     2106                         +           ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
    21072107             diss_s_u(k,tn) = - ABS(v_comp(k)) * (                             &
    21082108                           10.0_wp * ( u(k,j,i)   - u(k,j-1,i)   )             &
     
    21112111
    21122112          ENDDO
    2113          
     2113
    21142114       ENDIF
    21152115!
    21162116!--    Compute leftside fluxes for the respective boundary of PE
    21172117       IF ( i == i_omp  .OR.  i == nxlu )  THEN
    2118        
     2118
    21192119          DO  k = nzb+1, nzb_max_l
    21202120
     
    21372137                              ) *                                              &
    21382138                                          ( u(k,j,i+2) + u(k,j,i-3) )          &
    2139                                            )                                 
    2140                                                                              
     2139                                           )
     2140
    21412141             diss_l_u(k,j,tn) = - ABS( u_comp_l ) * (                          &
    21422142                              ( 10.0_wp * ibit2 * adv_mom_5                    &
     
    21692169
    21702170          ENDDO
    2171          
     2171
    21722172       ENDIF
    21732173!
     
    21882188                                    ( u(k,j,i+1) + u(k,j,i)   )                &
    21892189              -      (  8.0_wp * ibit2 * adv_mom_5                             &
    2190                   +              ibit1 * adv_mom_3                             & 
     2190                  +              ibit1 * adv_mom_3                             &
    21912191                     ) *                                                       &
    21922192                                    ( u(k,j,i+2) + u(k,j,i-1) )                &
     
    21942194                     ) *                                                       &
    21952195                                    ( u(k,j,i+3) + u(k,j,i-2) )                &
    2196                                            )                                 
    2197                                                                              
     2196                                           )
     2197
    21982198          diss_r(k) = - ABS( u_comp(k) - gu ) * (                              &
    21992199                     ( 10.0_wp * ibit2 * adv_mom_5                             &
     
    22292229                     ) *                                                       &
    22302230                                    ( u(k,j+3,i) + u(k,j-2,i) )                &
    2231                                   )                                           
    2232                                                                              
     2231                                  )
     2232
    22332233          diss_n(k) = - ABS ( v_comp(k) ) * (                                  &
    22342234                     ( 10.0_wp * ibit5 * adv_mom_5                             &
     
    22532253                         37.0_wp * ( u(k,j,i+1) + u(k,j,i)   )                 &
    22542254                       -  8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) )                 &
    2255                        +           ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5   
     2255                       +           ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
    22562256          diss_r(k) = - ABS( u_comp(k) - gu ) * (                              &
    22572257                         10.0_wp * ( u(k,j,i+1) - u(k,j,i)   )                 &
    22582258                       -  5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) )                 &
    2259                        +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5   
    2260                                                                                
    2261           v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv                           
     2259                       +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
     2260
     2261          v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
    22622262          flux_n(k) = v_comp(k) * (                                            &
    22632263                         37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                 &
    22642264                       -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                 &
    2265                        +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5   
     2265                       +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
    22662266          diss_n(k) = - ABS( v_comp(k) ) * (                                   &
    22672267                         10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                 &
     
    22712271       ENDDO
    22722272!
    2273 !--    Now, compute vertical fluxes. Split loop into a part treating the 
     2273!--    Now, compute vertical fluxes. Split loop into a part treating the
    22742274!--    lowest grid points with indirect indexing, a main loop without
    22752275!--    indirect indexing, and a loop for the uppermost grip points with
    22762276!--    indirect indexing. This allows better vectorization for the main loop.
    2277 !--    First, compute the flux at model surface, which need has to be 
     2277!--    First, compute the flux at model surface, which need has to be
    22782278!--    calculated explicitly for the tendency at
    22792279!--    the first w-level. For topography wall this is done implicitely by
     
    22822282       diss_t(nzb) = 0.0_wp
    22832283       w_comp(nzb) = 0.0_wp
    2284        
     2284
    22852285       DO  k = nzb+1, nzb+1
    22862286!
     
    23092309                     ) *                                                       &
    23102310                                ( u(k_ppp,j,i) + u(k_mm,j,i) )                 &
    2311                                                   )                           
    2312                                                                                
     2311                                                  )
     2312
    23132313          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                   &
    23142314                     ( 10.0_wp * ibit8 * adv_mom_5                             &
     
    23262326                                                           )
    23272327       ENDDO
    2328        
     2328
    23292329       DO  k = nzb+2, nzt-2
    23302330
     
    23442344                     ) *                                                       &
    23452345                                ( u(k+2,j,i) + u(k-1,j,i)   )                  &
    2346               +      (           ibit8 * adv_mom_5                             & 
     2346              +      (           ibit8 * adv_mom_5                             &
    23472347                     ) *                                                       &
    23482348                                ( u(k+3,j,i) + u(k-2,j,i) )                    &
     
    23642364                                                           )
    23652365       ENDDO
    2366        
     2366
    23672367       DO  k = nzt-1, nzt-symmetry_flag
    23682368!
     
    24082408                                                           )
    24092409       ENDDO
    2410        
     2410
    24112411!
    24122412!--    Set resolved/turbulent flux at model top to zero (w-level). In case that
    24132413!--    a symmetric behavior between bottom and top shall be guaranteed (closed
    2414 !--    channel flow), the flux at nzt is also set to zero. 
     2414!--    channel flow), the flux at nzt is also set to zero.
    24152415       IF ( symmetry_flag == 1 ) THEN
    24162416          flux_t(nzt) = 0.0_wp
     
    24212421       diss_t(nzt+1) = 0.0_wp
    24222422       w_comp(nzt+1) = 0.0_wp
    2423        
     2423
    24242424       DO  k = nzb+1, nzb_max_l
    24252425
    24262426          flux_d    = flux_t(k-1)
    24272427          diss_d    = diss_t(k-1)
    2428          
     2428
    24292429          ibit2 = REAL( IBITS(advc_flags_m(k,j,i),2,1), KIND = wp )
    24302430          ibit1 = REAL( IBITS(advc_flags_m(k,j,i),1,1), KIND = wp )
    24312431          ibit0 = REAL( IBITS(advc_flags_m(k,j,i),0,1), KIND = wp )
    2432          
     2432
    24332433          ibit5 = REAL( IBITS(advc_flags_m(k,j,i),5,1), KIND = wp )
    24342434          ibit4 = REAL( IBITS(advc_flags_m(k,j,i),4,1), KIND = wp )
    24352435          ibit3 = REAL( IBITS(advc_flags_m(k,j,i),3,1), KIND = wp )
    2436          
     2436
    24372437          ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
    24382438          ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
     
    24642464                   + REAL( IBITS(advc_flags_m(k-1,j,i),7,1), KIND = wp )       &
    24652465                   + REAL( IBITS(advc_flags_m(k-1,j,i),8,1), KIND = wp )       &
    2466                                       )                                        & 
     2466                                      )                                        &
    24672467                  ) * drho_air(k) * ddzw(k)                                    &
    2468                 ) * 0.5_wp                                                     
    2469                                                                                
     2468                ) * 0.5_wp
     2469
    24702470          tend(k,j,i) = tend(k,j,i) - (                                        &
    24712471                            ( flux_r(k) + diss_r(k)                            &
     
    25042504                  ) *   weight_substep(intermediate_timestep_count)
    25052505       ENDDO
    2506        
     2506
    25072507       DO  k = nzb_max_l+1, nzt
    25082508
     
    25162516               +  ( v_comp(k) + gv  - ( v(k,j,i)   + v(k,j,i-1) ) ) * ddy      &
    25172517               +  ( w_comp(k)   * rho_air_zw(k)                                &
    2518                  -  w_comp(k-1) * rho_air_zw(k-1)                              & 
     2518                 -  w_comp(k-1) * rho_air_zw(k-1)                              &
    25192519                  ) * drho_air(k) * ddzw(k)                                    &
    25202520                ) * 0.5_wp
     
    25282528                          -   ( flux_d    + diss_d    )                        &
    25292529                                                    ) * drho_air(k) * ddzw(k)  &
    2530                                        ) + div * u(k,j,i)
     2530                                      ) + div * u(k,j,i)
    25312531
    25322532          flux_l_u(k,j,tn) = flux_r(k)
     
    25782578       INTEGER(iwp)  ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
    25792579       INTEGER(iwp)  ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
    2580        INTEGER(iwp)  ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 
     2580       INTEGER(iwp)  ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
    25812581       INTEGER(iwp)  ::  tn        !< number of OpenMP thread
    2582        
     2582
    25832583       REAL(wp)      ::  ibit9    !< flag indicating 1st-order scheme along x-direction
    25842584       REAL(wp)      ::  ibit10   !< flag indicating 3rd-order scheme along x-direction
    25852585       REAL(wp)      ::  ibit11   !< flag indicating 5th-order scheme along x-direction
    2586        REAL(wp)      ::  ibit12   !< flag indicating 1st-order scheme along y-direction 
     2586       REAL(wp)      ::  ibit12   !< flag indicating 1st-order scheme along y-direction
    25872587       REAL(wp)      ::  ibit13   !< flag indicating 3rd-order scheme along y-direction
    25882588       REAL(wp)      ::  ibit14   !< flag indicating 3rd-order scheme along y-direction
     
    25962596       REAL(wp)      ::  gv       !< Galilei-transformation velocity along y
    25972597       REAL(wp)      ::  v_comp_l !< advection velocity along y on leftmost grid point on subdomain
    2598        
     2598
    25992599       REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
    26002600       REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
     
    26072607       REAL(wp), DIMENSION(nzb:nzt+1)  ::  w_comp !< advection velocity along z
    26082608!
    2609 !--    Used local modified copy of nzb_max (used to degrade order of 
     2609!--    Used local modified copy of nzb_max (used to degrade order of
    26102610!--    discretization) at non-cyclic boundaries. Modify only at relevant points
    2611 !--    instead of the entire subdomain. This should lead to better 
     2611!--    instead of the entire subdomain. This should lead to better
    26122612!--    load balance between boundary and non-boundary PEs.
    26132613       IF( ( bc_dirichlet_l  .OR.  bc_radiation_l )  .AND.  i <= nxl + 2  .OR. &
     
    26192619          nzb_max_l = nzb_max
    26202620       END IF
    2621        
     2621
    26222622       gu = 2.0_wp * u_gtrans
    26232623       gv = 2.0_wp * v_gtrans
    26242624
    2625 !       
     2625!
    26262626!--    Compute leftside fluxes for the respective boundary.
    26272627       IF ( i == i_omp )  THEN
     
    26342634
    26352635             u_comp(k)        = u(k,j-1,i) + u(k,j,i) - gu
    2636              flux_l_v(k,j,tn) = u_comp(k) * (                                 &
    2637                               ( 37.0_wp * ibit11 * adv_mom_5                  &
    2638                            +     7.0_wp * ibit10 * adv_mom_3                  &
    2639                            +              ibit9  * adv_mom_1                  &
    2640                               ) *                                             &
    2641                                         ( v(k,j,i)   + v(k,j,i-1) )           &
    2642                        -      (  8.0_wp * ibit11 * adv_mom_5                  &
    2643                            +              ibit10 * adv_mom_3                  &
    2644                               ) *                                             &
    2645                                         ( v(k,j,i+1) + v(k,j,i-2) )           &
    2646                        +      (           ibit11 * adv_mom_5                  &
    2647                               ) *                                             &
    2648                                         ( v(k,j,i+2) + v(k,j,i-3) )           &
     2636             flux_l_v(k,j,tn) = u_comp(k) * (                                  &
     2637                              ( 37.0_wp * ibit11 * adv_mom_5                   &
     2638                           +     7.0_wp * ibit10 * adv_mom_3                   &
     2639                           +              ibit9  * adv_mom_1                   &
     2640                              ) *                                              &
     2641                                        ( v(k,j,i)   + v(k,j,i-1) )            &
     2642                       -      (  8.0_wp * ibit11 * adv_mom_5                   &
     2643                           +              ibit10 * adv_mom_3                   &
     2644                              ) *                                              &
     2645                                        ( v(k,j,i+1) + v(k,j,i-2) )            &
     2646                       +      (           ibit11 * adv_mom_5                   &
     2647                              ) *                                              &
     2648                                        ( v(k,j,i+2) + v(k,j,i-3) )            &
    26492649                                            )
    26502650
    2651              diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * (                        &
    2652                               ( 10.0_wp * ibit11 * adv_mom_5                  &
    2653                            +     3.0_wp * ibit10 * adv_mom_3                  &
    2654                            +              ibit9  * adv_mom_1                  &
    2655                               ) *                                             &
    2656                                         ( v(k,j,i)   - v(k,j,i-1) )           &
    2657                        -      (  5.0_wp * ibit11 * adv_mom_5                  &
    2658                            +              ibit10 * adv_mom_3                  &
    2659                               ) *                                             &
    2660                                         ( v(k,j,i+1) - v(k,j,i-2) )           &
    2661                        +      (           ibit11 * adv_mom_5                  &
    2662                               ) *                                             &
    2663                                         ( v(k,j,i+2) - v(k,j,i-3) )           &
     2651             diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * (                         &
     2652                              ( 10.0_wp * ibit11 * adv_mom_5                   &
     2653                           +     3.0_wp * ibit10 * adv_mom_3                   &
     2654                           +              ibit9  * adv_mom_1                   &
     2655                              ) *                                              &
     2656                                        ( v(k,j,i)   - v(k,j,i-1) )            &
     2657                       -      (  5.0_wp * ibit11 * adv_mom_5                   &
     2658                           +              ibit10 * adv_mom_3                   &
     2659                              ) *                                              &
     2660                                        ( v(k,j,i+1) - v(k,j,i-2) )            &
     2661                       +      (           ibit11 * adv_mom_5                   &
     2662                              ) *                                              &
     2663                                        ( v(k,j,i+2) - v(k,j,i-3) )            &
    26642664                                                     )
    26652665
     
    26692669
    26702670             u_comp(k)        = u(k,j-1,i) + u(k,j,i) - gu
    2671              flux_l_v(k,j,tn) = u_comp(k) * (                                 &
    2672                              37.0_wp * ( v(k,j,i)   + v(k,j,i-1)   )          &
    2673                            -  8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) )            &
     2671             flux_l_v(k,j,tn) = u_comp(k) * (                                  &
     2672                             37.0_wp * ( v(k,j,i)   + v(k,j,i-1)   )           &
     2673                           -  8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) )             &
    26742674                           +           ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
    2675              diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * (                        &
    2676                              10.0_wp * ( v(k,j,i)   - v(k,j,i-1)   )          &
    2677                            -  5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) )            &
     2675             diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * (                         &
     2676                             10.0_wp * ( v(k,j,i)   - v(k,j,i-1)   )           &
     2677                           -  5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) )             &
    26782678                           +           ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
    26792679
    26802680          ENDDO
    2681          
     2681
    26822682       ENDIF
    26832683!
    26842684!--    Compute southside fluxes for the respective boundary.
    26852685       IF ( j == nysv )  THEN
    2686        
     2686
    26872687          DO  k = nzb+1, nzb_max_l
    26882688
     
    26922692
    26932693             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
    2694              flux_s_v(k,tn) = v_comp_l * (                                    &
    2695                             ( 37.0_wp * ibit14 * adv_mom_5                    &
    2696                          +     7.0_wp * ibit13 * adv_mom_3                    &
    2697                          +              ibit12 * adv_mom_1                    &
    2698                             ) *                                               &
    2699                                         ( v(k,j,i)   + v(k,j-1,i) )           &
    2700                      -      (  8.0_wp * ibit14 * adv_mom_5                    &
    2701                          +              ibit13 * adv_mom_3                    &
    2702                             ) *                                               &
    2703                                         ( v(k,j+1,i) + v(k,j-2,i) )           &
    2704                      +      (           ibit14 * adv_mom_5                    &
    2705                             ) *                                               &
    2706                                         ( v(k,j+2,i) + v(k,j-3,i) )           &
     2694             flux_s_v(k,tn) = v_comp_l * (                                     &
     2695                            ( 37.0_wp * ibit14 * adv_mom_5                     &
     2696                         +     7.0_wp * ibit13 * adv_mom_3                     &
     2697                         +              ibit12 * adv_mom_1                     &
     2698                            ) *                                                &
     2699                                        ( v(k,j,i)   + v(k,j-1,i) )            &
     2700                     -      (  8.0_wp * ibit14 * adv_mom_5                     &
     2701                         +              ibit13 * adv_mom_3                     &
     2702                            ) *                                                &
     2703                                        ( v(k,j+1,i) + v(k,j-2,i) )            &
     2704                     +      (           ibit14 * adv_mom_5                     &
     2705                            ) *                                                &
     2706                                        ( v(k,j+2,i) + v(k,j-3,i) )            &
    27072707                                         )
    27082708
    2709              diss_s_v(k,tn) = - ABS( v_comp_l ) * (                           &
    2710                             ( 10.0_wp * ibit14 * adv_mom_5                    &
    2711                          +     3.0_wp * ibit13 * adv_mom_3                    &
    2712                          +              ibit12 * adv_mom_1                    &
    2713                             ) *                                               &
    2714                                         ( v(k,j,i)   - v(k,j-1,i) )           &
    2715                      -      (  5.0_wp * ibit14 * adv_mom_5                    &
    2716                          +              ibit13 * adv_mom_3                    &
    2717                             ) *                                               &
    2718                                         ( v(k,j+1,i) - v(k,j-2,i) )           &
    2719                      +      (           ibit14 * adv_mom_5                    &
    2720                             ) *                                               &
    2721                                         ( v(k,j+2,i) - v(k,j-3,i) )           &
     2709             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                            &
     2710                            ( 10.0_wp * ibit14 * adv_mom_5                     &
     2711                         +     3.0_wp * ibit13 * adv_mom_3                     &
     2712                         +              ibit12 * adv_mom_1                     &
     2713                            ) *                                                &
     2714                                        ( v(k,j,i)   - v(k,j-1,i) )            &
     2715                     -      (  5.0_wp * ibit14 * adv_mom_5                     &
     2716                         +              ibit13 * adv_mom_3                     &
     2717                            ) *                                                &
     2718                                        ( v(k,j+1,i) - v(k,j-2,i) )            &
     2719                     +      (           ibit14 * adv_mom_5                     &
     2720                            ) *                                                &
     2721                                        ( v(k,j+2,i) - v(k,j-3,i) )            &
    27222722                                                  )
    27232723
     
    27272727
    27282728             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
    2729              flux_s_v(k,tn) = v_comp_l * (                                    &
    2730                            37.0_wp * ( v(k,j,i)   + v(k,j-1,i)   )            &
    2731                          -  8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) )              &
     2729             flux_s_v(k,tn) = v_comp_l * (                                     &
     2730                           37.0_wp * ( v(k,j,i)   + v(k,j-1,i)   )             &
     2731                         -  8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) )               &
    27322732                         +           ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
    2733              diss_s_v(k,tn) = - ABS( v_comp_l ) * (                           &
    2734                            10.0_wp * ( v(k,j,i)   - v(k,j-1,i)   )            &
    2735                          -  5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) )              &
     2733             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                            &
     2734                           10.0_wp * ( v(k,j,i)   - v(k,j-1,i)   )             &
     2735                         -  5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) )               &
    27362736                         +           ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
    27372737
    27382738          ENDDO
    2739          
     2739
    27402740       ENDIF
    27412741!
     
    27472747          ibit10 = REAL( IBITS(advc_flags_m(k,j,i),10,1), KIND = wp )
    27482748          ibit9  = REAL( IBITS(advc_flags_m(k,j,i),9,1),  KIND = wp )
    2749  
     2749
    27502750          u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
    2751           flux_r(k) = u_comp(k) * (                                           &
    2752                      ( 37.0_wp * ibit11 * adv_mom_5                           &
    2753                   +     7.0_wp * ibit10 * adv_mom_3                           &
    2754                   +              ibit9  * adv_mom_1                           &
    2755                      ) *                                                      &
    2756                                     ( v(k,j,i+1) + v(k,j,i)   )               &
    2757               -      (  8.0_wp * ibit11 * adv_mom_5                           &
    2758                   +              ibit10 * adv_mom_3                           &
    2759                      ) *                                                      &
    2760                                     ( v(k,j,i+2) + v(k,j,i-1) )               &
    2761               +      (           ibit11 * adv_mom_5                           &
    2762                      ) *                                                      &
    2763                                     ( v(k,j,i+3) + v(k,j,i-2) )               &
     2751          flux_r(k) = u_comp(k) * (                                            &
     2752                     ( 37.0_wp * ibit11 * adv_mom_5                            &
     2753                  +     7.0_wp * ibit10 * adv_mom_3                            &
     2754                  +              ibit9  * adv_mom_1                            &
     2755                     ) *                                                       &
     2756                                    ( v(k,j,i+1) + v(k,j,i)   )                &
     2757              -      (  8.0_wp * ibit11 * adv_mom_5                            &
     2758                  +              ibit10 * adv_mom_3                            &
     2759                     ) *                                                       &
     2760                                    ( v(k,j,i+2) + v(k,j,i-1) )                &
     2761              +      (           ibit11 * adv_mom_5                            &
     2762                     ) *                                                       &
     2763                                    ( v(k,j,i+3) + v(k,j,i-2) )                &
    27642764                                  )
    27652765
    2766           diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    2767                      ( 10.0_wp * ibit11 * adv_mom_5                           &
    2768                   +     3.0_wp * ibit10 * adv_mom_3                           &
    2769                   +              ibit9  * adv_mom_1                           &
    2770                      ) *                                                      &
    2771                                     ( v(k,j,i+1) - v(k,j,i)  )                &
    2772               -      (  5.0_wp * ibit11 * adv_mom_5                           &
    2773                   +              ibit10 * adv_mom_3                           &
    2774                      ) *                                                      &
    2775                                     ( v(k,j,i+2) - v(k,j,i-1) )               &
    2776               +      (           ibit11 * adv_mom_5                           &
    2777                      ) *                                                      &
    2778                                     ( v(k,j,i+3) - v(k,j,i-2) )               &
     2766          diss_r(k) = - ABS( u_comp(k) ) * (                                   &
     2767                     ( 10.0_wp * ibit11 * adv_mom_5                            &
     2768                  +     3.0_wp * ibit10 * adv_mom_3                            &
     2769                  +              ibit9  * adv_mom_1                            &
     2770                     ) *                                                       &
     2771                                    ( v(k,j,i+1) - v(k,j,i)  )                 &
     2772              -      (  5.0_wp * ibit11 * adv_mom_5                            &
     2773                  +              ibit10 * adv_mom_3                            &
     2774                     ) *                                                       &
     2775                                    ( v(k,j,i+2) - v(k,j,i-1) )                &
     2776              +      (           ibit11 * adv_mom_5                            &
     2777                     ) *                                                       &
     2778                                    ( v(k,j,i+3) - v(k,j,i-2) )                &
    27792779                                           )
    27802780
     
    27852785
    27862786          v_comp(k) = v(k,j+1,i) + v(k,j,i)
    2787           flux_n(k) = ( v_comp(k) - gv ) * (                                  &
    2788                      ( 37.0_wp * ibit14 * adv_mom_5                           &
    2789                   +     7.0_wp * ibit13 * adv_mom_3                           &
    2790                   +              ibit12 * adv_mom_1                           &
    2791                      ) *                                                      &
    2792                                     ( v(k,j+1,i) + v(k,j,i)   )               &
    2793               -      (  8.0_wp * ibit14 * adv_mom_5                           &
    2794                   +              ibit13 * adv_mom_3                           &
    2795                      ) *                                                      &
    2796                                     ( v(k,j+2,i) + v(k,j-1,i) )               &
    2797               +      (           ibit14 * adv_mom_5                           &
    2798                      ) *                                                      &
    2799                                     ( v(k,j+3,i) + v(k,j-2,i) )               &
     2787          flux_n(k) = ( v_comp(k) - gv ) * (                                   &
     2788                     ( 37.0_wp * ibit14 * adv_mom_5                            &
     2789                  +     7.0_wp * ibit13 * adv_mom_3                            &
     2790                  +              ibit12 * adv_mom_1                            &
     2791                     ) *                                                       &
     2792                                    ( v(k,j+1,i) + v(k,j,i)   )                &
     2793              -      (  8.0_wp * ibit14 * adv_mom_5                            &
     2794                  +              ibit13 * adv_mom_3                            &
     2795                     ) *                                                       &
     2796                                    ( v(k,j+2,i) + v(k,j-1,i) )                &
     2797              +      (           ibit14 * adv_mom_5                            &
     2798                     ) *                                                       &
     2799                                    ( v(k,j+3,i) + v(k,j-2,i) )                &
    28002800                                           )
    28012801
    2802           diss_n(k) = - ABS( v_comp(k) - gv ) * (                             &
    2803                      ( 10.0_wp * ibit14 * adv_mom_5                           &
    2804                   +     3.0_wp * ibit13 * adv_mom_3                           &
    2805                   +              ibit12 * adv_mom_1                           &
    2806                      ) *                                                      &
    2807                                     ( v(k,j+1,i) - v(k,j,i)   )               &
    2808               -      (  5.0_wp * ibit14 * adv_mom_5                           &
    2809                   +              ibit13 * adv_mom_3                           &
    2810                      ) *                                                      &
    2811                                     ( v(k,j+2,i) - v(k,j-1,i) )               &
    2812               +      (           ibit14 * adv_mom_5                           &
    2813                      ) *                                                      &
    2814                                     ( v(k,j+3,i) - v(k,j-2,i) )               &
     2802          diss_n(k) = - ABS( v_comp(k) - gv ) * (                              &
     2803                     ( 10.0_wp * ibit14 * adv_mom_5                            &
     2804                  +     3.0_wp * ibit13 * adv_mom_3                            &
     2805                  +              ibit12 * adv_mom_1                            &
     2806                     ) *                                                       &
     2807                                    ( v(k,j+1,i) - v(k,j,i)   )                &
     2808              -      (  5.0_wp * ibit14 * adv_mom_5                            &
     2809                  +              ibit13 * adv_mom_3                            &
     2810                     ) *                                                       &
     2811                                    ( v(k,j+2,i) - v(k,j-1,i) )                &
     2812              +      (           ibit14 * adv_mom_5                            &
     2813                     ) *                                                       &
     2814                                    ( v(k,j+3,i) - v(k,j-2,i) )                &
    28152815                                                )
    28162816       ENDDO
     
    28192819
    28202820          u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
    2821           flux_r(k) = u_comp(k) * (                                           &
    2822                       37.0_wp * ( v(k,j,i+1) + v(k,j,i)   )                   &
    2823                     -  8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) )                   &
     2821          flux_r(k) = u_comp(k) * (                                            &
     2822                      37.0_wp * ( v(k,j,i+1) + v(k,j,i)   )                    &
     2823                    -  8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) )                    &
    28242824                    +           ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
    28252825
    2826           diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    2827                       10.0_wp * ( v(k,j,i+1) - v(k,j,i) )                     &
    2828                     -  5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) )                   &
     2826          diss_r(k) = - ABS( u_comp(k) ) * (                                   &
     2827                      10.0_wp * ( v(k,j,i+1) - v(k,j,i) )                      &
     2828                    -  5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) )                    &
    28292829                    +           ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
    28302830
    28312831
    28322832          v_comp(k) = v(k,j+1,i) + v(k,j,i)
    2833           flux_n(k) = ( v_comp(k) - gv ) * (                                  &
    2834                       37.0_wp * ( v(k,j+1,i) + v(k,j,i)   )                   &
    2835                     -  8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) )                   &
     2833          flux_n(k) = ( v_comp(k) - gv ) * (                                   &
     2834                      37.0_wp * ( v(k,j+1,i) + v(k,j,i)   )                    &
     2835                    -  8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) )                    &
    28362836                      +         ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
    28372837
    2838           diss_n(k) = - ABS( v_comp(k) - gv ) * (                             &
    2839                       10.0_wp * ( v(k,j+1,i) - v(k,j,i)   )                   &
    2840                     -  5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) )                   &
     2838          diss_n(k) = - ABS( v_comp(k) - gv ) * (                              &
     2839                      10.0_wp * ( v(k,j+1,i) - v(k,j,i)   )                    &
     2840                    -  5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) )                    &
    28412841                    +           ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
    28422842       ENDDO
    28432843!
    2844 !--    Now, compute vertical fluxes. Split loop into a part treating the 
     2844!--    Now, compute vertical fluxes. Split loop into a part treating the
    28452845!--    lowest grid points with indirect indexing, a main loop without
    28462846!--    indirect indexing, and a loop for the uppermost grip points with
    28472847!--    indirect indexing. This allows better vectorization for the main loop.
    2848 !--    First, compute the flux at model surface, which need has to be 
     2848!--    First, compute the flux at model surface, which need has to be
    28492849!--    calculated explicitly for the tendency at
    28502850!--    the first w-level. For topography wall this is done implicitely by
     
    28532853       diss_t(nzb) = 0.0_wp
    28542854       w_comp(nzb) = 0.0_wp
    2855        
     2855
    28562856       DO  k = nzb+1, nzb+1
    28572857!
     
    28672867
    28682868          w_comp(k) = w(k,j-1,i) + w(k,j,i)
    2869           flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    2870                      ( 37.0_wp * ibit17 * adv_mom_5                           &
    2871                   +     7.0_wp * ibit16 * adv_mom_3                           &
    2872                   +              ibit15 * adv_mom_1                           &
    2873                      ) *                                                      &
    2874                                 ( v(k+1,j,i)   + v(k,j,i)    )                &
    2875               -      (  8.0_wp * ibit17 * adv_mom_5                           &
    2876                   +              ibit16 * adv_mom_3                           &
    2877                      ) *                                                      &
    2878                                 ( v(k_pp,j,i)  + v(k-1,j,i)  )                &
    2879               +      (           ibit17 * adv_mom_5                           &
    2880                      ) *                                                      &
    2881                                 ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
     2869          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                            &
     2870                     ( 37.0_wp * ibit17 * adv_mom_5                            &
     2871                  +     7.0_wp * ibit16 * adv_mom_3                            &
     2872                  +              ibit15 * adv_mom_1                            &
     2873                     ) *                                                       &
     2874                                ( v(k+1,j,i)   + v(k,j,i)    )                 &
     2875              -      (  8.0_wp * ibit17 * adv_mom_5                            &
     2876                  +              ibit16 * adv_mom_3                            &
     2877                     ) *                                                       &
     2878                                ( v(k_pp,j,i)  + v(k-1,j,i)  )                 &
     2879              +      (           ibit17 * adv_mom_5                            &
     2880                     ) *                                                       &
     2881                                ( v(k_ppp,j,i) + v(k_mm,j,i) )                 &
    28822882                                                  )
    28832883
    2884           diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    2885                      ( 10.0_wp * ibit17 * adv_mom_5                           &
    2886                   +     3.0_wp * ibit16 * adv_mom_3                           &
    2887                   +              ibit15 * adv_mom_1                           &
    2888                      ) *                                                      &
    2889                                 ( v(k+1,j,i)   - v(k,j,i)    )                &
    2890               -      (  5.0_wp * ibit17 * adv_mom_5                           &
    2891                   +              ibit16 * adv_mom_3                           &
    2892                      ) *                                                      &
    2893                                 ( v(k_pp,j,i)  - v(k-1,j,i)  )                &
    2894               +      (           ibit17 * adv_mom_5                           &
    2895                      ) *                                                      &
    2896                                 ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
     2884          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                   &
     2885                     ( 10.0_wp * ibit17 * adv_mom_5                            &
     2886                  +     3.0_wp * ibit16 * adv_mom_3                            &
     2887                  +              ibit15 * adv_mom_1                            &
     2888                     ) *                                                       &
     2889                                ( v(k+1,j,i)   - v(k,j,i)    )                 &
     2890              -      (  5.0_wp * ibit17 * adv_mom_5                            &
     2891                  +              ibit16 * adv_mom_3                            &
     2892                     ) *                                                       &
     2893                                ( v(k_pp,j,i)  - v(k-1,j,i)  )                 &
     2894              +      (           ibit17 * adv_mom_5                            &
     2895                     ) *                                                       &
     2896                                ( v(k_ppp,j,i) - v(k_mm,j,i) )                 &
    28972897                                                           )
    28982898       ENDDO
    2899        
     2899
    29002900       DO  k = nzb+2, nzt-2
    29012901
     
    29052905
    29062906          w_comp(k) = w(k,j-1,i) + w(k,j,i)
    2907           flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    2908                      ( 37.0_wp * ibit17 * adv_mom_5                           &
    2909                   +     7.0_wp * ibit16 * adv_mom_3                           &
    2910                   +              ibit15 * adv_mom_1                           &
    2911                      ) *                                                      &
    2912                                 ( v(k+1,j,i) + v(k,j,i)   )                   &
    2913               -      (  8.0_wp * ibit17 * adv_mom_5                           &
    2914                   +              ibit16 * adv_mom_3                           &
    2915                      ) *                                                      &
    2916                                 ( v(k+2,j,i) + v(k-1,j,i) )                   &
    2917               +      (           ibit17 * adv_mom_5                           &
    2918                      ) *                                                      &
    2919                                 ( v(k+3,j,i) + v(k-2,j,i) )                   &
     2907          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                            &
     2908                     ( 37.0_wp * ibit17 * adv_mom_5                            &
     2909                  +     7.0_wp * ibit16 * adv_mom_3                            &
     2910                  +              ibit15 * adv_mom_1                            &
     2911                     ) *                                                       &
     2912                                ( v(k+1,j,i) + v(k,j,i)   )                    &
     2913              -      (  8.0_wp * ibit17 * adv_mom_5                            &
     2914                  +              ibit16 * adv_mom_3                            &
     2915                     ) *                                                       &
     2916                                ( v(k+2,j,i) + v(k-1,j,i) )                    &
     2917              +      (           ibit17 * adv_mom_5                            &
     2918                     ) *                                                       &
     2919                                ( v(k+3,j,i) + v(k-2,j,i) )                    &
    29202920                                                  )
    29212921
    2922           diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    2923                      ( 10.0_wp * ibit17 * adv_mom_5                           &
    2924                   +     3.0_wp * ibit16 * adv_mom_3                           &
    2925                   +              ibit15 * adv_mom_1                           &
    2926                      ) *                                                      &
    2927                                 ( v(k+1,j,i) - v(k,j,i)   )                   &
    2928               -      (  5.0_wp * ibit17 * adv_mom_5                           &
    2929                   +              ibit16 * adv_mom_3                           &
    2930                      ) *                                                      &
     2922          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                   &
     2923                     ( 10.0_wp * ibit17 * adv_mom_5                            &
     2924                  +     3.0_wp * ibit16 * adv_mom_3                            &
     2925                  +              ibit15 * adv_mom_1                            &
     2926                     ) *                                                       &
     2927                                ( v(k+1,j,i) - v(k,j,i)   )                    &
     2928              -      (  5.0_wp * ibit17 * adv_mom_5                            &
     2929                  +              ibit16 * adv_mom_3                            &
     2930                     ) *                                                       &
    29312931                                ( v(k+2,j,i) - v(k-1,j,i) )                    &
    2932               +      (           ibit17 * adv_mom_5                           &
    2933                      ) *                                                      &
    2934                                 ( v(k+3,j,i) - v(k-2,j,i) )                   &
     2932              +      (           ibit17 * adv_mom_5                            &
     2933                     ) *                                                       &
     2934                                ( v(k+3,j,i) - v(k-2,j,i) )                    &
    29352935                                                           )
    29362936       ENDDO
    2937        
     2937
    29382938       DO  k = nzt-1, nzt-symmetry_flag
    29392939!
     
    29492949
    29502950          w_comp(k) = w(k,j-1,i) + w(k,j,i)
    2951           flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    2952                      ( 37.0_wp * ibit17 * adv_mom_5                           &
    2953                   +     7.0_wp * ibit16 * adv_mom_3                           &
    2954                   +              ibit15 * adv_mom_1                           &
    2955                      ) *                                                      &
    2956                                 ( v(k+1,j,i)   + v(k,j,i)    )                &
    2957               -      (  8.0_wp * ibit17 * adv_mom_5                           &
    2958                   +              ibit16 * adv_mom_3                           &
    2959                      ) *                                                      &
    2960                                 ( v(k_pp,j,i)  + v(k-1,j,i)  )                &
    2961               +      (           ibit17 * adv_mom_5                           &
    2962                      ) *                                                      &
    2963                                 ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
     2951          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                            &
     2952                     ( 37.0_wp * ibit17 * adv_mom_5                            &
     2953                  +     7.0_wp * ibit16 * adv_mom_3                            &
     2954                  +              ibit15 * adv_mom_1                            &
     2955                     ) *                                                       &
     2956                                ( v(k+1,j,i)   + v(k,j,i)    )                 &
     2957              -      (  8.0_wp * ibit17 * adv_mom_5                            &
     2958                  +              ibit16 * adv_mom_3                            &
     2959                     ) *                                                       &
     2960                                ( v(k_pp,j,i)  + v(k-1,j,i)  )                 &
     2961              +      (           ibit17 * adv_mom_5                            &
     2962                     ) *                                                       &
     2963                                ( v(k_ppp,j,i) + v(k_mm,j,i) )                 &
    29642964                                                  )
    29652965
    2966           diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    2967                      ( 10.0_wp * ibit17 * adv_mom_5                           &
    2968                   +     3.0_wp * ibit16 * adv_mom_3                           &
    2969                   +              ibit15 * adv_mom_1                           &
    2970                      ) *                                                      &
    2971                                 ( v(k+1,j,i)   - v(k,j,i)    )                &
    2972               -      (  5.0_wp * ibit17 * adv_mom_5                           &
    2973                   +              ibit16 * adv_mom_3                           &
    2974                      ) *                                                      &
    2975                                 ( v(k_pp,j,i)  - v(k-1,j,i)  )                &
    2976               +      (           ibit17 * adv_mom_5                           &
    2977                      ) *                                                      &
    2978                                 ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
     2966          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                   &
     2967                     ( 10.0_wp * ibit17 * adv_mom_5                            &
     2968                  +     3.0_wp * ibit16 * adv_mom_3                            &
     2969                  +              ibit15 * adv_mom_1                            &
     2970                     ) *                                                       &
     2971                                ( v(k+1,j,i)   - v(k,j,i)    )                 &
     2972              -      (  5.0_wp * ibit17 * adv_mom_5                            &
     2973                  +              ibit16 * adv_mom_3                            &
     2974                     ) *                                                       &
     2975                                ( v(k_pp,j,i)  - v(k-1,j,i)  )                 &
     2976              +      (           ibit17 * adv_mom_5                            &
     2977                     ) *                                                       &
     2978                                ( v(k_ppp,j,i) - v(k_mm,j,i) )                 &
    29792979                                                           )
    29802980       ENDDO
    2981        
     2981
    29822982!
    29832983!--    Set resolved/turbulent flux at model top to zero (w-level). In case that
    29842984!--    a symmetric behavior between bottom and top shall be guaranteed (closed
    2985 !--    channel flow), the flux at nzt is also set to zero. 
     2985!--    channel flow), the flux at nzt is also set to zero.
    29862986       IF ( symmetry_flag == 1 ) THEN
    29872987          flux_t(nzt) = 0.0_wp
     
    29922992       diss_t(nzt+1) = 0.0_wp
    29932993       w_comp(nzt+1) = 0.0_wp
    2994        
     2994
    29952995       DO  k = nzb+1, nzb_max_l
    29962996
    29972997          flux_d    = flux_t(k-1)
    29982998          diss_d    = diss_t(k-1)
    2999          
     2999
    30003000          ibit11 = REAL( IBITS(advc_flags_m(k,j,i),11,1), KIND = wp )
    30013001          ibit10 = REAL( IBITS(advc_flags_m(k,j,i),10,1), KIND = wp )
    30023002          ibit9  = REAL( IBITS(advc_flags_m(k,j,i),9,1),  KIND = wp )
    3003          
     3003
    30043004          ibit14 = REAL( IBITS(advc_flags_m(k,j,i),14,1), KIND = wp )
    30053005          ibit13 = REAL( IBITS(advc_flags_m(k,j,i),13,1), KIND = wp )
    30063006          ibit12 = REAL( IBITS(advc_flags_m(k,j,i),12,1), KIND = wp )
    3007          
     3007
    30083008          ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp )
    30093009          ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp )
     
    30133013!--       correction is needed to overcome numerical instabilities introduced
    30143014!--       by a not sufficient reduction of divergences near topography.
    3015           div = ( ( ( u_comp(k)     + gu )                                    &
    3016                                        * ( ibit9 + ibit10 + ibit11 )          &
    3017                   - ( u(k,j-1,i) + u(k,j,i) )                                 &
    3018                                        * (                                    &
    3019                         REAL( IBITS(advc_flags_m(k,j,i-1),9,1),  KIND = wp )  &
    3020                       + REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp )  &
    3021                       + REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp )  &
    3022                                          )                                    &
    3023                   ) * ddx                                                     &
    3024                +  ( v_comp(k)                                                 &
    3025                                        * ( ibit12 + ibit13 + ibit14 )         &
    3026                 - ( v(k,j,i)     + v(k,j-1,i) )                               &
    3027                                        * (                                    &
    3028                         REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp )  &
    3029                       + REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp )  &
    3030                       + REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp )  &
    3031                                          )                                    &
    3032                   ) * ddy                                                     &
    3033                +  ( w_comp(k)   * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )&
    3034                 -   w_comp(k-1) * rho_air_zw(k-1)                             &
    3035                                        * (                                    &
    3036                         REAL( IBITS(advc_flags_m(k-1,j,i),15,1), KIND = wp )  &
    3037                       + REAL( IBITS(advc_flags_m(k-1,j,i),16,1), KIND = wp )  &
    3038                       + REAL( IBITS(advc_flags_m(k-1,j,i),17,1), KIND = wp )  &
    3039                                          )                                    &
    3040                    ) * drho_air(k) * ddzw(k)                                  &
     3015          div = ( ( ( u_comp(k)     + gu )                                     &
     3016                                       * ( ibit9 + ibit10 + ibit11 )           &
     3017                  - ( u(k,j-1,i) + u(k,j,i) )                                  &
     3018                                       * (                                     &
     3019                        REAL( IBITS(advc_flags_m(k,j,i-1),9,1),  KIND = wp )   &
     3020                      + REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp )   &
     3021                      + REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp )   &
     3022                                         )                                     &
     3023                  ) * ddx                                                      &
     3024               +  ( v_comp(k)                                                  &
     3025                                       * ( ibit12 + ibit13 + ibit14 )          &
     3026                - ( v(k,j,i)     + v(k,j-1,i) )                                &
     3027                                       * (                                     &
     3028                        REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp )   &
     3029                      + REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp )   &
     3030                      + REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp )   &
     3031                                         )                                     &
     3032                  ) * ddy                                                      &
     3033               +  ( w_comp(k)   * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 ) &
     3034                -   w_comp(k-1) * rho_air_zw(k-1)                              &
     3035                                       * (                                     &
     3036                        REAL( IBITS(advc_flags_m(k-1,j,i),15,1), KIND = wp )   &
     3037                      + REAL( IBITS(advc_flags_m(k-1,j,i),16,1), KIND = wp )   &
     3038                      + REAL( IBITS(advc_flags_m(k-1,j,i),17,1), KIND = wp )   &
     3039                                         )                                     &
     3040                   ) * drho_air(k) * ddzw(k)                                   &
    30413041                ) * 0.5_wp
    30423042
    3043           tend(k,j,i) = tend(k,j,i) - (                                       &
    3044                          ( flux_r(k) + diss_r(k)                              &
    3045                        -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx      &
    3046                        + ( flux_n(k) + diss_n(k)                              &
    3047                        -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy      &
    3048                        + ( ( flux_t(k) + diss_t(k) )                          &
    3049                        -   ( flux_d    + diss_d    )                          &
    3050                                                    ) * drho_air(k) * ddzw(k)  &
     3043          tend(k,j,i) = tend(k,j,i) - (                                        &
     3044                         ( flux_r(k) + diss_r(k)                               &
     3045                       -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx       &
     3046                       + ( flux_n(k) + diss_n(k)                               &
     3047                       -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy       &
     3048                       + ( ( flux_t(k) + diss_t(k) )                           &
     3049                       -   ( flux_d    + diss_d    )                           &
     3050                                                   ) * drho_air(k) * ddzw(k)   &
    30513051                                      ) + v(k,j,i) * div
    30523052
     
    30783078
    30793079       ENDDO
    3080        
     3080
    30813081       DO  k = nzb_max_l+1, nzt
    30823082
     
    31523152       INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
    31533153       INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
    3154        INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 
     3154       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
    31553155       INTEGER(iwp) ::  tn        !< number of OpenMP thread
    3156        
     3156
    31573157       REAL(wp)    ::  ibit18  !< flag indicating 1st-order scheme along x-direction
    31583158       REAL(wp)    ::  ibit19  !< flag indicating 3rd-order scheme along x-direction
     
    31693169       REAL(wp)    ::  gu      !< Galilei-transformation velocity along x
    31703170       REAL(wp)    ::  gv      !< Galilei-transformation velocity along y
    3171        
     3171
    31723172       REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
    31733173       REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
     
    31803180       REAL(wp), DIMENSION(nzb:nzt+1)  ::  w_comp !< advection velocity along z
    31813181!
    3182 !--    Used local modified copy of nzb_max (used to degrade order of 
     3182!--    Used local modified copy of nzb_max (used to degrade order of
    31833183!--    discretization) at non-cyclic boundaries. Modify only at relevant points
    3184 !--    instead of the entire subdomain. This should lead to better 
     3184!--    instead of the entire subdomain. This should lead to better
    31853185!--    load balance between boundary and non-boundary PEs.
    31863186       IF( ( bc_dirichlet_l  .OR.  bc_radiation_l )  .AND.  i <= nxl + 2  .OR. &
     
    31923192          nzb_max_l = nzb_max
    31933193       END IF
    3194        
     3194
    31953195       gu = 2.0_wp * u_gtrans
    31963196       gv = 2.0_wp * v_gtrans
     
    32053205
    32063206             v_comp(k)      = v(k+1,j,i) + v(k,j,i) - gv
    3207              flux_s_w(k,tn) = v_comp(k) * (                                   &
    3208                             ( 37.0_wp * ibit23 * adv_mom_5                    &
    3209                          +     7.0_wp * ibit22 * adv_mom_3                    &
    3210                          +              ibit21 * adv_mom_1                    &
    3211                             ) *                                               &
    3212                                         ( w(k,j,i)   + w(k,j-1,i) )           &
    3213                      -      (  8.0_wp * ibit23 * adv_mom_5                    &
    3214                          +              ibit22 * adv_mom_3                    &
    3215                             ) *                                               &
    3216                                         ( w(k,j+1,i) + w(k,j-2,i) )           &
    3217                      +      (           ibit23 * adv_mom_5                    &
    3218                             ) *                                               &
    3219                                         ( w(k,j+2,i) + w(k,j-3,i) )           &
     3207             flux_s_w(k,tn) = v_comp(k) * (                                    &
     3208                            ( 37.0_wp * ibit23 * adv_mom_5                     &
     3209                         +     7.0_wp * ibit22 * adv_mom_3                     &
     3210                         +              ibit21 * adv_mom_1                     &
     3211                            ) *                                                &
     3212                                        ( w(k,j,i)   + w(k,j-1,i) )            &
     3213                     -      (  8.0_wp * ibit23 * adv_mom_5                     &
     3214                         +              ibit22 * adv_mom_3                     &
     3215                            ) *                                                &
     3216                                        ( w(k,j+1,i) + w(k,j-2,i) )            &
     3217                     +      (           ibit23 * adv_mom_5                     &
     3218                            ) *                                                &
     3219                                        ( w(k,j+2,i) + w(k,j-3,i) )            &
    32203220                                          )
    32213221
    3222              diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                          &
    3223                             ( 10.0_wp * ibit23 * adv_mom_5                    &
    3224                          +     3.0_wp * ibit22 * adv_mom_3                    &
    3225                          +              ibit21 * adv_mom_1                    &
    3226                             ) *                                               &
    3227                                         ( w(k,j,i)   - w(k,j-1,i) )           &
    3228                      -      (  5.0_wp * ibit23 * adv_mom_5                    &
    3229                          +              ibit22 * adv_mom_3                    &
    3230                             ) *                                               &
    3231                                         ( w(k,j+1,i) - w(k,j-2,i) )           &
    3232                      +      (           ibit23 * adv_mom_5                    &
    3233                             ) *                                               &
    3234                                         ( w(k,j+2,i) - w(k,j-3,i) )           &
     3222             diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                           &
     3223                            ( 10.0_wp * ibit23 * adv_mom_5                     &
     3224                         +     3.0_wp * ibit22 * adv_mom_3                     &
     3225                         +              ibit21 * adv_mom_1                     &
     3226                            ) *                                                &
     3227                                        ( w(k,j,i)   - w(k,j-1,i) )            &
     3228                     -      (  5.0_wp * ibit23 * adv_mom_5                     &
     3229                         +              ibit22 * adv_mom_3                     &
     3230                            ) *                                                &
     3231                                        ( w(k,j+1,i) - w(k,j-2,i) )            &
     3232                     +      (           ibit23 * adv_mom_5                     &
     3233                            ) *                                                &
     3234                                        ( w(k,j+2,i) - w(k,j-3,i) )            &
    32353235                                                   )
    32363236
     
    32403240
    32413241             v_comp(k)      = v(k+1,j,i) + v(k,j,i) - gv
    3242              flux_s_w(k,tn) = v_comp(k) * (                                   &
    3243                            37.0_wp * ( w(k,j,i)   + w(k,j-1,i) )              &
    3244                          -  8.0_wp * ( w(k,j+1,i) + w(k,j-2,i) )              &
     3242             flux_s_w(k,tn) = v_comp(k) * (                                    &
     3243                           37.0_wp * ( w(k,j,i)   + w(k,j-1,i) )               &
     3244                         -  8.0_wp * ( w(k,j+1,i) + w(k,j-2,i) )               &
    32453245                         +           ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
    3246              diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                          &
    3247                            10.0_wp * ( w(k,j,i)   - w(k,j-1,i) )              &
    3248                          -  5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) )              &
     3246             diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                           &
     3247                           10.0_wp * ( w(k,j,i)   - w(k,j-1,i) )               &
     3248                         -  5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) )               &
    32493249                         +           ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
    32503250
     
    32633263
    32643264             u_comp(k)        = u(k+1,j,i) + u(k,j,i) - gu
    3265              flux_l_w(k,j,tn) = u_comp(k) * (                                 &
    3266                              ( 37.0_wp * ibit20 * adv_mom_5                   &
    3267                           +     7.0_wp * ibit19 * adv_mom_3                   &
    3268                           +              ibit18 * adv_mom_1                   &
    3269                              ) *                                              &
    3270                                         ( w(k,j,i)   + w(k,j,i-1) )           &
    3271                       -      (  8.0_wp * ibit20 * adv_mom_5                   &
    3272                           +              ibit19 * adv_mom_3                   &
    3273                              ) *                                              &
    3274                                         ( w(k,j,i+1) + w(k,j,i-2) )           &
    3275                       +      (           ibit20 * adv_mom_5                   &
    3276                              ) *                                              &
    3277                                         ( w(k,j,i+2) + w(k,j,i-3) )           &
     3265             flux_l_w(k,j,tn) = u_comp(k) * (                                  &
     3266                             ( 37.0_wp * ibit20 * adv_mom_5                    &
     3267                          +     7.0_wp * ibit19 * adv_mom_3                    &
     3268                          +              ibit18 * adv_mom_1                    &
     3269                             ) *                                               &
     3270                                        ( w(k,j,i)   + w(k,j,i-1) )            &
     3271                      -      (  8.0_wp * ibit20 * adv_mom_5                    &
     3272                          +              ibit19 * adv_mom_3                    &
     3273                             ) *                                               &
     3274                                        ( w(k,j,i+1) + w(k,j,i-2) )            &
     3275                      +      (           ibit20 * adv_mom_5                    &
     3276                             ) *                                               &
     3277                                        ( w(k,j,i+2) + w(k,j,i-3) )            &
    32783278                                            )
    32793279
    3280              diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * (                        &
    3281                              ( 10.0_wp * ibit20 * adv_mom_5                   &
    3282                           +     3.0_wp * ibit19 * adv_mom_3                   &
    3283                           +              ibit18 * adv_mom_1                   &
    3284                              ) *                                              &
    3285                                         ( w(k,j,i)   - w(k,j,i-1) )           &
    3286                       -      (  5.0_wp * ibit20 * adv_mom_5                   &
    3287                           +              ibit19 * adv_mom_3                   &
    3288                              ) *                                              &
    3289                                         ( w(k,j,i+1) - w(k,j,i-2) )           &
    3290                       +      (           ibit20 * adv_mom_5                   &
    3291                              ) *                                              &
    3292                                         ( w(k,j,i+2) - w(k,j,i-3) )           &
     3280             diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * (                         &
     3281                             ( 10.0_wp * ibit20 * adv_mom_5                    &
     3282                          +     3.0_wp * ibit19 * adv_mom_3                    &
     3283                          +              ibit18 * adv_mom_1                    &
     3284                             ) *                                               &
     3285                                        ( w(k,j,i)   - w(k,j,i-1) )            &
     3286                      -      (  5.0_wp * ibit20 * adv_mom_5                    &
     3287                          +              ibit19 * adv_mom_3                    &
     3288                             ) *                                               &
     3289                                        ( w(k,j,i+1) - w(k,j,i-2) )            &
     3290                      +      (           ibit20 * adv_mom_5                    &
     3291                             ) *                                               &
     3292                                        ( w(k,j,i+2) - w(k,j,i-3) )            &
    32933293                                                     )
    32943294
     
    32983298
    32993299             u_comp(k)        = u(k+1,j,i) + u(k,j,i) - gu
    3300              flux_l_w(k,j,tn) = u_comp(k) * (                                 &
    3301                             37.0_wp * ( w(k,j,i)   + w(k,j,i-1) )             &
    3302                           -  8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) )             &
     3300             flux_l_w(k,j,tn) = u_comp(k) * (                                  &
     3301                            37.0_wp * ( w(k,j,i)   + w(k,j,i-1) )              &
     3302                          -  8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) )              &
    33033303                          +           ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
    3304              diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * (                        &
    3305                             10.0_wp * ( w(k,j,i)   - w(k,j,i-1) )             &
    3306                           -  5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) )             &
    3307                           +           ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 
     3304             diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * (                         &
     3305                            10.0_wp * ( w(k,j,i)   - w(k,j,i-1) )              &
     3306                          -  5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) )              &
     3307                          +           ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
    33083308
    33093309          ENDDO
     
    33203320
    33213321          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
    3322           flux_r(k) = u_comp(k) * (                                           &
    3323                      ( 37.0_wp * ibit20 * adv_mom_5                           &
    3324                   +     7.0_wp * ibit19 * adv_mom_3                           &
    3325                   +              ibit18 * adv_mom_1                           &
    3326                      ) *                                                      &
    3327                                     ( w(k,j,i+1) + w(k,j,i)   )               &
    3328               -      (  8.0_wp * ibit20 * adv_mom_5                           &
    3329                   +              ibit19 * adv_mom_3                           &
    3330                      ) *                                                      &
    3331                                     ( w(k,j,i+2) + w(k,j,i-1) )               &
    3332               +      (           ibit20 * adv_mom_5                           &
    3333                      ) *                                                      &
    3334                                     ( w(k,j,i+3) + w(k,j,i-2) )               &
     3322          flux_r(k) = u_comp(k) * (                                            &
     3323                     ( 37.0_wp * ibit20 * adv_mom_5                            &
     3324                  +     7.0_wp * ibit19 * adv_mom_3                            &
     3325                  +              ibit18 * adv_mom_1                            &
     3326                     ) *                                                       &
     3327                                    ( w(k,j,i+1) + w(k,j,i)   )                &
     3328              -      (  8.0_wp * ibit20 * adv_mom_5                            &
     3329                  +              ibit19 * adv_mom_3                            &
     3330                     ) *                                                       &
     3331                                    ( w(k,j,i+2) + w(k,j,i-1) )                &
     3332              +      (           ibit20 * adv_mom_5                            &
     3333                     ) *                                                       &
     3334                                    ( w(k,j,i+3) + w(k,j,i-2) )                &
    33353335                                  )
    33363336
    3337           diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    3338                      ( 10.0_wp * ibit20 * adv_mom_5                           &
    3339                   +     3.0_wp * ibit19 * adv_mom_3                           &
    3340                   +              ibit18 * adv_mom_1                           &
    3341                      ) *                                                      &
    3342                                     ( w(k,j,i+1) - w(k,j,i)   )               &
    3343               -      (  5.0_wp * ibit20 * adv_mom_5                           &
    3344                   +              ibit19 * adv_mom_3                           &
    3345                      ) *                                                      &
    3346                                     ( w(k,j,i+2) - w(k,j,i-1) )               &
    3347               +      (           ibit20 * adv_mom_5                           &
    3348                      ) *                                                      &
    3349                                     ( w(k,j,i+3) - w(k,j,i-2) )               &
     3337          diss_r(k) = - ABS( u_comp(k) ) * (                                   &
     3338                     ( 10.0_wp * ibit20 * adv_mom_5                            &
     3339                  +     3.0_wp * ibit19 * adv_mom_3                            &
     3340                  +              ibit18 * adv_mom_1                            &
     3341                     ) *                                                       &
     3342                                    ( w(k,j,i+1) - w(k,j,i)   )                &
     3343              -      (  5.0_wp * ibit20 * adv_mom_5                            &
     3344                  +              ibit19 * adv_mom_3                            &
     3345                     ) *                                                       &
     3346                                    ( w(k,j,i+2) - w(k,j,i-1) )                &
     3347              +      (           ibit20 * adv_mom_5                            &
     3348                     ) *                                                       &
     3349                                    ( w(k,j,i+3) - w(k,j,i-2) )                &
    33503350                                           )
    33513351
     
    33553355
    33563356          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
    3357           flux_n(k) = v_comp(k) * (                                           &
    3358                      ( 37.0_wp * ibit23 * adv_mom_5                           &
    3359                   +     7.0_wp * ibit22 * adv_mom_3                           &
    3360                   +              ibit21 * adv_mom_1                           &
    3361                      ) *                                                      &
    3362                                     ( w(k,j+1,i) + w(k,j,i)   )               &
    3363               -      (  8.0_wp * ibit23 * adv_mom_5                           &
    3364                   +              ibit22 * adv_mom_3                           &
    3365                      ) *                                                      &
    3366                                     ( w(k,j+2,i) + w(k,j-1,i) )               &
    3367               +      (           ibit23 * adv_mom_5                           &
    3368                      ) *                                                      &
    3369                                     ( w(k,j+3,i) + w(k,j-2,i) )               &
     3357          flux_n(k) = v_comp(k) * (                                            &
     3358                     ( 37.0_wp * ibit23 * adv_mom_5                            &
     3359                  +     7.0_wp * ibit22 * adv_mom_3                            &
     3360                  +              ibit21 * adv_mom_1                            &
     3361                     ) *                                                       &
     3362                                    ( w(k,j+1,i) + w(k,j,i)   )                &
     3363              -      (  8.0_wp * ibit23 * adv_mom_5                            &
     3364                  +              ibit22 * adv_mom_3                            &
     3365                     ) *                                                       &
     3366                                    ( w(k,j+2,i) + w(k,j-1,i) )                &
     3367              +      (           ibit23 * adv_mom_5                            &
     3368                     ) *                                                       &
     3369                                    ( w(k,j+3,i) + w(k,j-2,i) )                &
    33703370                                  )
    33713371
    3372           diss_n(k) = - ABS( v_comp(k) ) * (                                  &
    3373                      ( 10.0_wp * ibit23 * adv_mom_5                           &
    3374                   +     3.0_wp * ibit22 * adv_mom_3                           &
    3375                   +              ibit21 * adv_mom_1                           &
    3376                      ) *                                                      &
    3377                                     ( w(k,j+1,i) - w(k,j,i)   )               &
    3378               -      (  5.0_wp * ibit23 * adv_mom_5                           &
    3379                   +              ibit22 * adv_mom_3                           &
    3380                      ) *                                                      &
    3381                                    ( w(k,j+2,i)  - w(k,j-1,i) )               &
    3382               +      (           ibit23 * adv_mom_5                           &
    3383                      ) *                                                      &
    3384                                    ( w(k,j+3,i)  - w(k,j-2,i) )               &
     3372          diss_n(k) = - ABS( v_comp(k) ) * (                                   &
     3373                     ( 10.0_wp * ibit23 * adv_mom_5                            &
     3374                  +     3.0_wp * ibit22 * adv_mom_3                            &
     3375                  +              ibit21 * adv_mom_1                            &
     3376                     ) *                                                       &
     3377                                    ( w(k,j+1,i) - w(k,j,i)   )                &
     3378              -      (  5.0_wp * ibit23 * adv_mom_5                            &
     3379                  +              ibit22 * adv_mom_3                            &
     3380                     ) *                                                       &
     3381                                   ( w(k,j+2,i)  - w(k,j-1,i) )                &
     3382              +      (           ibit23 * adv_mom_5                            &
     3383                     ) *                                                       &
     3384                                   ( w(k,j+3,i)  - w(k,j-2,i) )                &
    33853385                                           )
    33863386       ENDDO
     
    33893389
    33903390          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
    3391           flux_r(k) = u_comp(k) * (                                           &
    3392                       37.0_wp * ( w(k,j,i+1) + w(k,j,i)   )                   &
    3393                     -  8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) )                   &
     3391          flux_r(k) = u_comp(k) * (                                            &
     3392                      37.0_wp * ( w(k,j,i+1) + w(k,j,i)   )                    &
     3393                    -  8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) )                    &
    33943394                    +           ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
    33953395
    3396           diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    3397                       10.0_wp * ( w(k,j,i+1) - w(k,j,i)   )                   &
    3398                     -  5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) )                   &
     3396          diss_r(k) = - ABS( u_comp(k) ) * (                                   &
     3397                      10.0_wp * ( w(k,j,i+1) - w(k,j,i)   )                    &
     3398                    -  5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) )                    &
    33993399                    +           ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
    34003400
    34013401          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
    3402           flux_n(k) = v_comp(k) * (                                           &
    3403                       37.0_wp * ( w(k,j+1,i) + w(k,j,i)   )                   &
    3404                     -  8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) )                   &
     3402          flux_n(k) = v_comp(k) * (                                            &
     3403                      37.0_wp * ( w(k,j+1,i) + w(k,j,i)   )                    &
     3404                    -  8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) )                    &
    34053405                    +           ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
    34063406
    3407           diss_n(k) = - ABS( v_comp(k) ) * (                                  &
    3408                       10.0_wp * ( w(k,j+1,i) - w(k,j,i)   )                   &
    3409                     -  5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) )                   &
     3407          diss_n(k) = - ABS( v_comp(k) ) * (                                   &
     3408                      10.0_wp * ( w(k,j+1,i) - w(k,j,i)   )                    &
     3409                    -  5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) )                    &
    34103410                    +           ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
    34113411       ENDDO
    34123412
    34133413!
    3414 !--    Now, compute vertical fluxes. Split loop into a part treating the 
     3414!--    Now, compute vertical fluxes. Split loop into a part treating the
    34153415!--    lowest grid points with indirect indexing, a main loop without
    34163416!--    indirect indexing, and a loop for the uppermost grip points with
    34173417!--    indirect indexing. This allows better vectorization for the main loop.
    3418 !--    First, compute the flux at model surface, which need has to be 
     3418!--    First, compute the flux at model surface, which need has to be
    34193419!--    calculated explicitly for the tendency at
    34203420!--    the first w-level. For topography wall this is done implicitely by
    3421 !--    advc_flags_m.
     3421!--    advc_flags_m. First, compute flux at lowest level, located at z=dz/2.
    34223422       k         = nzb + 1
    34233423       w_comp(k) = w(k,j,i) + w(k-1,j,i)
     
    34263426       diss_t(0) = -ABS(w_comp(k)) * rho_air(k)                                &
    34273427                 * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
    3428        
     3428
    34293429       DO  k = nzb+1, nzb+1
    34303430!
     
    34403440
    34413441          w_comp(k) = w(k+1,j,i) + w(k,j,i)
    3442           flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
    3443                      ( 37.0_wp * ibit26 * adv_mom_5                           &
    3444                   +     7.0_wp * ibit25 * adv_mom_3                           &
    3445                   +              ibit24 * adv_mom_1                           &
    3446                      ) *                                                      &
    3447                                 ( w(k+1,j,i)   + w(k,j,i)    )                &
    3448               -      (  8.0_wp * ibit26 * adv_mom_5                           &
    3449                   +              ibit25 * adv_mom_3                           &
    3450                      ) *                                                      &
    3451                                 ( w(k_pp,j,i)  + w(k-1,j,i)  )                &
    3452               +      (           ibit26 * adv_mom_5                           &
    3453                      ) *                                                      &
    3454                                 ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
     3442          flux_t(k) = w_comp(k) * rho_air(k+1) * (                             &
     3443                     ( 37.0_wp * ibit26 * adv_mom_5                            &
     3444                  +     7.0_wp * ibit25 * adv_mom_3                            &
     3445                  +              ibit24 * adv_mom_1                            &
     3446                     ) *                                                       &
     3447                                ( w(k+1,j,i)   + w(k,j,i)    )                 &
     3448              -      (  8.0_wp * ibit26 * adv_mom_5                            &
     3449                  +              ibit25 * adv_mom_3                            &
     3450                     ) *                                                       &
     3451                                ( w(k_pp,j,i)  + w(k-1,j,i)  )                 &
     3452              +      (           ibit26 * adv_mom_5                            &
     3453                     ) *                                                       &
     3454                                ( w(k_ppp,j,i) + w(k_mm,j,i) )                 &
    34553455                                                 )
    34563456
    3457           diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
    3458                      ( 10.0_wp * ibit26 * adv_mom_5                           &
    3459                   +     3.0_wp * ibit25 * adv_mom_3                           &
    3460                   +              ibit24 * adv_mom_1                           &
    3461                      ) *                                                      &
    3462                                 ( w(k+1,j,i)   - w(k,j,i)    )                &
    3463               -      (  5.0_wp * ibit26 * adv_mom_5                           &
    3464                   +              ibit25 * adv_mom_3                           &
    3465                      ) *                                                      &
    3466                                 ( w(k_pp,j,i)  - w(k-1,j,i)  )                &
    3467               +      (           ibit26 * adv_mom_5                           &
    3468                      ) *                                                      &
    3469                                 ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
     3457          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                    &
     3458                     ( 10.0_wp * ibit26 * adv_mom_5                            &
     3459                  +     3.0_wp * ibit25 * adv_mom_3                            &
     3460                  +              ibit24 * adv_mom_1                            &
     3461                     ) *                                                       &
     3462                                ( w(k+1,j,i)   - w(k,j,i)    )                 &
     3463              -      (  5.0_wp * ibit26 * adv_mom_5                            &
     3464                  +              ibit25 * adv_mom_3                            &
     3465                     ) *                                                       &
     3466                                ( w(k_pp,j,i)  - w(k-1,j,i)  )                 &
     3467              +      (           ibit26 * adv_mom_5                            &
     3468                     ) *                                                       &
     3469                                ( w(k_ppp,j,i) - w(k_mm,j,i) )                 &
    34703470                                                          )
    34713471       ENDDO
    3472        
     3472
    34733473       DO  k = nzb+2, nzt-2
    3474        
     3474
    34753475          ibit26 = REAL( IBITS(advc_flags_m(k,j,i),26,1), KIND = wp )
    34763476          ibit25 = REAL( IBITS(advc_flags_m(k,j,i),25,1), KIND = wp )
     
    34783478
    34793479          w_comp(k) = w(k+1,j,i) + w(k,j,i)
    3480           flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
    3481                      ( 37.0_wp * ibit26 * adv_mom_5                           &
    3482                   +     7.0_wp * ibit25 * adv_mom_3                           &
    3483                   +              ibit24 * adv_mom_1            &n