Changeset 4109


Ignore:
Timestamp:
Jul 22, 2019 5:00:34 PM (5 years ago)
Author:
suehring
Message:

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

Location:
palm/trunk/SOURCE
Files:
11 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/advec_ws.f90

    r4079 r4109  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! - Separate initialization of advection flags for momentum and scalars. In this
     23!   context, resort the bits and do some minor formatting.
     24! - Make flag initialization for scalars more flexible, introduce an
     25!   arguemnt list to indicate non-cyclic boundaries (required for decycled
     26!   scalars such as chemical species or aerosols)
     27! - Introduce extended 'degradation zones', where horizontal advection of
     28!   passive scalars is discretized by first-order scheme at all grid points
     29!   that in the vicinity of buildings (<= 3 grid points). Even though no
     30!   building is within the numerical stencil, first-order scheme is used.
     31!   At fourth and fifth grid point the order of the horizontal advection scheme
     32!   is successively upgraded.
     33!   These extended degradation zones are used to avoid stationary numerical
     34!   oscillations, which are responsible for high concentration maxima that may
     35!   appear under shear-free stable conditions.
     36! - Change interface for scalar advection routine.
     37! - Bugfix, avoid uninitialized value sk_num in vector version of scalar
     38!   advection
    2339!
    2440! Former revisions:
     
    106122! Change in file header (GPL part)
    107123! Implement advection for TKE-dissipation in case of RANS-TKE-e closure (TG)
    108 ! Allocate advc_flags_1/2 within ws_init_flags instead of init_grid
     124! Allocate advc_flags_m/2 within ws_init_flags instead of init_grid
    109125! Change argument list for exchange_horiz_2d_int (MS)
    110126!
     
    120136!
    121137! 2232 2017-05-30 17:47:52Z suehring
    122 ! Rename wall_flags_0 and wall_flags_00 into advc_flags_1 and advc_flags_2,
     138! Rename wall_flags_0 and wall_flags_00 into advc_flags_m and advc_flags_m,
    123139! respectively.
    124 ! Set advc_flags_1/2 on basis of wall_flags_0/00 instead of nzb_s/u/v/w_inner.
    125 ! Setting advc_flags_1/2 also for downward-facing walls
     140! Set advc_flags_m/2 on basis of wall_flags_0/00 instead of nzb_s/u/v/w_inner.
     141! Setting advc_flags_m/2 also for downward-facing walls
    126142!
    127143! 2200 2017-04-11 11:37:51Z suehring
     
    246262! vector version.
    247263! Degradation of the applied order of scheme is now steered by multiplying with
    248 ! Integer advc_flags_1. 2nd order scheme, WS3 and WS5 are calculated on each
     264! Integer advc_flags_m. 2nd order scheme, WS3 and WS5 are calculated on each
    249265! grid point and mulitplied with the appropriate flag.
    250266! 2nd order numerical dissipation term changed. Now the appropriate 2nd order
     
    267283!
    268284! 411 2009-12-11 12:31:43 Z suehring
     285!
     286!
     287!
     288! @author Matthias Suehring
     289!
    269290!
    270291! Description:
     
    282303!>
    283304!> @todo Implement monotonic flux limiter also for vector version.
     305!> @todo Move 3d arrays advc_flag, advc_flags_m from modules to advec_ws
     306!> @todo Move arrays flux_l_x from modules to advec_ws
    284307!------------------------------------------------------------------------------!
    285308 MODULE advec_ws
     
    299322
    300323    USE control_parameters,                                                    &
    301         ONLY:  humidity, loop_optimization, passive_scalar,                    &
    302                rans_tke_e, ws_scheme_mom, ws_scheme_sca,                       &
    303                momentum_advec, scalar_advec,                                   &
    304                bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,                 &
    305                bc_dirichlet_s, bc_radiation_l, bc_radiation_n,                 &
    306                bc_radiation_r, bc_radiation_s, intermediate_timestep_count,    &
    307                u_gtrans, v_gtrans, dt_3d
     324        ONLY:  air_chemistry,                                                  &
     325               bc_dirichlet_l,                                                 &
     326               bc_dirichlet_n,                                                 &
     327               bc_dirichlet_r,                                                 &
     328               bc_dirichlet_s,                                                 &
     329               bc_radiation_l,                                                 &
     330               bc_radiation_n,                                                 &
     331               bc_radiation_r,                                                 &
     332               bc_radiation_s,                                                 &
     333               humidity,                                                       &
     334               loop_optimization,                                              &
     335               passive_scalar,                                                 &
     336               rans_tke_e,                                                     &
     337               momentum_advec,                                                 &
     338               salsa,                                                          &
     339               scalar_advec,                                                   &
     340               intermediate_timestep_count,                                    &
     341               u_gtrans,                                                       &
     342               v_gtrans,                                                       &
     343               ws_scheme_mom,                                                  &
     344               ws_scheme_sca,                                                  &
     345               dt_3d
    308346
    309347    USE indices,                                                               &
    310         ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
    311                nzb, nzb_max, nzt, advc_flags_1, advc_flags_2, wall_flags_0
     348        ONLY:  advc_flags_m,                                                   &
     349               advc_flags_s,                                                   &
     350               nbgp,                                                           &
     351               nx,                                                             &
     352               nxl,                                                            &
     353               nxlg,                                                           &
     354               nxlu,                                                           &
     355               nxr,                                                            &
     356               nxrg,                                                           &
     357               ny,                                                             &
     358               nyn,                                                            &
     359               nyng,                                                           &
     360               nys,                                                            &
     361               nysg,                                                           &
     362               nysv,                                                           &
     363               nzb,                                                            &
     364               nzb_max,                                                        &
     365               nzt,                                                            &
     366               wall_flags_0
    312367
    313368    USE grid_variables,                                                        &
     
    338393    PRIVATE
    339394    PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, ws_init,          &
    340              ws_init_flags, ws_statistics
     395             ws_init_flags_momentum, ws_init_flags_scalar, ws_statistics
    341396
    342397    INTERFACE ws_init
    343398       MODULE PROCEDURE ws_init
    344     END INTERFACE ws_init
    345 
    346     INTERFACE ws_init_flags
    347        MODULE PROCEDURE ws_init_flags
    348     END INTERFACE ws_init_flags
     399    END INTERFACE ws_init         
     400             
     401    INTERFACE ws_init_flags_momentum
     402       MODULE PROCEDURE ws_init_flags_momentum
     403    END INTERFACE ws_init_flags_momentum
     404   
     405    INTERFACE ws_init_flags_scalar
     406       MODULE PROCEDURE ws_init_flags_scalar
     407    END INTERFACE ws_init_flags_scalar
    349408
    350409    INTERFACE ws_statistics
     
    490549! Description:
    491550! ------------
    492 !> Initialization of flags for WS-scheme used to degrade the order of the scheme
    493 !> near walls.
     551!> Initialization of flags to control the order of the advection scheme near
     552!> solid walls and non-cyclic inflow boundaries, where the order is sucessively
     553!> degraded.
    494554!------------------------------------------------------------------------------!
    495     SUBROUTINE ws_init_flags
     555    SUBROUTINE ws_init_flags_momentum
    496556
    497557
     
    505565       LOGICAL      ::  flag_set !< steering variable for advection flags
    506566   
    507        ALLOCATE( advc_flags_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                  &
    508                  advc_flags_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    509        advc_flags_1 = 0
    510        advc_flags_2 = 0
    511        IF ( scalar_advec == 'ws-scheme' )  THEN
    512 !
    513 !--       Set flags to steer the degradation of the advection scheme in advec_ws
    514 !--       near topography, inflow- and outflow boundaries as well as bottom and
    515 !--       top of model domain. advc_flags_1 remains zero for all non-prognostic
    516 !--       grid points.
    517           DO  i = nxl, nxr
    518              DO  j = nys, nyn
    519                 DO  k = nzb+1, nzt
    520 !
    521 !--                scalar - x-direction
    522 !--                WS1 (0), WS3 (1), WS5 (2)
    523                    IF ( ( .NOT. BTEST(wall_flags_0(k,j,i+1),0)                 &
    524                     .OR.  .NOT. BTEST(wall_flags_0(k,j,i+2),0)                 &   
    525                     .OR.  .NOT. BTEST(wall_flags_0(k,j,i-1),0) )               &
    526                       .OR.  ( ( bc_dirichlet_l .OR. bc_radiation_l )           &
    527                             .AND.  i == nxl   )                                &
    528                       .OR.  ( ( bc_dirichlet_r .OR. bc_radiation_r )           &
    529                             .AND.  i == nxr   ) )                              &
    530                    THEN
    531                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 0 )
    532                    ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+3),0)             &
    533                        .AND.        BTEST(wall_flags_0(k,j,i+1),0)             &
    534                        .AND.        BTEST(wall_flags_0(k,j,i+2),0)             &
    535                        .AND.        BTEST(wall_flags_0(k,j,i-1),0)             &
    536                             )                       .OR.                       &
    537                             ( .NOT. BTEST(wall_flags_0(k,j,i-2),0)             &
    538                        .AND.        BTEST(wall_flags_0(k,j,i+1),0)             &
    539                        .AND.        BTEST(wall_flags_0(k,j,i+2),0)             &
    540                        .AND.        BTEST(wall_flags_0(k,j,i-1),0)             &
    541                             )                                                  &
    542                                                     .OR.                       &
    543                             ( ( bc_dirichlet_r .OR. bc_radiation_r )           &
    544                               .AND. i == nxr-1 )    .OR.                       &
    545                             ( ( bc_dirichlet_l .OR. bc_radiation_l )           &
    546                               .AND. i == nxlu  ) )                             & ! why not nxl+1
    547                    THEN
    548                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 1 )
    549                    ELSEIF ( BTEST(wall_flags_0(k,j,i+1),0)                     &
    550                       .AND. BTEST(wall_flags_0(k,j,i+2),0)                     &
    551                       .AND. BTEST(wall_flags_0(k,j,i+3),0)                     &
    552                       .AND. BTEST(wall_flags_0(k,j,i-1),0)                     &
    553                       .AND. BTEST(wall_flags_0(k,j,i-2),0) )                   &
    554                    THEN
    555                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 2 )
    556                    ENDIF
    557 !
    558 !--                scalar - y-direction
    559 !--                WS1 (3), WS3 (4), WS5 (5)
    560                    IF ( ( .NOT. BTEST(wall_flags_0(k,j+1,i),0)                 &
    561                     .OR.  .NOT. BTEST(wall_flags_0(k,j+2,i),0)                 &   
    562                     .OR.  .NOT. BTEST(wall_flags_0(k,j-1,i),0))                &
    563                       .OR.  ( ( bc_dirichlet_s .OR. bc_radiation_s )           &
    564                             .AND.  j == nys   )                                &
    565                       .OR.  ( ( bc_dirichlet_n .OR. bc_radiation_n )           &
    566                             .AND.  j == nyn   ) )                              &
    567                    THEN
    568                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 3 )
    569 !
    570 !--                WS3
    571                    ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+3,i),0)             &
    572                        .AND.        BTEST(wall_flags_0(k,j+1,i),0)             &
    573                        .AND.        BTEST(wall_flags_0(k,j+2,i),0)             &
    574                        .AND.        BTEST(wall_flags_0(k,j-1,i),0)             &
    575                             )                       .OR.                       &
    576                             ( .NOT. BTEST(wall_flags_0(k,j-2,i),0)             &
    577                        .AND.        BTEST(wall_flags_0(k,j+1,i),0)             &
    578                        .AND.        BTEST(wall_flags_0(k,j+2,i),0)             &
    579                        .AND.        BTEST(wall_flags_0(k,j-1,i),0)             &
    580                             )                                                  &
    581                                                     .OR.                       &
    582                             ( ( bc_dirichlet_s .OR. bc_radiation_s )           &
    583                               .AND. j == nysv  )    .OR.                       & ! why not nys+1
    584                             ( ( bc_dirichlet_n .OR. bc_radiation_n )           &
    585                               .AND. j == nyn-1 ) )                             &         
    586                    THEN
    587                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 4 )
    588 !
    589 !--                WS5
    590                    ELSEIF ( BTEST(wall_flags_0(k,j+1,i),0)                     &
    591                       .AND. BTEST(wall_flags_0(k,j+2,i),0)                     &
    592                       .AND. BTEST(wall_flags_0(k,j+3,i),0)                     &
    593                       .AND. BTEST(wall_flags_0(k,j-1,i),0)                     &
    594                       .AND. BTEST(wall_flags_0(k,j-2,i),0) )                   &
    595                    THEN
    596                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 5 )
    597                    ENDIF
    598 !
    599 !--                scalar - z-direction
    600 !--                WS1 (6), WS3 (7), WS5 (8)
    601                    IF ( k == nzb+1 )  THEN
    602                       k_mm = nzb
    603                    ELSE
    604                       k_mm = k - 2
    605                    ENDIF
    606                    IF ( k > nzt-1 )  THEN
    607                       k_pp = nzt+1
    608                    ELSE
    609                       k_pp = k + 2
    610                    ENDIF
    611                    IF ( k > nzt-2 )  THEN
    612                       k_ppp = nzt+1
    613                    ELSE
    614                       k_ppp = k + 3
    615                    ENDIF
    616 
    617                    flag_set = .FALSE.
    618                    IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),0)  .AND.            &
    619                               BTEST(wall_flags_0(k,j,i),0)    .OR.             &
    620                         .NOT. BTEST(wall_flags_0(k_pp,j,i),0) .AND.            &                             
    621                               BTEST(wall_flags_0(k,j,i),0)    .OR.             &
    622                         k == nzt )                                             &
    623                    THEN
    624                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 6 )
    625                       flag_set = .TRUE.
    626                    ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),0)    .OR.    &
    627                               .NOT. BTEST(wall_flags_0(k_ppp,j,i),0) ) .AND.   &
    628                                   BTEST(wall_flags_0(k-1,j,i),0)  .AND.        &
    629                                   BTEST(wall_flags_0(k,j,i),0)    .AND.        &
    630                                   BTEST(wall_flags_0(k+1,j,i),0)  .AND.        &
    631                                   BTEST(wall_flags_0(k_pp,j,i),0) .OR.         &   
    632                             k == nzt - 1 )                                     &
    633                    THEN
    634                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 7 )
    635                       flag_set = .TRUE.
    636                    ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),0)                    &
    637                      .AND.  BTEST(wall_flags_0(k-1,j,i),0)                     &
    638                      .AND.  BTEST(wall_flags_0(k,j,i),0)                       &
    639                      .AND.  BTEST(wall_flags_0(k+1,j,i),0)                     &
    640                      .AND.  BTEST(wall_flags_0(k_pp,j,i),0)                    &
    641                      .AND.  BTEST(wall_flags_0(k_ppp,j,i),0)                   &
    642                      .AND. .NOT. flag_set )                                    &
    643                    THEN
    644                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 8 )
    645                    ENDIF
    646 
    647                 ENDDO
     567       advc_flags_m = 0
     568
     569!
     570!--    Set advc_flags_m to steer the degradation of the advection scheme in advec_ws
     571!--    near topography, inflow- and outflow boundaries as well as bottom and
     572!--    top of model domain. advc_flags_m remains zero for all non-prognostic
     573!--    grid points.
     574!--    u-component
     575       DO  i = nxl, nxr
     576          DO  j = nys, nyn
     577             DO  k = nzb+1, nzt
     578
     579!--             At first, set flags to WS1.
     580!--             Since fluxes are swapped in advec_ws.f90, this is necessary to
     581!--             in order to handle the left/south flux.
     582!--             near vertical walls.
     583                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 )
     584                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 )
     585!
     586!--             u component - x-direction
     587!--             WS1 (0), WS3 (1), WS5 (2)
     588                IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),1)  .OR.                &
     589                         ( ( bc_dirichlet_l .OR. bc_radiation_l )              &
     590                           .AND. i <= nxlu  )    .OR.                          &
     591                         ( ( bc_dirichlet_r .OR. bc_radiation_r )              &
     592                           .AND. i == nxr   ) )                                &
     593                THEN                                                           
     594                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 )     
     595                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),1)  .AND.         &
     596                                 BTEST(wall_flags_0(k,j,i+1),1)  .OR.          &
     597                           .NOT. BTEST(wall_flags_0(k,j,i-1),1) )              &
     598                                                     .OR.                      &
     599                         ( ( bc_dirichlet_r .OR. bc_radiation_r )              &
     600                           .AND. i == nxr-1 )    .OR.                          &
     601                         ( ( bc_dirichlet_l .OR. bc_radiation_l )              &
     602                           .AND. i == nxlu+1) )                                &
     603                THEN                                                           
     604                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 1 )       
     605!                                                                             
     606!--                Clear flag for WS1                                         
     607                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 )       
     608                ELSEIF ( BTEST(wall_flags_0(k,j,i+1),1)  .AND.                 &
     609                         BTEST(wall_flags_0(k,j,i+2),1)  .AND.                 &
     610                         BTEST(wall_flags_0(k,j,i-1),1) )                      &
     611                THEN                                                           
     612                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 2 )       
     613!                                                                             
     614!--                Clear flag for WS1                                         
     615                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 )       
     616                ENDIF                                                         
     617!                                                                             
     618!--             u component - y-direction                                     
     619!--             WS1 (3), WS3 (4), WS5 (5)                                     
     620                IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),1)   .OR.               &
     621                         ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
     622                           .AND. j == nys   )    .OR.                          &
     623                         ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
     624                           .AND. j == nyn   ) )                                &
     625                THEN                                                           
     626                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 )       
     627                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),1)  .AND.         &
     628                                 BTEST(wall_flags_0(k,j+1,i),1)  .OR.          &
     629                           .NOT. BTEST(wall_flags_0(k,j-1,i),1) )              &
     630                                                     .OR.                      &
     631                         ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
     632                           .AND. j == nysv  )    .OR.                          &
     633                         ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
     634                           .AND. j == nyn-1 ) )                                &
     635                THEN                                                           
     636                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 4 )       
     637!                                                                             
     638!--                Clear flag for WS1                                         
     639                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 )       
     640                ELSEIF ( BTEST(wall_flags_0(k,j+1,i),1)  .AND.                 &
     641                         BTEST(wall_flags_0(k,j+2,i),1)  .AND.                 &
     642                         BTEST(wall_flags_0(k,j-1,i),1) )                      &
     643                THEN                                                           
     644                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 5 )       
     645!                                                                             
     646!--                Clear flag for WS1                                         
     647                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 )       
     648                ENDIF                                                         
     649!                                                                             
     650!--             u component - z-direction                                     
     651!--             WS1 (6), WS3 (7), WS5 (8)                                     
     652                IF ( k == nzb+1 )  THEN                                       
     653                   k_mm = nzb                                                 
     654                ELSE                                                           
     655                   k_mm = k - 2                                               
     656                ENDIF                                                         
     657                IF ( k > nzt-1 )  THEN                                         
     658                   k_pp = nzt+1                                               
     659                ELSE                                                           
     660                   k_pp = k + 2                                               
     661                ENDIF                                                         
     662                IF ( k > nzt-2 )  THEN                                         
     663                   k_ppp = nzt+1                                               
     664                ELSE                                                           
     665                   k_ppp = k + 3                                               
     666                ENDIF                                                         
     667                                                                               
     668                flag_set = .FALSE.                                             
     669                IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),1)  .AND.               &
     670                           BTEST(wall_flags_0(k,j,i),1)    .OR.                &
     671                     .NOT. BTEST(wall_flags_0(k_pp,j,i),1) .AND.               &                             
     672                           BTEST(wall_flags_0(k,j,i),1)    .OR.                &
     673                     k == nzt )                                                &
     674                THEN                                                           
     675                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 6 )       
     676                   flag_set = .TRUE.                                           
     677                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),1)    .OR.       &
     678                           .NOT. BTEST(wall_flags_0(k_ppp,j,i),1) ) .AND.      &
     679                               BTEST(wall_flags_0(k-1,j,i),1)  .AND.           &
     680                               BTEST(wall_flags_0(k,j,i),1)    .AND.           &
     681                               BTEST(wall_flags_0(k+1,j,i),1)  .AND.           &
     682                               BTEST(wall_flags_0(k_pp,j,i),1) .OR.            &
     683                               k == nzt - 1 )                                  &
     684                THEN                                                           
     685                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 )       
     686                   flag_set = .TRUE.                                           
     687                ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),1)  .AND.                &
     688                         BTEST(wall_flags_0(k-1,j,i),1)   .AND.                &
     689                         BTEST(wall_flags_0(k,j,i),1)     .AND.                &
     690                         BTEST(wall_flags_0(k+1,j,i),1)   .AND.                &
     691                         BTEST(wall_flags_0(k_pp,j,i),1)  .AND.                &
     692                         BTEST(wall_flags_0(k_ppp,j,i),1) .AND.                &
     693                         .NOT. flag_set )                                      &
     694                THEN
     695                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 8 )
     696                ENDIF
     697
    648698             ENDDO
    649699          ENDDO
    650        ENDIF
    651 
    652        IF ( momentum_advec == 'ws-scheme' )  THEN
    653 !
    654 !--       Set advc_flags_1 to steer the degradation of the advection scheme in advec_ws
    655 !--       near topography, inflow- and outflow boundaries as well as bottom and
    656 !--       top of model domain. advc_flags_1 remains zero for all non-prognostic
    657 !--       grid points.
    658           DO  i = nxl, nxr
    659              DO  j = nys, nyn
    660                 DO  k = nzb+1, nzt
    661 
    662 !--                At first, set flags to WS1.
    663 !--                Since fluxes are swapped in advec_ws.f90, this is necessary to
    664 !--                in order to handle the left/south flux.
    665 !--                near vertical walls.
    666                    advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 9 )
    667                    advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 12 )
    668 !
    669 !--                u component - x-direction
    670 !--                WS1 (9), WS3 (10), WS5 (11)
    671                    IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),1)  .OR.             &
    672                             ( ( bc_dirichlet_l .OR. bc_radiation_l )           &
    673                               .AND. i <= nxlu  )    .OR.                       &
    674                             ( ( bc_dirichlet_r .OR. bc_radiation_r )           &
    675                               .AND. i == nxr   ) )                             &
    676                    THEN
    677                        advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 9 )
    678                    ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),1)  .AND.      &
    679                                     BTEST(wall_flags_0(k,j,i+1),1)  .OR.       &
    680                               .NOT. BTEST(wall_flags_0(k,j,i-1),1) )           &
    681                                                         .OR.                   &
    682                             ( ( bc_dirichlet_r .OR. bc_radiation_r )           &
    683                               .AND. i == nxr-1 )    .OR.                       &
    684                             ( ( bc_dirichlet_l .OR. bc_radiation_l )           &
    685                               .AND. i == nxlu+1) )                             &
    686                    THEN
    687                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 10 )
    688 !
    689 !--                   Clear flag for WS1
    690                       advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 9 )
    691                    ELSEIF ( BTEST(wall_flags_0(k,j,i+1),1)  .AND.              &
    692                             BTEST(wall_flags_0(k,j,i+2),1)  .AND.              &
    693                             BTEST(wall_flags_0(k,j,i-1),1) )                   &
    694                    THEN   
    695                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 11 )
    696 !
    697 !--                   Clear flag for WS1
    698                       advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 9 )
    699                    ENDIF
    700 !
    701 !--                u component - y-direction
    702 !--                WS1 (12), WS3 (13), WS5 (14)
    703                    IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),1)   .OR.            &
    704                             ( ( bc_dirichlet_s .OR. bc_radiation_s )           &
    705                               .AND. j == nys   )    .OR.                       &
    706                             ( ( bc_dirichlet_n .OR. bc_radiation_n )           &
    707                               .AND. j == nyn   ) )                             &
    708                    THEN
    709                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 12 )
    710                    ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),1)  .AND.      &
    711                                     BTEST(wall_flags_0(k,j+1,i),1)  .OR.       &
    712                               .NOT. BTEST(wall_flags_0(k,j-1,i),1) )           &
    713                                                         .OR.                   &
    714                             ( ( bc_dirichlet_s .OR. bc_radiation_s )           &
    715                               .AND. j == nysv  )    .OR.                       &
    716                             ( ( bc_dirichlet_n .OR. bc_radiation_n )           &
    717                               .AND. j == nyn-1 ) )                             &
    718                    THEN
    719                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 13 )
    720 !
    721 !--                   Clear flag for WS1
    722                       advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 12 )
    723                    ELSEIF ( BTEST(wall_flags_0(k,j+1,i),1)  .AND.              &
    724                             BTEST(wall_flags_0(k,j+2,i),1)  .AND.              &
    725                             BTEST(wall_flags_0(k,j-1,i),1) )                   &
    726                    THEN
    727                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 14 )
    728 !
    729 !--                   Clear flag for WS1
    730                       advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 12 )
    731                    ENDIF
    732 !
    733 !--                u component - z-direction
    734 !--                WS1 (15), WS3 (16), WS5 (17)
    735                    IF ( k == nzb+1 )  THEN
    736                       k_mm = nzb
    737                    ELSE
    738                       k_mm = k - 2
    739                    ENDIF
    740                    IF ( k > nzt-1 )  THEN
    741                       k_pp = nzt+1
    742                    ELSE
    743                       k_pp = k + 2
    744                    ENDIF
    745                    IF ( k > nzt-2 )  THEN
    746                       k_ppp = nzt+1
    747                    ELSE
    748                       k_ppp = k + 3
    749                    ENDIF                   
    750 
    751                    flag_set = .FALSE.
    752                    IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),1)  .AND.            &
    753                               BTEST(wall_flags_0(k,j,i),1)    .OR.             &
    754                         .NOT. BTEST(wall_flags_0(k_pp,j,i),1) .AND.            &                             
    755                               BTEST(wall_flags_0(k,j,i),1)    .OR.             &
    756                         k == nzt )                                             &
    757                    THEN
    758                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 15 )
    759                       flag_set = .TRUE.                     
    760                    ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),1)    .OR.    &
    761                               .NOT. BTEST(wall_flags_0(k_ppp,j,i),1) ) .AND.   &
    762                                   BTEST(wall_flags_0(k-1,j,i),1)  .AND.        &
    763                                   BTEST(wall_flags_0(k,j,i),1)    .AND.        &
    764                                   BTEST(wall_flags_0(k+1,j,i),1)  .AND.        &
    765                                   BTEST(wall_flags_0(k_pp,j,i),1) .OR.         &
    766                                   k == nzt - 1 )                               &
    767                    THEN
    768                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 16 )
    769                       flag_set = .TRUE.
    770                    ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),1)  .AND.             &
    771                             BTEST(wall_flags_0(k-1,j,i),1)   .AND.             &
    772                             BTEST(wall_flags_0(k,j,i),1)     .AND.             &
    773                             BTEST(wall_flags_0(k+1,j,i),1)   .AND.             &
    774                             BTEST(wall_flags_0(k_pp,j,i),1)  .AND.             &
    775                             BTEST(wall_flags_0(k_ppp,j,i),1) .AND.             &
    776                             .NOT. flag_set )                                   &
    777                    THEN
    778                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 17 )
    779                    ENDIF
    780 
    781                 ENDDO
     700       ENDDO
     701!
     702!--    v-component
     703       DO  i = nxl, nxr
     704          DO  j = nys, nyn
     705             DO  k = nzb+1, nzt
     706!
     707!--             At first, set flags to WS1.
     708!--             Since fluxes are swapped in advec_ws.f90, this is necessary to
     709!--             in order to handle the left/south flux.
     710                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9  )
     711                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 )
     712!
     713!--             v component - x-direction
     714!--             WS1 (9), WS3 (10), WS5 (11)
     715                IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),2)  .OR.                &
     716                         ( ( bc_dirichlet_l .OR. bc_radiation_l )              &
     717                           .AND. i == nxl   )    .OR.                          &
     718                         ( ( bc_dirichlet_r .OR. bc_radiation_r )              &
     719                           .AND. i == nxr   ) )                                &
     720               THEN                                                           
     721                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 )       
     722!                                                                             
     723!--             WS3                                                           
     724                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),2)   .AND.        &
     725                                 BTEST(wall_flags_0(k,j,i+1),2) ) .OR.         &
     726                           .NOT. BTEST(wall_flags_0(k,j,i-1),2)                &
     727                                                 .OR.                          &
     728                         ( ( bc_dirichlet_r .OR. bc_radiation_r )              &
     729                           .AND. i == nxr-1 )    .OR.                          &
     730                         ( ( bc_dirichlet_l .OR. bc_radiation_l )              &
     731                           .AND. i == nxlu  ) )                                &
     732                THEN                                                           
     733                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 10 )     
     734!                                                                             
     735!--                Clear flag for WS1                                         
     736                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 )       
     737                ELSEIF ( BTEST(wall_flags_0(k,j,i+1),2)  .AND.                 &
     738                         BTEST(wall_flags_0(k,j,i+2),2)  .AND.                 &
     739                         BTEST(wall_flags_0(k,j,i-1),2) )                      &
     740                THEN                                                           
     741                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 11 )     
     742!                                                                             
     743!--                Clear flag for WS1                                         
     744                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 )       
     745                ENDIF                                                         
     746!                                                                             
     747!--             v component - y-direction                                     
     748!--             WS1 (12), WS3 (13), WS5 (14)                                   
     749                IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),2) .OR.                 &
     750                         ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
     751                           .AND. j <= nysv  )    .OR.                          &
     752                         ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
     753                           .AND. j == nyn   ) )                                &
     754                THEN                                                           
     755                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 )     
     756                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),2)  .AND.         &
     757                                 BTEST(wall_flags_0(k,j+1,i),2)  .OR.          &
     758                           .NOT. BTEST(wall_flags_0(k,j-1,i),2) )              &
     759                                                     .OR.                      &
     760                         ( (  bc_dirichlet_s .OR. bc_radiation_s )             &
     761                           .AND. j == nysv+1)    .OR.                          &
     762                         ( (  bc_dirichlet_n .OR. bc_radiation_n )             &
     763                           .AND. j == nyn-1 ) )                                &
     764                THEN                                                           
     765                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 13 )     
     766!                                                                             
     767!--                Clear flag for WS1                                         
     768                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 )     
     769                ELSEIF ( BTEST(wall_flags_0(k,j+1,i),2)  .AND.                 &
     770                         BTEST(wall_flags_0(k,j+2,i),2)  .AND.                 &
     771                         BTEST(wall_flags_0(k,j-1,i),2) )                      &
     772                THEN
     773                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 14 )
     774!
     775!--                Clear flag for WS1
     776                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 )
     777                ENDIF
     778!
     779!--             v component - z-direction
     780!--             WS1 (15), WS3 (16), WS5 (17)
     781                IF ( k == nzb+1 )  THEN
     782                   k_mm = nzb
     783                ELSE
     784                   k_mm = k - 2
     785                ENDIF
     786                IF ( k > nzt-1 )  THEN
     787                   k_pp = nzt+1
     788                ELSE
     789                   k_pp = k + 2
     790                ENDIF
     791                IF ( k > nzt-2 )  THEN
     792                   k_ppp = nzt+1
     793                ELSE
     794                   k_ppp = k + 3
     795                ENDIF 
     796               
     797                flag_set = .FALSE.
     798                IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),2)  .AND.               &
     799                           BTEST(wall_flags_0(k,j,i),2)    .OR.                &
     800                     .NOT. BTEST(wall_flags_0(k_pp,j,i),2) .AND.               &                             
     801                           BTEST(wall_flags_0(k,j,i),2)    .OR.                &
     802                     k == nzt )                                                &
     803                THEN                                                           
     804                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 15 )     
     805                   flag_set = .TRUE.                                           
     806                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),2)    .OR.       &
     807                           .NOT. BTEST(wall_flags_0(k_ppp,j,i),2) ) .AND.      &
     808                               BTEST(wall_flags_0(k-1,j,i),2)  .AND.           &
     809                               BTEST(wall_flags_0(k,j,i),2)    .AND.           &
     810                               BTEST(wall_flags_0(k+1,j,i),2)  .AND.           &
     811                               BTEST(wall_flags_0(k_pp,j,i),2)  .OR.           &
     812                               k == nzt - 1 )                                  &
     813                THEN                                                           
     814                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 )     
     815                   flag_set = .TRUE.                                           
     816                ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),2)  .AND.                &
     817                         BTEST(wall_flags_0(k-1,j,i),2)   .AND.                &
     818                         BTEST(wall_flags_0(k,j,i),2)     .AND.                &
     819                         BTEST(wall_flags_0(k+1,j,i),2)   .AND.                &
     820                         BTEST(wall_flags_0(k_pp,j,i),2)  .AND.                &
     821                         BTEST(wall_flags_0(k_ppp,j,i),2) .AND.                &
     822                         .NOT. flag_set )                                      &
     823                THEN
     824                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 17 )
     825                ENDIF
     826
    782827             ENDDO
    783828          ENDDO
    784 
    785           DO  i = nxl, nxr
    786              DO  j = nys, nyn
    787                 DO  k = nzb+1, nzt
    788 !
    789 !--                At first, set flags to WS1.
    790 !--                Since fluxes are swapped in advec_ws.f90, this is necessary to
    791 !--                in order to handle the left/south flux.
    792                    advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 18 )
    793                    advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 21 )
    794 !
    795 !--                v component - x-direction
    796 !--                WS1 (18), WS3 (19), WS5 (20)
    797                    IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),2)  .OR.             &
    798                             ( ( bc_dirichlet_l .OR. bc_radiation_l )           &
    799                               .AND. i == nxl   )    .OR.                       &
    800                             ( ( bc_dirichlet_r .OR. bc_radiation_r )           &
    801                               .AND. i == nxr   ) )                             &
    802                   THEN
    803                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 18 )
    804 !
    805 !--                WS3
    806                    ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),2)   .AND.     &
    807                                     BTEST(wall_flags_0(k,j,i+1),2) ) .OR.      &
    808                               .NOT. BTEST(wall_flags_0(k,j,i-1),2)             &
    809                                                     .OR.                       &
    810                             ( ( bc_dirichlet_r .OR. bc_radiation_r )           &
    811                               .AND. i == nxr-1 )    .OR.                       &
    812                             ( ( bc_dirichlet_l .OR. bc_radiation_l )           &
    813                               .AND. i == nxlu  ) )                             &
    814                    THEN
    815                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 19 )
    816 !
    817 !--                   Clear flag for WS1
    818                       advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 18 )
    819                    ELSEIF ( BTEST(wall_flags_0(k,j,i+1),2)  .AND.              &
    820                             BTEST(wall_flags_0(k,j,i+2),2)  .AND.              &
    821                             BTEST(wall_flags_0(k,j,i-1),2) )                   &
    822                    THEN
    823                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 20 )
    824 !
    825 !--                   Clear flag for WS1
    826                       advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 18 )
    827                    ENDIF
    828 !
    829 !--                v component - y-direction
    830 !--                WS1 (21), WS3 (22), WS5 (23)
    831                    IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),2) .OR.              &
    832                             ( ( bc_dirichlet_s .OR. bc_radiation_s )           &
    833                               .AND. j <= nysv  )    .OR.                       &
    834                             ( ( bc_dirichlet_n .OR. bc_radiation_n )           &
    835                               .AND. j == nyn   ) )                             &
    836                    THEN
    837                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 21 )
    838                    ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),2)  .AND.      &
    839                                     BTEST(wall_flags_0(k,j+1,i),2)  .OR.       &
    840                               .NOT. BTEST(wall_flags_0(k,j-1,i),2) )           &
    841                                                         .OR.                   &
    842                             ( (  bc_dirichlet_s .OR. bc_radiation_s )          &
    843                               .AND. j == nysv+1)    .OR.                       &
    844                             ( (  bc_dirichlet_n .OR. bc_radiation_n )          &
    845                               .AND. j == nyn-1 ) )                             &
    846                    THEN
    847                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 22 )
    848 !
    849 !--                   Clear flag for WS1
    850                       advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 21 )
    851                    ELSEIF ( BTEST(wall_flags_0(k,j+1,i),2)  .AND.              &
    852                             BTEST(wall_flags_0(k,j+2,i),2)  .AND.              &
    853                             BTEST(wall_flags_0(k,j-1,i),2) )                   &
    854                    THEN
    855                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 23 )
    856 !
    857 !--                   Clear flag for WS1
    858                       advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 21 )
    859                    ENDIF
    860 !
    861 !--                v component - z-direction
    862 !--                WS1 (24), WS3 (25), WS5 (26)
    863                    IF ( k == nzb+1 )  THEN
    864                       k_mm = nzb
    865                    ELSE
    866                       k_mm = k - 2
    867                    ENDIF
    868                    IF ( k > nzt-1 )  THEN
    869                       k_pp = nzt+1
    870                    ELSE
    871                       k_pp = k + 2
    872                    ENDIF
    873                    IF ( k > nzt-2 )  THEN
    874                       k_ppp = nzt+1
    875                    ELSE
    876                       k_ppp = k + 3
    877                    ENDIF 
    878                    
    879                    flag_set = .FALSE.
    880                    IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),2)  .AND.            &
    881                               BTEST(wall_flags_0(k,j,i),2)    .OR.             &
    882                         .NOT. BTEST(wall_flags_0(k_pp,j,i),2) .AND.            &                             
    883                               BTEST(wall_flags_0(k,j,i),2)    .OR.             &
    884                         k == nzt )                                             &
    885                    THEN
    886                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 24 )
    887                       flag_set = .TRUE.
    888                    ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),2)    .OR.    &
    889                               .NOT. BTEST(wall_flags_0(k_ppp,j,i),2) ) .AND.   &
    890                                   BTEST(wall_flags_0(k-1,j,i),2)  .AND.        &
    891                                   BTEST(wall_flags_0(k,j,i),2)    .AND.        &
    892                                   BTEST(wall_flags_0(k+1,j,i),2)  .AND.        &
    893                                   BTEST(wall_flags_0(k_pp,j,i),2)  .OR.        &
    894                                   k == nzt - 1 )                               &
    895                    THEN
    896                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 25 )
    897                       flag_set = .TRUE.
    898                    ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),2)  .AND.             &
    899                             BTEST(wall_flags_0(k-1,j,i),2)   .AND.             &
    900                             BTEST(wall_flags_0(k,j,i),2)     .AND.             &
    901                             BTEST(wall_flags_0(k+1,j,i),2)   .AND.             &
    902                             BTEST(wall_flags_0(k_pp,j,i),2)  .AND.             &
    903                             BTEST(wall_flags_0(k_ppp,j,i),2) .AND.             &
    904                             .NOT. flag_set )                                   &
    905                    THEN
    906                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 26 )
    907                    ENDIF
    908 
    909                 ENDDO
     829       ENDDO
     830!
     831!--    w - component
     832       DO  i = nxl, nxr
     833          DO  j = nys, nyn
     834             DO  k = nzb+1, nzt
     835!
     836!--             At first, set flags to WS1.
     837!--             Since fluxes are swapped in advec_ws.f90, this is necessary to
     838!--             in order to handle the left/south flux.
     839                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 )
     840                advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 )
     841!
     842!--             w component - x-direction
     843!--             WS1 (18), WS3 (19), WS5 (20)
     844                IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),3) .OR.                 &
     845                         ( (  bc_dirichlet_l .OR. bc_radiation_l )             &
     846                           .AND. i == nxl   )    .OR.                          &
     847                         ( (  bc_dirichlet_r .OR. bc_radiation_r )             &
     848                           .AND. i == nxr   ) )                                &
     849                THEN                                                           
     850                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 )     
     851                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),3)  .AND.         &
     852                                 BTEST(wall_flags_0(k,j,i+1),3)  .OR.          &
     853                           .NOT. BTEST(wall_flags_0(k,j,i-1),3) )              &
     854                                                     .OR.                      &
     855                         ( ( bc_dirichlet_r .OR. bc_radiation_r )              &
     856                           .AND. i == nxr-1 )    .OR.                          &
     857                         ( ( bc_dirichlet_l .OR.  bc_radiation_l )             &
     858                           .AND. i == nxlu  ) )                                &
     859                THEN                                                           
     860                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 19 )     
     861!                                                                             
     862!--                Clear flag for WS1                                         
     863                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 )     
     864                ELSEIF ( BTEST(wall_flags_0(k,j,i+1),3)  .AND.                 &
     865                         BTEST(wall_flags_0(k,j,i+2),3)  .AND.                 &
     866                         BTEST(wall_flags_0(k,j,i-1),3) )                      &
     867                THEN                                                           
     868                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i),20 )       
     869!                                                                             
     870!--                Clear flag for WS1                                         
     871                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 )     
     872                ENDIF                                                         
     873!                                                                             
     874!--             w component - y-direction                                     
     875!--             WS1 (21), WS3 (22), WS5 (23)                                   
     876                IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),3) .OR.                 &
     877                         ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
     878                           .AND. j == nys   )    .OR.                          &
     879                         ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
     880                           .AND. j == nyn   ) )                                &
     881                THEN                                                           
     882                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 )     
     883                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),3)  .AND.         &
     884                                 BTEST(wall_flags_0(k,j+1,i),3)  .OR.          &
     885                           .NOT. BTEST(wall_flags_0(k,j-1,i),3) )              &
     886                                                     .OR.                      &
     887                         ( ( bc_dirichlet_s .OR. bc_radiation_s )              &
     888                           .AND. j == nysv  )    .OR.                          &
     889                         ( ( bc_dirichlet_n .OR. bc_radiation_n )              &
     890                           .AND. j == nyn-1 ) )                                &
     891                THEN                                                           
     892                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 22 )     
     893!                                                                             
     894!--                Clear flag for WS1                                         
     895                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 )     
     896                ELSEIF ( BTEST(wall_flags_0(k,j+1,i),3)  .AND.                 &
     897                         BTEST(wall_flags_0(k,j+2,i),3)  .AND.                 &
     898                         BTEST(wall_flags_0(k,j-1,i),3) )                      &
     899                THEN
     900                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 23 )
     901!
     902!--                Clear flag for WS1
     903                   advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 )
     904                ENDIF
     905!
     906!--             w component - z-direction
     907!--             WS1 (24), WS3 (25), WS5 (26)
     908                flag_set = .FALSE.
     909                IF ( k == nzb+1 )  THEN
     910                   k_mm = nzb
     911                ELSE
     912                   k_mm = k - 2
     913                ENDIF
     914                IF ( k > nzt-1 )  THEN
     915                   k_pp = nzt+1
     916                ELSE
     917                   k_pp = k + 2
     918                ENDIF
     919                IF ( k > nzt-2 )  THEN
     920                   k_ppp = nzt+1
     921                ELSE
     922                   k_ppp = k + 3
     923                ENDIF 
     924               
     925                IF ( ( .NOT. BTEST(wall_flags_0(k-1,j,i),3)  .AND.             &
     926                       .NOT. BTEST(wall_flags_0(k,j,i),3)    .AND.             &
     927                             BTEST(wall_flags_0(k+1,j,i),3) )  .OR.            &
     928                     ( .NOT. BTEST(wall_flags_0(k-1,j,i),3)  .AND.             &
     929                             BTEST(wall_flags_0(k,j,i),3) )  .OR.              &
     930                     ( .NOT. BTEST(wall_flags_0(k+1,j,i),3)  .AND.             &
     931                             BTEST(wall_flags_0(k,j,i),3) )  .OR.              &       
     932                     k == nzt )                                                &
     933                THEN
     934!
     935!--                Please note, at k == nzb_w_inner(j,i) a flag is explictely
     936!--                set, although this is not a prognostic level. However,
     937!--                contrary to the advection of u,v and s this is necessary
     938!--                because flux_t(nzb_w_inner(j,i)) is used for the tendency
     939!--                at k == nzb_w_inner(j,i)+1.
     940                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 24 )
     941                   flag_set = .TRUE.
     942                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),3)     .OR.      &
     943                           .NOT. BTEST(wall_flags_0(k_ppp,j,i),3) ) .AND.      &
     944                                 BTEST(wall_flags_0(k-1,j,i),3)  .AND.         &
     945                                 BTEST(wall_flags_0(k,j,i),3)    .AND.         &
     946                                 BTEST(wall_flags_0(k+1,j,i),3)  .OR.          &
     947                         k == nzt - 1 )                                        &
     948                THEN                                                           
     949                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 25 )     
     950                   flag_set = .TRUE.                                           
     951                ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),3)  .AND.                &
     952                         BTEST(wall_flags_0(k-1,j,i),3)   .AND.                &
     953                         BTEST(wall_flags_0(k,j,i),3)     .AND.                &
     954                         BTEST(wall_flags_0(k+1,j,i),3)   .AND.                &
     955                         BTEST(wall_flags_0(k_pp,j,i),3)  .AND.                &
     956                         BTEST(wall_flags_0(k_ppp,j,i),3) .AND.                &
     957                         .NOT. flag_set )                                      &
     958                THEN
     959                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 26 )
     960                ENDIF
     961
    910962             ENDDO
    911963          ENDDO
    912           DO  i = nxl, nxr
    913              DO  j = nys, nyn
    914                 DO  k = nzb+1, nzt
    915 !
    916 !--                At first, set flags to WS1.
    917 !--                Since fluxes are swapped in advec_ws.f90, this is necessary to
    918 !--                in order to handle the left/south flux.
    919                    advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 27 )
    920                    advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 30 )
    921 !
    922 !--                w component - x-direction
    923 !--                WS1 (27), WS3 (28), WS5 (29)
    924                    IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),3) .OR.              &
    925                             ( (  bc_dirichlet_l .OR. bc_radiation_l )          &
    926                               .AND. i == nxl   )    .OR.                       &
    927                             ( (  bc_dirichlet_r .OR. bc_radiation_r )          &
    928                               .AND. i == nxr   ) )                             &
    929                    THEN
    930                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 27 )
    931                    ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),3)  .AND.      &
    932                                     BTEST(wall_flags_0(k,j,i+1),3)  .OR.       &
    933                               .NOT. BTEST(wall_flags_0(k,j,i-1),3) )           &
    934                                                         .OR.                   &
    935                             ( ( bc_dirichlet_r .OR. bc_radiation_r )           &
    936                               .AND. i == nxr-1 )    .OR.                       &
    937                             ( ( bc_dirichlet_l .OR.  bc_radiation_l )          &
    938                               .AND. i == nxlu  ) )                             &
    939                    THEN
    940                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 28 )
    941 !   
    942 !--                   Clear flag for WS1
    943                       advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 27 )
    944                    ELSEIF ( BTEST(wall_flags_0(k,j,i+1),3)  .AND.              &
    945                             BTEST(wall_flags_0(k,j,i+2),3)  .AND.              &
    946                             BTEST(wall_flags_0(k,j,i-1),3) )                   &
    947                    THEN
    948                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i),29 )
    949 !   
    950 !--                   Clear flag for WS1
    951                       advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 27 )
     964       ENDDO
     965!
     966!--    Exchange ghost points for advection flags
     967       CALL exchange_horiz_int( advc_flags_m, nys, nyn, nxl, nxr, nzt, nbgp )
     968!
     969!--    Set boundary flags at inflow and outflow boundary in case of
     970!--    non-cyclic boundary conditions.
     971       IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
     972          advc_flags_m(:,:,nxl-1) = advc_flags_m(:,:,nxl)
     973       ENDIF
     974
     975       IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
     976         advc_flags_m(:,:,nxr+1) = advc_flags_m(:,:,nxr)
     977       ENDIF
     978
     979       IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
     980          advc_flags_m(:,nyn+1,:) = advc_flags_m(:,nyn,:)
     981       ENDIF
     982
     983       IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
     984          advc_flags_m(:,nys-1,:) = advc_flags_m(:,nys,:)
     985       ENDIF
     986
     987    END SUBROUTINE ws_init_flags_momentum
     988
     989
     990!------------------------------------------------------------------------------!
     991! Description:
     992! ------------
     993!> Initialization of flags to control the order of the advection scheme near
     994!> solid walls and non-cyclic inflow boundaries, where the order is sucessively
     995!> degraded.
     996!------------------------------------------------------------------------------!
     997    SUBROUTINE ws_init_flags_scalar( non_cyclic_l, non_cyclic_n, non_cyclic_r, &
     998                                     non_cyclic_s, advc_flag, extensive_degrad )
     999
     1000
     1001       INTEGER(iwp) ::  i     !< index variable along x
     1002       INTEGER(iwp) ::  j     !< index variable along y
     1003       INTEGER(iwp) ::  k     !< index variable along z
     1004       INTEGER(iwp) ::  k_mm  !< dummy index along z
     1005       INTEGER(iwp) ::  k_pp  !< dummy index along z
     1006       INTEGER(iwp) ::  k_ppp !< dummy index along z
     1007       
     1008       INTEGER(iwp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::&
     1009                                                  advc_flag !< flag array to control order of scalar advection
     1010       
     1011       LOGICAL ::  flag_set     !< steering variable for advection flags
     1012       LOGICAL ::  non_cyclic_l !< flag that indicates non-cyclic boundary on the left
     1013       LOGICAL ::  non_cyclic_n !< flag that indicates non-cyclic boundary on the north
     1014       LOGICAL ::  non_cyclic_r !< flag that indicates non-cyclic boundary on the right
     1015       LOGICAL ::  non_cyclic_s !< flag that indicates non-cyclic boundary on the south
     1016       
     1017       LOGICAL, OPTIONAL ::  extensive_degrad !< flag indicating that extensive degradation is required, e.g. for
     1018                                              !< passive scalars nearby topography along the horizontal directions,
     1019                                              !< as no monotonic limiter can be applied there
     1020!
     1021!--    Set flags to steer the degradation of the advection scheme in advec_ws
     1022!--    near topography, inflow- and outflow boundaries as well as bottom and
     1023!--    top of model domain. advc_flags_m remains zero for all non-prognostic
     1024!--    grid points.
     1025       DO  i = nxl, nxr
     1026          DO  j = nys, nyn
     1027             DO  k = nzb+1, nzt
     1028                IF ( .NOT.  BTEST(wall_flags_0(k,j,i),0) )  CYCLE
     1029!
     1030!--             scalar - x-direction
     1031!--             WS1 (0), WS3 (1), WS5 (2)
     1032                IF ( ( .NOT. BTEST(wall_flags_0(k,j,i+1),0)                    &
     1033                 .OR.  .NOT. BTEST(wall_flags_0(k,j,i+2),0)                    &   
     1034                 .OR.  .NOT. BTEST(wall_flags_0(k,j,i-1),0) )                  &
     1035                   .OR.  ( non_cyclic_l  .AND.  i == 0  )                      &
     1036                   .OR.  ( non_cyclic_r  .AND.  i == nx ) )  THEN           
     1037                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )             
     1038                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+3),0)                &
     1039                    .AND.        BTEST(wall_flags_0(k,j,i+1),0)                &
     1040                    .AND.        BTEST(wall_flags_0(k,j,i+2),0)                &
     1041                    .AND.        BTEST(wall_flags_0(k,j,i-1),0)                &
     1042                         )                       .OR.                          &
     1043                         ( .NOT. BTEST(wall_flags_0(k,j,i-2),0)                &
     1044                    .AND.        BTEST(wall_flags_0(k,j,i+1),0)                &
     1045                    .AND.        BTEST(wall_flags_0(k,j,i+2),0)                &
     1046                    .AND.        BTEST(wall_flags_0(k,j,i-1),0)                &
     1047                         )                                                     &
     1048                                                 .OR.                          &
     1049                         ( non_cyclic_r  .AND.  i == nx-1 )  .OR.              &
     1050                         ( non_cyclic_l  .AND.  i == 1    ) )  THEN           
     1051                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )             
     1052                ELSEIF ( BTEST(wall_flags_0(k,j,i+1),0)                        &
     1053                   .AND. BTEST(wall_flags_0(k,j,i+2),0)                        &
     1054                   .AND. BTEST(wall_flags_0(k,j,i+3),0)                        &
     1055                   .AND. BTEST(wall_flags_0(k,j,i-1),0)                        &
     1056                   .AND. BTEST(wall_flags_0(k,j,i-2),0) )                      &
     1057                THEN                                                           
     1058                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 )             
     1059                ENDIF                                                         
     1060!                                                                             
     1061!--             scalar - y-direction                                           
     1062!--             WS1 (3), WS3 (4), WS5 (5)                                     
     1063                IF ( ( .NOT. BTEST(wall_flags_0(k,j+1,i),0)                    &
     1064                 .OR.  .NOT. BTEST(wall_flags_0(k,j+2,i),0)                    &   
     1065                 .OR.  .NOT. BTEST(wall_flags_0(k,j-1,i),0))                   &
     1066                   .OR.  ( non_cyclic_s  .AND.  j == 0  )                      &
     1067                   .OR.  ( non_cyclic_n  .AND.  j == ny ) )  THEN                                                           
     1068                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )             
     1069!                                                                             
     1070!--             WS3                                                           
     1071                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+3,i),0)                &
     1072                    .AND.        BTEST(wall_flags_0(k,j+1,i),0)                &
     1073                    .AND.        BTEST(wall_flags_0(k,j+2,i),0)                &
     1074                    .AND.        BTEST(wall_flags_0(k,j-1,i),0)                &
     1075                         )                       .OR.                          &
     1076                         ( .NOT. BTEST(wall_flags_0(k,j-2,i),0)                &
     1077                    .AND.        BTEST(wall_flags_0(k,j+1,i),0)                &
     1078                    .AND.        BTEST(wall_flags_0(k,j+2,i),0)                &
     1079                    .AND.        BTEST(wall_flags_0(k,j-1,i),0)                &
     1080                         )                                                     &
     1081                                                 .OR.                          &
     1082                         ( non_cyclic_s  .AND.  j == 1    )  .OR.              &
     1083                         ( non_cyclic_n  .AND.  j == ny-1 ) )  THEN           
     1084                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )             
     1085!                                                                             
     1086!--             WS5                                                           
     1087                ELSEIF ( BTEST(wall_flags_0(k,j+1,i),0)                        &
     1088                   .AND. BTEST(wall_flags_0(k,j+2,i),0)                        &
     1089                   .AND. BTEST(wall_flags_0(k,j+3,i),0)                        &
     1090                   .AND. BTEST(wall_flags_0(k,j-1,i),0)                        &
     1091                   .AND. BTEST(wall_flags_0(k,j-2,i),0) )                      &
     1092                THEN                                                           
     1093                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 )             
     1094                ENDIF
     1095!
     1096!--             Near topography, set horizontal advection scheme to 1st order
     1097!--             for passive scalars, even if only one direction may be
     1098!--             blocked by topography. These locations will be identified
     1099!--             by wall_flags_0 bit 31. Note, since several modules define
     1100!--             advection flags but may apply different scalar boundary
     1101!--             conditions, bit 31 is temporarily stored on advc_flags.
     1102!--             Moreover, note that this extended degradtion for passive
     1103!--             scalars is not required for the vertical direction as there
     1104!--             the monotonic limiter can be applied.
     1105                IF ( PRESENT( extensive_degrad ) )  THEN
     1106                   IF ( extensive_degrad )  THEN
     1107!
     1108!--                   At all grid points that are within a three-grid point
     1109!--                   range to topography, set 1st-order scheme.
     1110                      IF( BTEST( advc_flag(k,j,i), 31 ) )  THEN
     1111!
     1112!--                      Clear flags that might indicate higher-order
     1113!--                      advection along x- and y-direction.
     1114                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
     1115                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
     1116                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
     1117                         advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
     1118!
     1119!--                      Set flags that indicate 1st-order advection along
     1120!--                      x- and y-direction.
     1121                         advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
     1122                         advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) 
     1123                      ENDIF
     1124!
     1125!--                   Adjacent to this extended degradation zone, successively
     1126!--                   upgrade the order of the scheme if this grid point isn't
     1127!--                   flagged with bit 31 (indicating extended degradation
     1128!--                   zone).
     1129                      IF ( .NOT. BTEST( advc_flag(k,j,i), 31 ) )  THEN
     1130!
     1131!--                      x-direction. First, clear all previous settings, than
     1132!--                      set flag for 3rd-order scheme.
     1133                         IF ( BTEST( advc_flag(k,j,i-1), 31 )  .AND.        &
     1134                              BTEST( advc_flag(k,j,i+1), 31 ) )  THEN
     1135                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
     1136                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
     1137                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
     1138                           
     1139                            advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
     1140                         ENDIF
     1141!
     1142!--                      x-direction. First, clear all previous settings, than
     1143!--                      set flag for 5rd-order scheme.
     1144                         IF ( .NOT. BTEST( advc_flag(k,j,i-1), 31 )  .AND.  &
     1145                                    BTEST( advc_flag(k,j,i-2), 31 )  .AND.  &
     1146                              .NOT. BTEST( advc_flag(k,j,i+1), 31 )  .AND.  &
     1147                                    BTEST( advc_flag(k,j,i+2), 31 ) )  THEN
     1148                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
     1149                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
     1150                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
     1151                           
     1152                            advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 )
     1153                         ENDIF
     1154!
     1155!--                      y-direction. First, clear all previous settings, than
     1156!--                      set flag for 3rd-order scheme.
     1157                         IF ( BTEST( advc_flag(k,j-1,i), 31 )  .AND.        &
     1158                              BTEST( advc_flag(k,j+1,i), 31 ) )  THEN
     1159                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
     1160                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
     1161                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
     1162                           
     1163                            advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
     1164                         ENDIF
     1165!
     1166!--                      y-direction. First, clear all previous settings, than
     1167!--                      set flag for 5rd-order scheme.
     1168                         IF ( .NOT. BTEST( advc_flag(k,j-1,i), 31 )  .AND.  &
     1169                                    BTEST( advc_flag(k,j-2,i), 31 )  .AND.  &
     1170                              .NOT. BTEST( advc_flag(k,j+1,i), 31 )  .AND.  &
     1171                                    BTEST( advc_flag(k,j+2,i), 31 ) )  THEN
     1172                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
     1173                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
     1174                            advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
     1175                           
     1176                            advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 )
     1177                         ENDIF
     1178                      ENDIF
     1179                   
    9521180                   ENDIF
    953 !
    954 !--                w component - y-direction
    955 !--                WS1 (30), WS3 (31), WS5 (32)
    956                    IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),3) .OR.              &
    957                             ( ( bc_dirichlet_s .OR. bc_radiation_s )           &
    958                               .AND. j == nys   )    .OR.                       &
    959                             ( ( bc_dirichlet_n .OR. bc_radiation_n )           &
    960                               .AND. j == nyn   ) )                             &
    961                    THEN
    962                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 30 )
    963                    ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),3)  .AND.      &
    964                                     BTEST(wall_flags_0(k,j+1,i),3)  .OR.       &
    965                               .NOT. BTEST(wall_flags_0(k,j-1,i),3) )           &
    966                                                         .OR.                   &
    967                             ( ( bc_dirichlet_s .OR. bc_radiation_s )           &
    968                               .AND. j == nysv  )    .OR.                       &
    969                             ( ( bc_dirichlet_n .OR. bc_radiation_n )           &
    970                               .AND. j == nyn-1 ) )                             &
    971                    THEN
    972                       advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 31 )
    973 !
    974 !--                   Clear flag for WS1
    975                       advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 30 )
    976                    ELSEIF ( BTEST(wall_flags_0(k,j+1,i),3)  .AND.              &
    977                             BTEST(wall_flags_0(k,j+2,i),3)  .AND.              &
    978                             BTEST(wall_flags_0(k,j-1,i),3) )                   &
    979                    THEN
    980                       advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 0 )
    981 !
    982 !--                   Clear flag for WS1
    983                       advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 30 )
     1181                   
     1182!
     1183!--                Near lateral boundary flags might be overwritten. Set
     1184!--                them again.
     1185!--                x-direction
     1186                   IF ( ( non_cyclic_l  .AND.  i == 0  )  .OR.                 &
     1187                        ( non_cyclic_r  .AND.  i == nx ) )  THEN
     1188                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
     1189                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
     1190                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
     1191                         
     1192                      advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 )
    9841193                   ENDIF
    985 !
    986 !--                w component - z-direction
    987 !--                WS1 (33), WS3 (34), WS5 (35)
    988                    flag_set = .FALSE.
    989                    IF ( k == nzb+1 )  THEN
    990                       k_mm = nzb
    991                    ELSE
    992                       k_mm = k - 2
     1194                   
     1195                   IF ( ( non_cyclic_l  .AND.  i == 1    )  .OR.               &
     1196                        ( non_cyclic_r  .AND.  i == nx-1 ) )  THEN
     1197                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 )
     1198                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 )
     1199                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 )
     1200                        
     1201                      advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 )
    9931202                   ENDIF
    994                    IF ( k > nzt-1 )  THEN
    995                       k_pp = nzt+1
    996                    ELSE
    997                       k_pp = k + 2
     1203!
     1204!--                y-direction
     1205                   IF ( ( non_cyclic_n  .AND.  j == 0  )  .OR.                 &
     1206                        ( non_cyclic_s  .AND.  j == ny ) )  THEN
     1207                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
     1208                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
     1209                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
     1210                         
     1211                      advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 )
    9981212                   ENDIF
    999                    IF ( k > nzt-2 )  THEN
    1000                       k_ppp = nzt+1
    1001                    ELSE
    1002                       k_ppp = k + 3
    1003                    ENDIF 
    10041213                   
    1005                    IF ( ( .NOT. BTEST(wall_flags_0(k-1,j,i),3)  .AND.          &
    1006                           .NOT. BTEST(wall_flags_0(k,j,i),3)    .AND.          &
    1007                                 BTEST(wall_flags_0(k+1,j,i),3) )  .OR.         &
    1008                         ( .NOT. BTEST(wall_flags_0(k-1,j,i),3)  .AND.          &
    1009                                 BTEST(wall_flags_0(k,j,i),3) )  .OR.           &
    1010                         ( .NOT. BTEST(wall_flags_0(k+1,j,i),3)  .AND.          &
    1011                                 BTEST(wall_flags_0(k,j,i),3) )  .OR.           &       
    1012                         k == nzt )                                             &
    1013                    THEN
    1014 !
    1015 !--                   Please note, at k == nzb_w_inner(j,i) a flag is explictely
    1016 !--                   set, although this is not a prognostic level. However,
    1017 !--                   contrary to the advection of u,v and s this is necessary
    1018 !--                   because flux_t(nzb_w_inner(j,i)) is used for the tendency
    1019 !--                   at k == nzb_w_inner(j,i)+1.
    1020                       advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 1 )
    1021                       flag_set = .TRUE.
    1022                    ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),3)     .OR.   &
    1023                               .NOT. BTEST(wall_flags_0(k_ppp,j,i),3) ) .AND.   &
    1024                                     BTEST(wall_flags_0(k-1,j,i),3)  .AND.      &
    1025                                     BTEST(wall_flags_0(k,j,i),3)    .AND.      &
    1026                                     BTEST(wall_flags_0(k+1,j,i),3)  .OR.       &
    1027                             k == nzt - 1 )                                     &
    1028                    THEN
    1029                       advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 2 )
    1030                       flag_set = .TRUE.
    1031                    ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),3)  .AND.             &
    1032                             BTEST(wall_flags_0(k-1,j,i),3)   .AND.             &
    1033                             BTEST(wall_flags_0(k,j,i),3)     .AND.             &
    1034                             BTEST(wall_flags_0(k+1,j,i),3)   .AND.             &
    1035                             BTEST(wall_flags_0(k_pp,j,i),3)  .AND.             &
    1036                             BTEST(wall_flags_0(k_ppp,j,i),3) .AND.             &
    1037                             .NOT. flag_set )                                   &
    1038                    THEN
    1039                       advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 3 )
     1214                   IF ( ( non_cyclic_n  .AND.  j == 1    )  .OR.               &
     1215                        ( non_cyclic_s  .AND.  j == ny-1 ) )  THEN
     1216                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 )
     1217                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 )
     1218                      advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 )
     1219                         
     1220                      advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 )
    10401221                   ENDIF
    1041 
    1042                 ENDDO
     1222                   
     1223                ENDIF
     1224               
     1225               
     1226!                                                                             
     1227!--             scalar - z-direction                                           
     1228!--             WS1 (6), WS3 (7), WS5 (8)                                     
     1229                IF ( k == nzb+1 )  THEN                                       
     1230                   k_mm = nzb                                                 
     1231                ELSE                                                           
     1232                   k_mm = k - 2                                               
     1233                ENDIF                                                         
     1234                IF ( k > nzt-1 )  THEN                                         
     1235                   k_pp = nzt+1                                               
     1236                ELSE                                                           
     1237                   k_pp = k + 2                                               
     1238                ENDIF                                                         
     1239                IF ( k > nzt-2 )  THEN                                         
     1240                   k_ppp = nzt+1                                               
     1241                ELSE                                                           
     1242                   k_ppp = k + 3                                               
     1243                ENDIF                                                         
     1244                                                                               
     1245                flag_set = .FALSE.                                             
     1246                IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),0)  .AND.               &
     1247                           BTEST(wall_flags_0(k,j,i),0)    .OR.                &
     1248                     .NOT. BTEST(wall_flags_0(k_pp,j,i),0) .AND.               &                             
     1249                           BTEST(wall_flags_0(k,j,i),0)    .OR.                &
     1250                     k == nzt )                                                &
     1251                THEN                                                           
     1252                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 6 )             
     1253                   flag_set = .TRUE.                                           
     1254                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),0)    .OR.       &
     1255                           .NOT. BTEST(wall_flags_0(k_ppp,j,i),0) ) .AND.      &
     1256                               BTEST(wall_flags_0(k-1,j,i),0)  .AND.           &
     1257                               BTEST(wall_flags_0(k,j,i),0)    .AND.           &
     1258                               BTEST(wall_flags_0(k+1,j,i),0)  .AND.           &
     1259                               BTEST(wall_flags_0(k_pp,j,i),0) .OR.            &   
     1260                         k == nzt - 1 )                                        &
     1261                THEN                                                           
     1262                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 7 )             
     1263                   flag_set = .TRUE.                                           
     1264                ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),0)                       &
     1265                  .AND.  BTEST(wall_flags_0(k-1,j,i),0)                        &
     1266                  .AND.  BTEST(wall_flags_0(k,j,i),0)                          &
     1267                  .AND.  BTEST(wall_flags_0(k+1,j,i),0)                        &
     1268                  .AND.  BTEST(wall_flags_0(k_pp,j,i),0)                       &
     1269                  .AND.  BTEST(wall_flags_0(k_ppp,j,i),0)                      &
     1270                  .AND. .NOT. flag_set )                                       &
     1271                THEN
     1272                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 8 )
     1273                ENDIF
     1274
    10431275             ENDDO
    10441276          ENDDO
    1045 
     1277       ENDDO
     1278!
     1279!--    Exchange 3D integer wall_flags.
     1280!
     1281!--    Exchange ghost points for advection flags
     1282       CALL exchange_horiz_int( advc_flag, nys, nyn, nxl, nxr, nzt, nbgp )
     1283!
     1284!--    Set boundary flags at inflow and outflow boundary in case of
     1285!--    non-cyclic boundary conditions.
     1286       IF ( non_cyclic_l )  THEN
     1287          advc_flag(:,:,nxl-1) = advc_flag(:,:,nxl)
    10461288       ENDIF
    10471289
    1048 
    1049 !
    1050 !--    Exchange 3D integer wall_flags.
    1051        IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme'     &
    1052           )  THEN 
    1053 !
    1054 !--       Exchange ghost points for advection flags
    1055           CALL exchange_horiz_int( advc_flags_1, nys, nyn, nxl, nxr, nzt, nbgp )
    1056           CALL exchange_horiz_int( advc_flags_2, nys, nyn, nxl, nxr, nzt, nbgp )
    1057 !
    1058 !--       Set boundary flags at inflow and outflow boundary in case of
    1059 !--       non-cyclic boundary conditions.
    1060           IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
    1061              advc_flags_1(:,:,nxl-1) = advc_flags_1(:,:,nxl)
    1062              advc_flags_2(:,:,nxl-1) = advc_flags_2(:,:,nxl)
    1063           ENDIF
    1064 
    1065           IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
    1066             advc_flags_1(:,:,nxr+1) = advc_flags_1(:,:,nxr)
    1067             advc_flags_2(:,:,nxr+1) = advc_flags_2(:,:,nxr)
    1068           ENDIF
    1069 
    1070           IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
    1071              advc_flags_1(:,nyn+1,:) = advc_flags_1(:,nyn,:)
    1072              advc_flags_2(:,nyn+1,:) = advc_flags_2(:,nyn,:)
    1073           ENDIF
    1074 
    1075           IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
    1076              advc_flags_1(:,nys-1,:) = advc_flags_1(:,nys,:)
    1077              advc_flags_2(:,nys-1,:) = advc_flags_2(:,nys,:)
    1078           ENDIF
     1290       IF ( non_cyclic_r )  THEN
     1291         advc_flag(:,:,nxr+1) = advc_flag(:,:,nxr)
     1292       ENDIF
     1293
     1294       IF ( non_cyclic_n )  THEN
     1295          advc_flag(:,nyn+1,:) = advc_flag(:,nyn,:)
     1296       ENDIF
     1297
     1298       IF ( non_cyclic_s )  THEN
     1299          advc_flag(:,nys-1,:) = advc_flag(:,nys,:)
     1300       ENDIF
    10791301 
    1080        ENDIF
    1081 
    1082 
    1083     END SUBROUTINE ws_init_flags
    1084 
    1085 
     1302
     1303
     1304    END SUBROUTINE ws_init_flags_scalar   
     1305   
    10861306!------------------------------------------------------------------------------!
    10871307! Description:
     
    11231343!> Scalar advection - Call for grid point i,j
    11241344!------------------------------------------------------------------------------!
    1125     SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local,            &
     1345    SUBROUTINE advec_s_ws_ij( advc_flag, i, j, sk, sk_char, swap_flux_y_local, &
    11261346                              swap_diss_y_local, swap_flux_x_local,            &
    1127                               swap_diss_x_local, i_omp, tn, flux_limitation )
     1347                              swap_diss_x_local, i_omp, tn,                    &
     1348                              non_cyclic_l, non_cyclic_n,                      &
     1349                              non_cyclic_r, non_cyclic_s,                      &
     1350                              flux_limitation )
    11281351
    11291352
     
    11401363       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 
    11411364       INTEGER(iwp) ::  tn        !< number of OpenMP thread
    1142 
     1365       
     1366       INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::   &
     1367                                                  advc_flag !< flag array to control order of scalar advection
     1368
     1369       LOGICAL           ::  non_cyclic_l    !< flag that indicates non-cyclic boundary on the left
     1370       LOGICAL           ::  non_cyclic_n    !< flag that indicates non-cyclic boundary on the north
     1371       LOGICAL           ::  non_cyclic_r    !< flag that indicates non-cyclic boundary on the right
     1372       LOGICAL           ::  non_cyclic_s    !< flag that indicates non-cyclic boundary on the south                                           
    11431373       LOGICAL, OPTIONAL ::  flux_limitation !< flag indicating flux limitation of the vertical advection
    11441374       LOGICAL           ::  limiter         !< control flag indicating the application of flux limitation
     
    11961426!--    instead of the entire subdomain. This should lead to better
    11971427!--    load balance between boundary and non-boundary PEs.
    1198        IF( ( bc_dirichlet_l  .OR.  bc_radiation_l )  .AND.  i <= nxl + 2  .OR. &
    1199            ( bc_dirichlet_r  .OR.  bc_radiation_r )  .AND.  i >= nxr - 2  .OR. &
    1200            ( bc_dirichlet_s  .OR.  bc_radiation_s )  .AND.  j <= nys + 2  .OR. &
    1201            ( bc_dirichlet_n  .OR.  bc_radiation_n )  .AND.  j >= nyn - 2 )  THEN
     1428       IF( non_cyclic_l  .AND.  i <= nxl + 2  .OR.                            &
     1429           non_cyclic_r  .AND.  i >= nxr - 2  .OR.                            &
     1430           non_cyclic_s  .AND.  j <= nys + 2  .OR.                            &
     1431           non_cyclic_n  .AND.  j >= nyn - 2 )  THEN
    12021432          nzb_max_l = nzt
    12031433       ELSE
     
    12151445          DO  k = nzb+1, nzb_max_l
    12161446
    1217              ibit5 = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )
    1218              ibit4 = REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )
    1219              ibit3 = REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )
     1447             ibit5 = REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp )
     1448             ibit4 = REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp )
     1449             ibit3 = REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp )
    12201450
    12211451             v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
     
    12761506          DO  k = nzb+1, nzb_max_l
    12771507
    1278              ibit2 = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )
    1279              ibit1 = REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )
    1280              ibit0 = REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )
     1508             ibit2 = REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp )
     1509             ibit1 = REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp )
     1510             ibit0 = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )
    12811511
    12821512             u_comp                     = u(k,j,i) - u_gtrans + u_stokes_zu(k)
     
    13401570!--       flux at the end.
    13411571
    1342           ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp )
    1343           ibit1 = REAL( IBITS(advc_flags_1(k,j,i),1,1), KIND = wp )
    1344           ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp )
     1572          ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
     1573          ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp )
     1574          ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp )
    13451575
    13461576          u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
     
    13751605                                       )
    13761606
    1377           ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp )
    1378           ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp )
    1379           ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp )
     1607          ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp )
     1608          ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp )
     1609          ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp )
    13801610
    13811611          v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
     
    14451675!--    calculated explicetely for the tendency at
    14461676!--    the first w-level. For topography wall this is done implicitely by
    1447 !--    advc_flags_1.
     1677!--    advc_flag.
    14481678       flux_t(nzb) = 0.0_wp
    14491679       diss_t(nzb) = 0.0_wp
    14501680       DO  k = nzb+1, nzb+2
    1451           ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
    1452           ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
    1453           ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
     1681          ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     1682          ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     1683          ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
    14541684!
    14551685!--       k index has to be modified near bottom and top, else array
     
    14911721       
    14921722       DO  k = nzb+3, nzt-2
    1493           ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
    1494           ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
    1495           ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
     1723          ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     1724          ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     1725          ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
    14961726
    14971727          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
     
    15261756       
    15271757       DO  k = nzt-1, nzt
    1528           ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
    1529           ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
    1530           ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
     1758          ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     1759          ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     1760          ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
    15311761!
    15321762!--       k index has to be modified near bottom and top, else array
     
    16751905          diss_d    = diss_t(k-1)
    16761906         
    1677           ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp )
    1678           ibit1 = REAL( IBITS(advc_flags_1(k,j,i),1,1), KIND = wp )
    1679           ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp )
     1907          ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
     1908          ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp )
     1909          ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp )
    16801910         
    1681           ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp )
    1682           ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp )
    1683           ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp )
     1911          ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp )
     1912          ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp )
     1913          ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp )
    16841914         
    1685           ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
    1686           ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
    1687           ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
     1915          ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     1916          ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     1917          ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
    16881918!
    16891919!--       Calculate the divergence of the velocity field. A respective
     
    16921922          div         =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )            &
    16931923                          - u(k,j,i)   * (                                    &
    1694                         REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )   &
    1695                       + REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )   &
    1696                       + REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )   &
     1924                        REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )      &
     1925                      + REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp )      &
     1926                      + REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp )      &
    16971927                                         )                                    &
    16981928                          ) * ddx                                             &
    16991929                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )            &
    17001930                          - v(k,j,i)   * (                                    &
    1701                         REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )   &
    1702                       + REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )   &
    1703                       + REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )   &
     1931                        REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp )      &
     1932                      + REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp )      &
     1933                      + REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp )      &
    17041934                                         )                                    &
    17051935                          ) * ddy                                             &
     
    17081938                          - w(k-1,j,i) * rho_air_zw(k-1) *                    &
    17091939                                         (                                    &
    1710                         REAL( IBITS(advc_flags_1(k-1,j,i),6,1), KIND = wp )   &
    1711                       + REAL( IBITS(advc_flags_1(k-1,j,i),7,1), KIND = wp )   &
    1712                       + REAL( IBITS(advc_flags_1(k-1,j,i),8,1), KIND = wp )   &
     1940                        REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp )      &
     1941                      + REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp )      &
     1942                      + REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp )      &
    17131943                                         )                                    &     
    17141944                          ) * drho_air(k) * ddzw(k)
     
    18972127       INTEGER(iwp) ::  tn        !< number of OpenMP thread
    18982128       
    1899        REAL(wp)    ::  ibit9    !< flag indicating 1st-order scheme along x-direction
    1900        REAL(wp)    ::  ibit10   !< flag indicating 3rd-order scheme along x-direction
    1901        REAL(wp)    ::  ibit11   !< flag indicating 5th-order scheme along x-direction
    1902        REAL(wp)    ::  ibit12   !< flag indicating 1st-order scheme along y-direction
    1903        REAL(wp)    ::  ibit13   !< flag indicating 3rd-order scheme along y-direction
    1904        REAL(wp)    ::  ibit14   !< flag indicating 5th-order scheme along y-direction
    1905        REAL(wp)    ::  ibit15   !< flag indicating 1st-order scheme along z-direction
    1906        REAL(wp)    ::  ibit16   !< flag indicating 3rd-order scheme along z-direction
    1907        REAL(wp)    ::  ibit17   !< flag indicating 5th-order scheme along z-direction
     2129       REAL(wp)    ::  ibit0   !< flag indicating 1st-order scheme along x-direction
     2130       REAL(wp)    ::  ibit1   !< flag indicating 3rd-order scheme along x-direction
     2131       REAL(wp)    ::  ibit2   !< flag indicating 5th-order scheme along x-direction
     2132       REAL(wp)    ::  ibit3   !< flag indicating 1st-order scheme along y-direction
     2133       REAL(wp)    ::  ibit4   !< flag indicating 3rd-order scheme along y-direction
     2134       REAL(wp)    ::  ibit5   !< flag indicating 5th-order scheme along y-direction
     2135       REAL(wp)    ::  ibit6   !< flag indicating 1st-order scheme along z-direction
     2136       REAL(wp)    ::  ibit7   !< flag indicating 3rd-order scheme along z-direction
     2137       REAL(wp)    ::  ibit8   !< flag indicating 5th-order scheme along z-direction
    19082138       REAL(wp)    ::  diss_d   !< artificial dissipation term at grid box bottom
    19092139       REAL(wp)    ::  div      !< diverence on u-grid
     
    19442174          DO  k = nzb+1, nzb_max_l
    19452175
    1946              ibit14 = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )
    1947              ibit13 = REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp )
    1948              ibit12 = REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp )
     2176             ibit5 = REAL( IBITS(advc_flags_m(k,j-1,i),5,1), KIND = wp )
     2177             ibit4 = REAL( IBITS(advc_flags_m(k,j-1,i),4,1), KIND = wp )
     2178             ibit3 = REAL( IBITS(advc_flags_m(k,j-1,i),3,1), KIND = wp )
    19492179
    19502180             v_comp(k)      = v(k,j,i) + v(k,j,i-1) - gv
    1951              flux_s_u(k,tn) = v_comp(k) * (                                   &
    1952                             ( 37.0_wp * ibit14 * adv_mom_5                    &
    1953                          +     7.0_wp * ibit13 * adv_mom_3                    &
    1954                          +              ibit12 * adv_mom_1                    &
    1955                             ) *                                               &
    1956                                         ( u(k,j,i)   + u(k,j-1,i) )           &
    1957                      -      (  8.0_wp * ibit14 * adv_mom_5                    &
    1958                          +              ibit13 * adv_mom_3                    &
    1959                             ) *                                               &
    1960                                         ( u(k,j+1,i) + u(k,j-2,i) )           &
    1961                      +      (           ibit14 * adv_mom_5                    &
    1962                             ) *                                               &
    1963                                         ( u(k,j+2,i) + u(k,j-3,i) )           &
     2181             flux_s_u(k,tn) = v_comp(k) * (                                    &
     2182                            ( 37.0_wp * ibit5 * adv_mom_5                      &
     2183                         +     7.0_wp * ibit4 * adv_mom_3                      &
     2184                         +              ibit3 * adv_mom_1                      &
     2185                            ) *                                                &
     2186                                        ( u(k,j,i)   + u(k,j-1,i) )            &
     2187                     -      (  8.0_wp * ibit5 * adv_mom_5                      &
     2188                         +              ibit4 * adv_mom_3                      &
     2189                            ) *                                                &
     2190                                        ( u(k,j+1,i) + u(k,j-2,i) )            &
     2191                     +      (           ibit5 * adv_mom_5                      &
     2192                            ) *                                                &
     2193                                        ( u(k,j+2,i) + u(k,j-3,i) )            &
    19642194                                          )
    19652195
    1966              diss_s_u(k,tn) = - ABS ( v_comp(k) ) * (                         &
    1967                             ( 10.0_wp * ibit14 * adv_mom_5                    &
    1968                          +     3.0_wp * ibit13 * adv_mom_3                    &
    1969                          +              ibit12 * adv_mom_1                    &
    1970                             ) *                                               &
    1971                                         ( u(k,j,i)   - u(k,j-1,i) )           &
    1972                      -      (  5.0_wp * ibit14 * adv_mom_5                    &
    1973                          +              ibit13 * adv_mom_3                    &
    1974                             ) *                                               &
    1975                                         ( u(k,j+1,i) - u(k,j-2,i) )           &
    1976                      +      (           ibit14 * adv_mom_5                    &
    1977                             ) *                                               &
    1978                                         ( u(k,j+2,i) - u(k,j-3,i) )           &
     2196             diss_s_u(k,tn) = - ABS ( v_comp(k) ) * (                          &
     2197                            ( 10.0_wp * ibit5 * adv_mom_5                      &
     2198                         +     3.0_wp * ibit4 * adv_mom_3                      &
     2199                         +              ibit3 * adv_mom_1                      &
     2200                            ) *                                                &
     2201                                        ( u(k,j,i)   - u(k,j-1,i) )            &
     2202                     -      (  5.0_wp * ibit5 * adv_mom_5                      &
     2203                         +              ibit4 * adv_mom_3                      &
     2204                            ) *                                                &
     2205                                        ( u(k,j+1,i) - u(k,j-2,i) )            &
     2206                     +      (           ibit5 * adv_mom_5                      &
     2207                            ) *                                                &
     2208                                        ( u(k,j+2,i) - u(k,j-3,i) )            &
    19792209                                                    )
    19802210
     
    19842214
    19852215             v_comp(k)      = v(k,j,i) + v(k,j,i-1) - gv
    1986              flux_s_u(k,tn) = v_comp(k) * (                                   &
    1987                            37.0_wp * ( u(k,j,i) + u(k,j-1,i)   )              &
    1988                          -  8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) )              &
    1989                          +           ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
    1990              diss_s_u(k,tn) = - ABS(v_comp(k)) * (                            &
    1991                            10.0_wp * ( u(k,j,i) - u(k,j-1,i)   )              &
    1992                          -  5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) )              &
     2216             flux_s_u(k,tn) = v_comp(k) * (                                    &
     2217                           37.0_wp * ( u(k,j,i)   + u(k,j-1,i)   )             &
     2218                         -  8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) )               &
     2219                         +           ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 
     2220             diss_s_u(k,tn) = - ABS(v_comp(k)) * (                             &
     2221                           10.0_wp * ( u(k,j,i)   - u(k,j-1,i)   )             &
     2222                         -  5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) )               &
    19932223                         +           ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
    19942224
     
    20022232          DO  k = nzb+1, nzb_max_l
    20032233
    2004              ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )
    2005              ibit10 = REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )
    2006              ibit9  = REAL( IBITS(advc_flags_1(k,j,i-1),9,1),  KIND = wp )
    2007 
    2008              u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
    2009              flux_l_u(k,j,tn) = u_comp_l * (                                  &
    2010                               ( 37.0_wp * ibit11 * adv_mom_5                  &
    2011                            +     7.0_wp * ibit10 * adv_mom_3                  &
    2012                            +              ibit9  * adv_mom_1                  &
    2013                               ) *                                             &
    2014                                           ( u(k,j,i)   + u(k,j,i-1) )         &
    2015                        -      (  8.0_wp * ibit11 * adv_mom_5                  &
    2016                            +              ibit10 * adv_mom_3                  &
    2017                               ) *                                             &
    2018                                           ( u(k,j,i+1) + u(k,j,i-2) )         &
    2019                        +      (           ibit11 * adv_mom_5                  &
    2020                               ) *                                             &
    2021                                           ( u(k,j,i+2) + u(k,j,i-3) )         &
    2022                                            )
    2023 
    2024               diss_l_u(k,j,tn) = - ABS( u_comp_l ) * (                        &
    2025                               ( 10.0_wp * ibit11 * adv_mom_5                  &
    2026                            +     3.0_wp * ibit10 * adv_mom_3                  &
    2027                            +              ibit9  * adv_mom_1                  &
    2028                               ) *                                             &
    2029                                         ( u(k,j,i)   - u(k,j,i-1) )           &
    2030                        -      (  5.0_wp * ibit11 * adv_mom_5                  &
    2031                            +              ibit10 * adv_mom_3                  &
    2032                               ) *                                             &
    2033                                         ( u(k,j,i+1) - u(k,j,i-2) )           &
    2034                        +      (           ibit11 * adv_mom_5                  &
    2035                               ) *                                             &
    2036                                         ( u(k,j,i+2) - u(k,j,i-3) )           &
    2037                                                      )
    2038 
    2039           ENDDO
    2040 
    2041           DO  k = nzb_max_l+1, nzt
     2234             ibit2 = REAL( IBITS(advc_flags_m(k,j,i-1),2,1), KIND = wp )
     2235             ibit1 = REAL( IBITS(advc_flags_m(k,j,i-1),1,1), KIND = wp )
     2236             ibit0 = REAL( IBITS(advc_flags_m(k,j,i-1),0,1), KIND = wp )
    20422237
    20432238             u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
    20442239             flux_l_u(k,j,tn) = u_comp_l * (                                   &
    2045                              37.0_wp * ( u(k,j,i) + u(k,j,i-1)   )             &
     2240                              ( 37.0_wp * ibit2 * adv_mom_5                    &
     2241                           +     7.0_wp * ibit1 * adv_mom_3                    &
     2242                           +              ibit0 * adv_mom_1                    &
     2243                              ) *                                              &
     2244                                          ( u(k,j,i)   + u(k,j,i-1) )          &
     2245                       -      (  8.0_wp * ibit2 * adv_mom_5                    &
     2246                           +              ibit1 * adv_mom_3                    &
     2247                              ) *                                              &
     2248                                          ( u(k,j,i+1) + u(k,j,i-2) )          &
     2249                       +      (           ibit2 * adv_mom_5                    &
     2250                              ) *                                              &
     2251                                          ( u(k,j,i+2) + u(k,j,i-3) )          &
     2252                                           )                                 
     2253                                                                             
     2254              diss_l_u(k,j,tn) = - ABS( u_comp_l ) * (                         &
     2255                              ( 10.0_wp * ibit2 * adv_mom_5                    &
     2256                           +     3.0_wp * ibit1 * adv_mom_3                    &
     2257                           +              ibit0  * adv_mom_1                   &
     2258                              ) *                                              &
     2259                                        ( u(k,j,i)   - u(k,j,i-1) )            &
     2260                       -      (  5.0_wp * ibit2 * adv_mom_5                    &
     2261                           +              ibit1 * adv_mom_3                    &
     2262                              ) *                                              &
     2263                                        ( u(k,j,i+1) - u(k,j,i-2) )            &
     2264                       +      (           ibit2 * adv_mom_5                    &
     2265                              ) *                                              &
     2266                                        ( u(k,j,i+2) - u(k,j,i-3) )            &
     2267                                                     )
     2268
     2269          ENDDO
     2270
     2271          DO  k = nzb_max_l+1, nzt
     2272
     2273             u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
     2274             flux_l_u(k,j,tn) = u_comp_l * (                                   &
     2275                             37.0_wp * ( u(k,j,i)   + u(k,j,i-1)   )           &
    20462276                           -  8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) )             &
    20472277                           +           ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
    20482278             diss_l_u(k,j,tn) = - ABS(u_comp_l) * (                            &
    2049                              10.0_wp * ( u(k,j,i) - u(k,j,i-1)   )             &
     2279                             10.0_wp * ( u(k,j,i)   - u(k,j,i-1)   )           &
    20502280                           -  5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) )             &
    20512281                           +           ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
     
    20592289       DO  k = nzb+1, nzb_max_l
    20602290
    2061           ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp )
    2062           ibit10 = REAL( IBITS(advc_flags_1(k,j,i),10,1), KIND = wp )
    2063           ibit9  = REAL( IBITS(advc_flags_1(k,j,i),9,1), KIND = wp )
     2291          ibit2 = REAL( IBITS(advc_flags_m(k,j,i),2,1), KIND = wp )
     2292          ibit1 = REAL( IBITS(advc_flags_m(k,j,i),1,1), KIND = wp )
     2293          ibit0 = REAL( IBITS(advc_flags_m(k,j,i),0,1), KIND = wp )
    20642294
    20652295          u_comp(k) = u(k,j,i+1) + u(k,j,i)
    2066           flux_r(k) = ( u_comp(k) - gu ) * (                                  &
    2067                      ( 37.0_wp * ibit11 * adv_mom_5                           &
    2068                   +     7.0_wp * ibit10 * adv_mom_3                           &
    2069                   +              ibit9  * adv_mom_1                           &
    2070                      ) *                                                      &
    2071                                     ( u(k,j,i+1) + u(k,j,i)   )               &
    2072               -      (  8.0_wp * ibit11 * adv_mom_5                           &
    2073                   +              ibit10 * adv_mom_3                           &
    2074                      ) *                                                      &
    2075                                     ( u(k,j,i+2) + u(k,j,i-1) )               &
    2076               +      (           ibit11 * adv_mom_5                           &
    2077                      ) *                                                      &
    2078                                     ( u(k,j,i+3) + u(k,j,i-2) )               &
    2079                                            )
    2080 
    2081           diss_r(k) = - ABS( u_comp(k) - gu ) * (                             &
    2082                      ( 10.0_wp * ibit11 * adv_mom_5                           &
    2083                   +     3.0_wp * ibit10 * adv_mom_3                           &
    2084                   +              ibit9  * adv_mom_1                           &
    2085                      ) *                                                      &
    2086                                     ( u(k,j,i+1) - u(k,j,i)   )               &
    2087               -      (  5.0_wp * ibit11 * adv_mom_5                           &
    2088                   +              ibit10 * adv_mom_3                           &
    2089                      ) *                                                      &
    2090                                     ( u(k,j,i+2) - u(k,j,i-1) )               &
    2091               +      (           ibit11 * adv_mom_5                           &
    2092                      ) *                                                      &
    2093                                     ( u(k,j,i+3) - u(k,j,i-2) )               &
     2296          flux_r(k) = ( u_comp(k) - gu ) * (                                   &
     2297                     ( 37.0_wp * ibit2 * adv_mom_5                             &
     2298                  +     7.0_wp * ibit1 * adv_mom_3                             &
     2299                  +              ibit0 * adv_mom_1                             &
     2300                     ) *                                                       &
     2301                                    ( u(k,j,i+1) + u(k,j,i)   )                &
     2302              -      (  8.0_wp * ibit2 * adv_mom_5                             &
     2303                  +              ibit1 * adv_mom_3                             &
     2304                     ) *                                                       &
     2305                                    ( u(k,j,i+2) + u(k,j,i-1) )                &
     2306              +      (           ibit2 * adv_mom_5                             &
     2307                     ) *                                                       &
     2308                                    ( u(k,j,i+3) + u(k,j,i-2) )                &
     2309                                           )                                 
     2310                                                                             
     2311          diss_r(k) = - ABS( u_comp(k) - gu ) * (                              &
     2312                     ( 10.0_wp * ibit2 * adv_mom_5                             &
     2313                  +     3.0_wp * ibit1 * adv_mom_3                             &
     2314                  +              ibit0 * adv_mom_1                             &
     2315                     ) *                                                       &
     2316                                    ( u(k,j,i+1) - u(k,j,i)   )                &
     2317              -      (  5.0_wp * ibit2 * adv_mom_5                             &
     2318                  +              ibit1 * adv_mom_3                             &
     2319                     ) *                                                       &
     2320                                    ( u(k,j,i+2) - u(k,j,i-1) )                &
     2321              +      (           ibit2 * adv_mom_5                             &
     2322                     ) *                                                       &
     2323                                    ( u(k,j,i+3) - u(k,j,i-2) )                &
    20942324                                                )
    20952325
    2096           ibit14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp )
    2097           ibit13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp )
    2098           ibit12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp )
     2326          ibit5 = REAL( IBITS(advc_flags_m(k,j,i),5,1), KIND = wp )
     2327          ibit4 = REAL( IBITS(advc_flags_m(k,j,i),4,1), KIND = wp )
     2328          ibit3 = REAL( IBITS(advc_flags_m(k,j,i),3,1), KIND = wp )
    20992329
    21002330          v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2101           flux_n(k) = v_comp(k) * (                                           &
    2102                      ( 37.0_wp * ibit14 * adv_mom_5                           &
    2103                   +     7.0_wp * ibit13 * adv_mom_3                           &
    2104                   +              ibit12 * adv_mom_1                           &
    2105                      ) *                                                      &
    2106                                     ( u(k,j+1,i) + u(k,j,i)   )               &
    2107               -      (  8.0_wp * ibit14 * adv_mom_5                           &
    2108                   +              ibit13 * adv_mom_3                           &
    2109                      ) *                                                      &
    2110                                     ( u(k,j+2,i) + u(k,j-1,i) )               &
    2111               +      (           ibit14 * adv_mom_5                           &
    2112                      ) *                                                      &
    2113                                     ( u(k,j+3,i) + u(k,j-2,i) )               &
    2114                                   )
    2115 
    2116           diss_n(k) = - ABS ( v_comp(k) ) * (                                 &
    2117                      ( 10.0_wp * ibit14 * adv_mom_5                           &
    2118                   +     3.0_wp * ibit13 * adv_mom_3                           &
    2119                   +              ibit12 * adv_mom_1                           &
    2120                      ) *                                                      &
    2121                                     ( u(k,j+1,i) - u(k,j,i)   )               &
    2122               -      (  5.0_wp * ibit14 * adv_mom_5                           &
    2123                   +              ibit13 * adv_mom_3                           &
    2124                      ) *                                                      &
    2125                                     ( u(k,j+2,i) - u(k,j-1,i) )               &
    2126               +      (           ibit14 * adv_mom_5                           &
    2127                      ) *                                                      &
    2128                                     ( u(k,j+3,i) - u(k,j-2,i) )               &
     2331          flux_n(k) = v_comp(k) * (                                            &
     2332                     ( 37.0_wp * ibit5 * adv_mom_5                             &
     2333                  +     7.0_wp * ibit4 * adv_mom_3                             &
     2334                  +              ibit3 * adv_mom_1                             &
     2335                     ) *                                                       &
     2336                                    ( u(k,j+1,i) + u(k,j,i)   )                &
     2337              -      (  8.0_wp * ibit5 * adv_mom_5                             &
     2338                  +              ibit4 * adv_mom_3                             &
     2339                     ) *                                                       &
     2340                                    ( u(k,j+2,i) + u(k,j-1,i) )                &
     2341              +      (           ibit5 * adv_mom_5                             &
     2342                     ) *                                                       &
     2343                                    ( u(k,j+3,i) + u(k,j-2,i) )                &
     2344                                  )                                           
     2345                                                                             
     2346          diss_n(k) = - ABS ( v_comp(k) ) * (                                  &
     2347                     ( 10.0_wp * ibit5 * adv_mom_5                             &
     2348                  +     3.0_wp * ibit4 * adv_mom_3                             &
     2349                  +              ibit3 * adv_mom_1                             &
     2350                     ) *                                                       &
     2351                                    ( u(k,j+1,i) - u(k,j,i)   )                &
     2352              -      (  5.0_wp * ibit5 * adv_mom_5                             &
     2353                  +              ibit4 * adv_mom_3                             &
     2354                     ) *                                                       &
     2355                                    ( u(k,j+2,i) - u(k,j-1,i) )                &
     2356              +      (           ibit5 * adv_mom_5                             &
     2357                     ) *                                                       &
     2358                                    ( u(k,j+3,i) - u(k,j-2,i) )                &
    21292359                                            )
    21302360       ENDDO
     
    21332363
    21342364          u_comp(k) = u(k,j,i+1) + u(k,j,i)
    2135           flux_r(k) = ( u_comp(k) - gu ) * (                                  &
    2136                          37.0_wp * ( u(k,j,i+1) + u(k,j,i)   )                &
    2137                        -  8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) )                &
    2138                        +           ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
    2139           diss_r(k) = - ABS( u_comp(k) - gu ) * (                             &
    2140                          10.0_wp * ( u(k,j,i+1) - u(k,j,i)   )                &
    2141                        -  5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) )                &
    2142                        +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
    2143 
    2144           v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
    2145           flux_n(k) = v_comp(k) * (                                           &
    2146                          37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                &
    2147                        -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                &
    2148                        +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
    2149           diss_n(k) = - ABS( v_comp(k) ) * (                                  &
    2150                          10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                &
    2151                        -  5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) )                &
     2365          flux_r(k) = ( u_comp(k) - gu ) * (                                   &
     2366                         37.0_wp * ( u(k,j,i+1) + u(k,j,i)   )                 &
     2367                       -  8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) )                 &
     2368                       +           ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5   
     2369          diss_r(k) = - ABS( u_comp(k) - gu ) * (                              &
     2370                         10.0_wp * ( u(k,j,i+1) - u(k,j,i)   )                 &
     2371                       -  5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) )                 &
     2372                       +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5   
     2373                                                                               
     2374          v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv                           
     2375          flux_n(k) = v_comp(k) * (                                            &
     2376                         37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                 &
     2377                       -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                 &
     2378                       +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5   
     2379          diss_n(k) = - ABS( v_comp(k) ) * (                                   &
     2380                         10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                 &
     2381                       -  5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) )                 &
    21522382                       +           ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
    21532383
     
    21612391!--    calculated explicetely for the tendency at
    21622392!--    the first w-level. For topography wall this is done implicitely by
    2163 !--    advc_flags_1.
     2393!--    advc_flags_m.
    21642394       flux_t(nzb) = 0.0_wp
    21652395       diss_t(nzb) = 0.0_wp
     
    21692399!--       k index has to be modified near bottom and top, else array
    21702400!--       subscripts will be exceeded.
    2171           ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
    2172           ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
    2173           ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
    2174 
    2175           k_ppp = k + 3 * ibit17
    2176           k_pp  = k + 2 * ( 1 - ibit15 )
    2177           k_mm  = k - 2 * ibit17
     2401          ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
     2402          ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
     2403          ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp )
     2404
     2405          k_ppp = k + 3 * ibit8
     2406          k_pp  = k + 2 * ( 1 - ibit6 )
     2407          k_mm  = k - 2 * ibit8
    21782408
    21792409          w_comp(k) = w(k,j,i) + w(k,j,i-1)
    2180           flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    2181                      ( 37.0_wp * ibit17 * adv_mom_5                           &
    2182                   +     7.0_wp * ibit16 * adv_mom_3                           &
    2183                   +              ibit15 * adv_mom_1                           &
    2184                      ) *                                                      &
    2185                                 ( u(k+1,j,i)  + u(k,j,i)     )                &
    2186               -      (  8.0_wp * ibit17 * adv_mom_5                           &
    2187                   +              ibit16 * adv_mom_3                           &
    2188                      ) *                                                      &
    2189                                 ( u(k_pp,j,i) + u(k-1,j,i)   )                &
    2190               +      (           ibit17 * adv_mom_5                           &
    2191                      ) *                                                      &
    2192                                 ( u(k_ppp,j,i) + u(k_mm,j,i) )                &
    2193                                                   )
    2194 
    2195           diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    2196                      ( 10.0_wp * ibit17 * adv_mom_5                           &
    2197                   +     3.0_wp * ibit16 * adv_mom_3                           &
    2198                   +              ibit15 * adv_mom_1                           &
    2199                      ) *                                                      &
    2200                                 ( u(k+1,j,i)   - u(k,j,i)    )                &
    2201               -      (  5.0_wp * ibit17 * adv_mom_5                           &
    2202                   +              ibit16 * adv_mom_3                           &
    2203                      ) *                                                      &
    2204                                 ( u(k_pp,j,i)  - u(k-1,j,i)  )                &
    2205               +      (           ibit17 * adv_mom_5                           &
    2206                      ) *                                                      &
    2207                                 ( u(k_ppp,j,i) - u(k_mm,j,i) )                &
     2410          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                            &
     2411                     ( 37.0_wp * ibit8 * adv_mom_5                             &
     2412                  +     7.0_wp * ibit7 * adv_mom_3                             &
     2413                  +              ibit6 * adv_mom_1                             &
     2414                     ) *                                                       &
     2415                                ( u(k+1,j,i)   + u(k,j,i)    )                 &
     2416              -      (  8.0_wp * ibit8 * adv_mom_5                             &
     2417                  +              ibit7 * adv_mom_3                             &
     2418                     ) *                                                       &
     2419                                ( u(k_pp,j,i)  + u(k-1,j,i)  )                 &
     2420              +      (           ibit8 * adv_mom_5                             &
     2421                     ) *                                                       &
     2422                                ( u(k_ppp,j,i) + u(k_mm,j,i) )                 &
     2423                                                  )                           
     2424                                                                               
     2425          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                   &
     2426                     ( 10.0_wp * ibit8 * adv_mom_5                             &
     2427                  +     3.0_wp * ibit7 * adv_mom_3                             &
     2428                  +              ibit6 * adv_mom_1                             &
     2429                     ) *                                                       &
     2430                                ( u(k+1,j,i)   - u(k,j,i)    )                 &
     2431              -      (  5.0_wp * ibit8 * adv_mom_5                             &
     2432                  +              ibit7 * adv_mom_3                             &
     2433                     ) *                                                       &
     2434                                ( u(k_pp,j,i)  - u(k-1,j,i)  )                 &
     2435              +      (           ibit8 * adv_mom_5                             &
     2436                     ) *                                                       &
     2437                                ( u(k_ppp,j,i) - u(k_mm,j,i) )                 &
    22082438                                                           )
    22092439       ENDDO
     
    22112441       DO  k = nzb+3, nzt-2
    22122442
    2213           ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
    2214           ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
    2215           ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
     2443          ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
     2444          ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
     2445          ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp )
    22162446
    22172447          w_comp(k) = w(k,j,i) + w(k,j,i-1)
    2218           flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    2219                      ( 37.0_wp * ibit17 * adv_mom_5                           &
    2220                   +     7.0_wp * ibit16 * adv_mom_3                           &
    2221                   +              ibit15 * adv_mom_1                           &
    2222                      ) *                                                      &
    2223                                 ( u(k+1,j,i)  + u(k,j,i)     )                &
    2224               -      (  8.0_wp * ibit17 * adv_mom_5                           &
    2225                   +              ibit16 * adv_mom_3                           &
    2226                      ) *                                                      &
    2227                                 ( u(k+2,j,i) + u(k-1,j,i)   )                 &
    2228               +      (           ibit17 * adv_mom_5                           &
    2229                      ) *                                                      &
    2230                                 ( u(k+3,j,i) + u(k-2,j,i) )                   &
     2448          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                            &
     2449                     ( 37.0_wp * ibit8 * adv_mom_5                             &
     2450                  +     7.0_wp * ibit7 * adv_mom_3                             &
     2451                  +              ibit6 * adv_mom_1                             &
     2452                     ) *                                                       &
     2453                                ( u(k+1,j,i)  + u(k,j,i)     )                 &
     2454              -      (  8.0_wp * ibit8 * adv_mom_5                             &
     2455                  +              ibit7 * adv_mom_3                             &
     2456                     ) *                                                       &
     2457                                ( u(k+2,j,i) + u(k-1,j,i)   )                  &
     2458              +      (           ibit8 * adv_mom_5                             &
     2459                     ) *                                                       &
     2460                                ( u(k+3,j,i) + u(k-2,j,i) )                    &
    22312461                                                  )
    22322462
    2233           diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    2234                      ( 10.0_wp * ibit17 * adv_mom_5                           &
    2235                   +     3.0_wp * ibit16 * adv_mom_3                           &
    2236                   +              ibit15 * adv_mom_1                           &
    2237                      ) *                                                      &
    2238                                 ( u(k+1,j,i)   - u(k,j,i)    )                &
    2239               -      (  5.0_wp * ibit17 * adv_mom_5                           &
    2240                   +              ibit16 * adv_mom_3                           &
    2241                      ) *                                                      &
    2242                                 ( u(k+2,j,i)  - u(k-1,j,i)  )                 &
    2243               +      (           ibit17 * adv_mom_5                           &
    2244                      ) *                                                      &
    2245                                 ( u(k+3,j,i) - u(k-2,j,i) )                   &
     2463          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                   &
     2464                     ( 10.0_wp * ibit8 * adv_mom_5                             &
     2465                  +     3.0_wp * ibit7 * adv_mom_3                             &
     2466                  +              ibit6 * adv_mom_1                             &
     2467                     ) *                                                       &
     2468                                ( u(k+1,j,i)  - u(k,j,i)    )                  &
     2469              -      (  5.0_wp * ibit8 * adv_mom_5                             &
     2470                  +              ibit7 * adv_mom_3                             &
     2471                     ) *                                                       &
     2472                                ( u(k+2,j,i)  - u(k-1,j,i)  )                  &
     2473              +      (           ibit8 * adv_mom_5                             &
     2474                     ) *                                                       &
     2475                                ( u(k+3,j,i) - u(k-2,j,i)   )                  &
    22462476                                                           )
    22472477       ENDDO
     
    22512481!--       k index has to be modified near bottom and top, else array
    22522482!--       subscripts will be exceeded.
    2253           ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
    2254           ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
    2255           ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
    2256 
    2257           k_ppp = k + 3 * ibit17
    2258           k_pp  = k + 2 * ( 1 - ibit15 )
    2259           k_mm  = k - 2 * ibit17
     2483          ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
     2484          ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
     2485          ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp )
     2486
     2487          k_ppp = k + 3 * ibit8
     2488          k_pp  = k + 2 * ( 1 - ibit6 )
     2489          k_mm  = k - 2 * ibit8
    22602490
    22612491          w_comp(k) = w(k,j,i) + w(k,j,i-1)
    2262           flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    2263                      ( 37.0_wp * ibit17 * adv_mom_5                           &
    2264                   +     7.0_wp * ibit16 * adv_mom_3                           &
    2265                   +              ibit15 * adv_mom_1                           &
    2266                      ) *                                                      &
    2267                                 ( u(k+1,j,i)  + u(k,j,i)     )                &
    2268               -      (  8.0_wp * ibit17 * adv_mom_5                           &
    2269                   +              ibit16 * adv_mom_3                           &
    2270                      ) *                                                      &
    2271                                 ( u(k_pp,j,i) + u(k-1,j,i)   )                &
    2272               +      (           ibit17 * adv_mom_5                           &
    2273                      ) *                                                      &
    2274                                 ( u(k_ppp,j,i) + u(k_mm,j,i) )                &
     2492          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                            &
     2493                     ( 37.0_wp * ibit8 * adv_mom_5                             &
     2494                  +     7.0_wp * ibit7 * adv_mom_3                             &
     2495                  +              ibit6 * adv_mom_1                             &
     2496                     ) *                                                       &
     2497                                ( u(k+1,j,i)   + u(k,j,i)    )                 &
     2498              -      (  8.0_wp * ibit8 * adv_mom_5                             &
     2499                  +              ibit7 * adv_mom_3                             &
     2500                     ) *                                                       &
     2501                                ( u(k_pp,j,i)  + u(k-1,j,i)  )                 &
     2502              +      (           ibit8 * adv_mom_5                             &
     2503                     ) *                                                       &
     2504                                ( u(k_ppp,j,i) + u(k_mm,j,i) )                 &
    22752505                                                  )
    22762506
    2277           diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    2278                      ( 10.0_wp * ibit17 * adv_mom_5                           &
    2279                   +     3.0_wp * ibit16 * adv_mom_3                           &
    2280                   +              ibit15 * adv_mom_1                           &
    2281                      ) *                                                      &
    2282                                 ( u(k+1,j,i)   - u(k,j,i)    )                &
    2283               -      (  5.0_wp * ibit17 * adv_mom_5                           &
    2284                   +              ibit16 * adv_mom_3                           &
    2285                      ) *                                                      &
    2286                                 ( u(k_pp,j,i)  - u(k-1,j,i)  )                &
    2287               +      (           ibit17 * adv_mom_5                           &
    2288                      ) *                                                      &
    2289                                 ( u(k_ppp,j,i) - u(k_mm,j,i) )                &
     2507          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                   &
     2508                     ( 10.0_wp * ibit8 * adv_mom_5                             &
     2509                  +     3.0_wp * ibit7 * adv_mom_3                             &
     2510                  +              ibit6 * adv_mom_1                             &
     2511                     ) *                                                       &
     2512                                ( u(k+1,j,i)   - u(k,j,i)    )                 &
     2513              -      (  5.0_wp * ibit8 * adv_mom_5                             &
     2514                  +              ibit7 * adv_mom_3                             &
     2515                     ) *                                                       &
     2516                                ( u(k_pp,j,i)  - u(k-1,j,i)  )                 &
     2517              +      (           ibit8 * adv_mom_5                             &
     2518                     ) *                                                       &
     2519                                ( u(k_ppp,j,i) - u(k_mm,j,i) )                 &
    22902520                                                           )
    22912521       ENDDO
     
    22962526          diss_d    = diss_t(k-1)
    22972527         
    2298           ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp )
    2299           ibit10 = REAL( IBITS(advc_flags_1(k,j,i),10,1), KIND = wp )
    2300           ibit9  = REAL( IBITS(advc_flags_1(k,j,i),9,1), KIND = wp )
     2528          ibit2 = REAL( IBITS(advc_flags_m(k,j,i),2,1), KIND = wp )
     2529          ibit1 = REAL( IBITS(advc_flags_m(k,j,i),1,1), KIND = wp )
     2530          ibit0 = REAL( IBITS(advc_flags_m(k,j,i),0,1), KIND = wp )
    23012531         
    2302           ibit14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp )
    2303           ibit13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp )
    2304           ibit12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp )
     2532          ibit5 = REAL( IBITS(advc_flags_m(k,j,i),5,1), KIND = wp )
     2533          ibit4 = REAL( IBITS(advc_flags_m(k,j,i),4,1), KIND = wp )
     2534          ibit3 = REAL( IBITS(advc_flags_m(k,j,i),3,1), KIND = wp )
    23052535         
    2306           ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
    2307           ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
    2308           ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
     2536          ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp )
     2537          ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp )
     2538          ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp )
    23092539!
    23102540!--       Calculate the divergence of the velocity field. A respective
    23112541!--       correction is needed to overcome numerical instabilities introduced
    23122542!--       by a not sufficient reduction of divergences near topography.
    2313           div = ( ( u_comp(k)       * ( ibit9 + ibit10 + ibit11 )             &
    2314                 - ( u(k,j,i)   + u(k,j,i-1)   )                               &
    2315                                     * (                                       &
    2316                      REAL( IBITS(advc_flags_1(k,j,i-1),9,1),  KIND = wp )     &
    2317                    + REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )     &
    2318                    + REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )     &
    2319                                       )                                       &
    2320                   ) * ddx                                                     &
    2321                +  ( ( v_comp(k) + gv ) * ( ibit12 + ibit13 + ibit14 )         &
    2322                   - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
    2323                                     * (                                       &
    2324                      REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp )     &
    2325                    + REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp )     &
    2326                    + REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )     &
    2327                                       )                                       &
    2328                   ) * ddy                                                     &
    2329                +  ( w_comp(k)   * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )&
    2330                 -   w_comp(k-1) * rho_air_zw(k-1)                             &
    2331                                     * (                                       &
    2332                      REAL( IBITS(advc_flags_1(k-1,j,i),15,1), KIND = wp )     &
    2333                    + REAL( IBITS(advc_flags_1(k-1,j,i),16,1), KIND = wp )     &
    2334                    + REAL( IBITS(advc_flags_1(k-1,j,i),17,1), KIND = wp )     &
    2335                                       )                                       
    2336                   ) * drho_air(k) * ddzw(k)                                   &
    2337                 ) * 0.5_wp
    2338 
    2339           tend(k,j,i) = tend(k,j,i) - (                                       &
    2340                             ( flux_r(k) + diss_r(k)                           &
    2341                           -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx     &
    2342                           + ( flux_n(k) + diss_n(k)                           &
    2343                           -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
    2344                           + ( ( flux_t(k) + diss_t(k) )                       &
    2345                           -   ( flux_d    + diss_d    )                       &
    2346                                                     ) * drho_air(k) * ddzw(k) &
     2543          div = ( ( u_comp(k)       * ( ibit0 + ibit1 + ibit2 )                &
     2544                - ( u(k,j,i)   + u(k,j,i-1)   )                                &
     2545                                    * (                                        &
     2546                     REAL( IBITS(advc_flags_m(k,j,i-1),0,1),  KIND = wp )      &
     2547                   + REAL( IBITS(advc_flags_m(k,j,i-1),1,1), KIND = wp )       &
     2548                   + REAL( IBITS(advc_flags_m(k,j,i-1),2,1), KIND = wp )       &
     2549                                      )                                        &
     2550                  ) * ddx                                                      &
     2551               +  ( ( v_comp(k) + gv ) * ( ibit3 + ibit4 + ibit5 )             &
     2552                  - ( v(k,j,i)   + v(k,j,i-1 )  )                              &
     2553                                    * (                                        &
     2554                     REAL( IBITS(advc_flags_m(k,j-1,i),3,1), KIND = wp )       &
     2555                   + REAL( IBITS(advc_flags_m(k,j-1,i),4,1), KIND = wp )       &
     2556                   + REAL( IBITS(advc_flags_m(k,j-1,i),5,1), KIND = wp )       &
     2557                                      )                                        &
     2558                  ) * ddy                                                      &
     2559               +  ( w_comp(k)   * rho_air_zw(k) * ( ibit6 + ibit7 + ibit8 )    &
     2560                -   w_comp(k-1) * rho_air_zw(k-1)                              &
     2561                                    * (                                        &
     2562                     REAL( IBITS(advc_flags_m(k-1,j,i),6,1), KIND = wp )       &
     2563                   + REAL( IBITS(advc_flags_m(k-1,j,i),7,1), KIND = wp )       &
     2564                   + REAL( IBITS(advc_flags_m(k-1,j,i),8,1), KIND = wp )       &
     2565                                      )                                        
     2566                  ) * drho_air(k) * ddzw(k)                                    &
     2567                ) * 0.5_wp                                                     
     2568                                                                               
     2569          tend(k,j,i) = tend(k,j,i) - (                                        &
     2570                            ( flux_r(k) + diss_r(k)                            &
     2571                          -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx      &
     2572                          + ( flux_n(k) + diss_n(k)                            &
     2573                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy      &
     2574                          + ( ( flux_t(k) + diss_t(k) )                        &
     2575                          -   ( flux_d    + diss_d    )                        &
     2576                                                    ) * drho_air(k) * ddzw(k)  &
    23472577                                       ) + div * u(k,j,i)
    23482578
     
    24502680       INTEGER(iwp)  ::  tn        !< number of OpenMP thread
    24512681       
    2452        REAL(wp)      ::  ibit18   !< flag indicating 1st-order scheme along x-direction
    2453        REAL(wp)      ::  ibit19   !< flag indicating 3rd-order scheme along x-direction
    2454        REAL(wp)      ::  ibit20   !< flag indicating 5th-order scheme along x-direction
    2455        REAL(wp)      ::  ibit21   !< flag indicating 1st-order scheme along y-direction
    2456        REAL(wp)      ::  ibit22   !< flag indicating 3rd-order scheme along y-direction
    2457        REAL(wp)      ::  ibit23   !< flag indicating 3rd-order scheme along y-direction
    2458        REAL(wp)      ::  ibit24   !< flag indicating 1st-order scheme along z-direction
    2459        REAL(wp)      ::  ibit25   !< flag indicating 3rd-order scheme along z-direction
    2460        REAL(wp)      ::  ibit26   !< flag indicating 3rd-order scheme along z-direction
     2682       REAL(wp)      ::  ibit9    !< flag indicating 1st-order scheme along x-direction
     2683       REAL(wp)      ::  ibit10   !< flag indicating 3rd-order scheme along x-direction
     2684       REAL(wp)      ::  ibit11   !< flag indicating 5th-order scheme along x-direction
     2685       REAL(wp)      ::  ibit12   !< flag indicating 1st-order scheme along y-direction
     2686       REAL(wp)      ::  ibit13   !< flag indicating 3rd-order scheme along y-direction
     2687       REAL(wp)      ::  ibit14   !< flag indicating 3rd-order scheme along y-direction
     2688       REAL(wp)      ::  ibit15   !< flag indicating 1st-order scheme along z-direction
     2689       REAL(wp)      ::  ibit16   !< flag indicating 3rd-order scheme along z-direction
     2690       REAL(wp)      ::  ibit17   !< flag indicating 3rd-order scheme along z-direction
    24612691       REAL(wp)      ::  diss_d   !< artificial dissipation term at grid box bottom
    24622692       REAL(wp)      ::  div      !< divergence on v-grid
     
    24982728          DO  k = nzb+1, nzb_max_l
    24992729
    2500              ibit20 = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )
    2501              ibit19 = REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp )
    2502              ibit18 = REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp )
     2730             ibit11 = REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp )
     2731             ibit10 = REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp )
     2732             ibit9  = REAL( IBITS(advc_flags_m(k,j,i-1),9,1), KIND = wp )
    25032733
    25042734             u_comp(k)        = u(k,j-1,i) + u(k,j,i) - gu
    25052735             flux_l_v(k,j,tn) = u_comp(k) * (                                 &
    2506                               ( 37.0_wp * ibit20 * adv_mom_5                  &
    2507                            +     7.0_wp * ibit19 * adv_mom_3                  &
    2508                            +              ibit18 * adv_mom_1                  &
     2736                              ( 37.0_wp * ibit11 * adv_mom_5                  &
     2737                           +     7.0_wp * ibit10 * adv_mom_3                  &
     2738                           +              ibit9 * adv_mom_1                  &
    25092739                              ) *                                             &
    25102740                                        ( v(k,j,i)   + v(k,j,i-1) )           &
    2511                        -      (  8.0_wp * ibit20 * adv_mom_5                  &
    2512                            +              ibit19 * adv_mom_3                  &
     2741                       -      (  8.0_wp * ibit11 * adv_mom_5                  &
     2742                           +              ibit10 * adv_mom_3                  &
    25132743                              ) *                                             &
    25142744                                        ( v(k,j,i+1) + v(k,j,i-2) )           &
    2515                        +      (           ibit20 * adv_mom_5                  &
     2745                       +      (           ibit11 * adv_mom_5                  &
    25162746                              ) *                                             &
    25172747                                        ( v(k,j,i+2) + v(k,j,i-3) )           &
     
    25192749
    25202750              diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * (                       &
    2521                               ( 10.0_wp * ibit20 * adv_mom_5                  &
    2522                            +     3.0_wp * ibit19 * adv_mom_3                  &
    2523                            +              ibit18 * adv_mom_1                  &
     2751                              ( 10.0_wp * ibit11 * adv_mom_5                  &
     2752                           +     3.0_wp * ibit10 * adv_mom_3                  &
     2753                           +              ibit9 * adv_mom_1                  &
    25242754                              ) *                                             &
    25252755                                        ( v(k,j,i)   - v(k,j,i-1) )           &
    2526                        -      (  5.0_wp * ibit20 * adv_mom_5                  &
    2527                            +              ibit19 * adv_mom_3                  &
     2756                       -      (  5.0_wp * ibit11 * adv_mom_5                  &
     2757                           +              ibit10 * adv_mom_3                  &
    25282758                              ) *                                             &
    25292759                                        ( v(k,j,i+1) - v(k,j,i-2) )           &
    2530                        +      (           ibit20 * adv_mom_5                  &
     2760                       +      (           ibit11 * adv_mom_5                  &
    25312761                              ) *                                             &
    25322762                                        ( v(k,j,i+2) - v(k,j,i-3) )           &
     
    25372767          DO  k = nzb_max_l+1, nzt
    25382768
    2539              u_comp(k)           = u(k,j-1,i) + u(k,j,i) - gu
     2769             u_comp(k)        = u(k,j-1,i) + u(k,j,i) - gu
    25402770             flux_l_v(k,j,tn) = u_comp(k) * (                                 &
    2541                              37.0_wp * ( v(k,j,i) + v(k,j,i-1)   )            &
     2771                             37.0_wp * ( v(k,j,i)   + v(k,j,i-1)   )          &
    25422772                           -  8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) )            &
    25432773                           +           ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
    25442774             diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * (                        &
    2545                              10.0_wp * ( v(k,j,i) - v(k,j,i-1)   )            &
     2775                             10.0_wp * ( v(k,j,i)   - v(k,j,i-1)   )          &
    25462776                           -  5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) )            &
    25472777                           +           ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
     
    25562786          DO  k = nzb+1, nzb_max_l
    25572787
    2558              ibit23 = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )
    2559              ibit22 = REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp )
    2560              ibit21 = REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )
     2788             ibit14 = REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp )
     2789             ibit13 = REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp )
     2790             ibit12 = REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp )
    25612791
    25622792             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
    25632793             flux_s_v(k,tn) = v_comp_l * (                                    &
    2564                             ( 37.0_wp * ibit23 * adv_mom_5                    &
    2565                          +     7.0_wp * ibit22 * adv_mom_3                    &
    2566                          +              ibit21 * adv_mom_1                    &
     2794                            ( 37.0_wp * ibit14 * adv_mom_5                    &
     2795                         +     7.0_wp * ibit13 * adv_mom_3                    &
     2796                         +              ibit12 * adv_mom_1                    &
    25672797                            ) *                                               &
    25682798                                        ( v(k,j,i)   + v(k,j-1,i) )           &
    2569                      -      (  8.0_wp * ibit23 * adv_mom_5                    &
    2570                          +              ibit22 * adv_mom_3                    &
     2799                     -      (  8.0_wp * ibit14 * adv_mom_5                    &
     2800                         +              ibit13 * adv_mom_3                    &
    25712801                            ) *                                               &
    25722802                                        ( v(k,j+1,i) + v(k,j-2,i) )           &
    2573                      +      (           ibit23 * adv_mom_5                    &
     2803                     +      (           ibit14 * adv_mom_5                    &
    25742804                            ) *                                               &
    25752805                                        ( v(k,j+2,i) + v(k,j-3,i) )           &
     
    25772807
    25782808             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                           &
    2579                             ( 10.0_wp * ibit23 * adv_mom_5                    &
    2580                          +     3.0_wp * ibit22 * adv_mom_3                    &
    2581                          +              ibit21 * adv_mom_1                    &
     2809                            ( 10.0_wp * ibit14 * adv_mom_5                    &
     2810                         +     3.0_wp * ibit13 * adv_mom_3                    &
     2811                         +              ibit12 * adv_mom_1                    &
    25822812                            ) *                                               &
    25832813                                        ( v(k,j,i)   - v(k,j-1,i) )           &
    2584                      -      (  5.0_wp * ibit23 * adv_mom_5                    &
    2585                          +              ibit22 * adv_mom_3                    &
     2814                     -      (  5.0_wp * ibit14 * adv_mom_5                    &
     2815                         +              ibit13 * adv_mom_3                    &
    25862816                            ) *                                               &
    25872817                                        ( v(k,j+1,i) - v(k,j-2,i) )           &
    2588                      +      (           ibit23 * adv_mom_5                    &
     2818                     +      (           ibit14 * adv_mom_5                    &
    25892819                            ) *                                               &
    25902820                                        ( v(k,j+2,i) - v(k,j-3,i) )           &
     
    25972827             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
    25982828             flux_s_v(k,tn) = v_comp_l * (                                    &
    2599                            37.0_wp * ( v(k,j,i) + v(k,j-1,i)   )              &
     2829                           37.0_wp * ( v(k,j,i)   + v(k,j-1,i)   )            &
    26002830                         -  8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) )              &
    26012831                         +           ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
    26022832             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                           &
    2603                            10.0_wp * ( v(k,j,i) - v(k,j-1,i)   )              &
     2833                           10.0_wp * ( v(k,j,i)   - v(k,j-1,i)   )            &
    26042834                         -  5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) )              &
    26052835                         +           ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
     
    26132843       DO  k = nzb+1, nzb_max_l
    26142844
    2615           ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp )
    2616           ibit19 = REAL( IBITS(advc_flags_1(k,j,i),19,1), KIND = wp )
    2617           ibit18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp )
     2845          ibit11 = REAL( IBITS(advc_flags_m(k,j,i),11,1), KIND = wp )
     2846          ibit10 = REAL( IBITS(advc_flags_m(k,j,i),10,1), KIND = wp )
     2847          ibit9  = REAL( IBITS(advc_flags_m(k,j,i),9,1), KIND = wp )
    26182848 
    26192849          u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
    26202850          flux_r(k) = u_comp(k) * (                                           &
    2621                      ( 37.0_wp * ibit20 * adv_mom_5                           &
    2622                   +     7.0_wp * ibit19 * adv_mom_3                           &
    2623                   +              ibit18 * adv_mom_1                           &
     2851                     ( 37.0_wp * ibit11 * adv_mom_5                           &
     2852                  +     7.0_wp * ibit10 * adv_mom_3                           &
     2853                  +              ibit9 * adv_mom_1                           &
    26242854                     ) *                                                      &
    26252855                                    ( v(k,j,i+1) + v(k,j,i)   )               &
    2626               -      (  8.0_wp * ibit20 * adv_mom_5                           &
    2627                   +              ibit19 * adv_mom_3                           &
     2856              -      (  8.0_wp * ibit11 * adv_mom_5                           &
     2857                  +              ibit10 * adv_mom_3                           &
    26282858                     ) *                                                      &
    26292859                                    ( v(k,j,i+2) + v(k,j,i-1) )               &
    2630               +      (           ibit20 * adv_mom_5                           &
     2860              +      (           ibit11 * adv_mom_5                           &
    26312861                     ) *                                                      &
    26322862                                    ( v(k,j,i+3) + v(k,j,i-2) )               &
     
    26342864
    26352865          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    2636                      ( 10.0_wp * ibit20 * adv_mom_5                           &
    2637                   +     3.0_wp * ibit19 * adv_mom_3                           &
    2638                   +              ibit18 * adv_mom_1                           &
     2866                     ( 10.0_wp * ibit11 * adv_mom_5                           &
     2867                  +     3.0_wp * ibit10 * adv_mom_3                           &
     2868                  +              ibit9 * adv_mom_1                           &
    26392869                     ) *                                                      &
    26402870                                    ( v(k,j,i+1) - v(k,j,i)  )                &
    2641               -      (  5.0_wp * ibit20 * adv_mom_5                           &
    2642                   +              ibit19 * adv_mom_3                           &
     2871              -      (  5.0_wp * ibit11 * adv_mom_5                           &
     2872                  +              ibit10 * adv_mom_3                           &
    26432873                     ) *                                                      &
    26442874                                    ( v(k,j,i+2) - v(k,j,i-1) )               &
    2645               +      (           ibit20 * adv_mom_5                           &
     2875              +      (           ibit11 * adv_mom_5                           &
    26462876                     ) *                                                      &
    26472877                                    ( v(k,j,i+3) - v(k,j,i-2) )               &
    26482878                                           )
    26492879
    2650           ibit23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp )
    2651           ibit22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp )
    2652           ibit21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp )
     2880          ibit14 = REAL( IBITS(advc_flags_m(k,j,i),14,1), KIND = wp )
     2881          ibit13 = REAL( IBITS(advc_flags_m(k,j,i),13,1), KIND = wp )
     2882          ibit12 = REAL( IBITS(advc_flags_m(k,j,i),12,1), KIND = wp )
    26532883
    26542884
    26552885          v_comp(k) = v(k,j+1,i) + v(k,j,i)
    26562886          flux_n(k) = ( v_comp(k) - gv ) * (                                  &
    2657                      ( 37.0_wp * ibit23 * adv_mom_5                           &
    2658                   +     7.0_wp * ibit22 * adv_mom_3                           &
    2659                   +              ibit21 * adv_mom_1                           &
     2887                     ( 37.0_wp * ibit14 * adv_mom_5                           &
     2888                  +     7.0_wp * ibit13 * adv_mom_3                           &
     2889                  +              ibit12 * adv_mom_1                           &
    26602890                     ) *                                                      &
    26612891                                    ( v(k,j+1,i) + v(k,j,i)   )               &
    2662               -      (  8.0_wp * ibit23 * adv_mom_5                           &
    2663                   +              ibit22 * adv_mom_3                           &
     2892              -      (  8.0_wp * ibit14 * adv_mom_5                           &
     2893                  +              ibit13 * adv_mom_3                           &
    26642894                     ) *                                                      &
    26652895                                    ( v(k,j+2,i) + v(k,j-1,i) )               &
    2666               +      (           ibit23 * adv_mom_5                           &
     2896              +      (           ibit14 * adv_mom_5                           &
    26672897                     ) *                                                      &
    26682898                                    ( v(k,j+3,i) + v(k,j-2,i) )               &
     
    26702900
    26712901          diss_n(k) = - ABS( v_comp(k) - gv ) * (                             &
    2672                      ( 10.0_wp * ibit23 * adv_mom_5                           &
    2673                   +     3.0_wp * ibit22 * adv_mom_3                           &
    2674                   +              ibit21 * adv_mom_1                           &
     2902                     ( 10.0_wp * ibit14 * adv_mom_5                           &
     2903                  +     3.0_wp * ibit13 * adv_mom_3                           &
     2904                  +              ibit12 * adv_mom_1                           &
    26752905                     ) *                                                      &
    26762906                                    ( v(k,j+1,i) - v(k,j,i)   )               &
    2677               -      (  5.0_wp * ibit23 * adv_mom_5                           &
    2678                   +              ibit22 * adv_mom_3                           &
     2907              -      (  5.0_wp * ibit14 * adv_mom_5                           &
     2908                  +              ibit13 * adv_mom_3                           &
    26792909                     ) *                                                      &
    26802910                                    ( v(k,j+2,i) - v(k,j-1,i) )               &
    2681               +      (           ibit23 * adv_mom_5                           &
     2911              +      (           ibit14 * adv_mom_5                           &
    26822912                     ) *                                                      &
    26832913                                    ( v(k,j+3,i) - v(k,j-2,i) )               &
     
    27182948!--    calculated explicetely for the tendency at
    27192949!--    the first w-level. For topography wall this is done implicitely by
    2720 !--    advc_flags_1.
     2950!--    advc_flags_m.
    27212951       flux_t(nzb) = 0.0_wp
    27222952       diss_t(nzb) = 0.0_wp
     
    27272957!--       k index has to be modified near bottom and top, else array
    27282958!--       subscripts will be exceeded.
    2729           ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
    2730           ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
    2731           ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
    2732 
    2733           k_ppp = k + 3 * ibit26
    2734           k_pp  = k + 2 * ( 1 - ibit24  )
    2735           k_mm  = k - 2 * ibit26
     2959          ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp )
     2960          ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp )
     2961          ibit15 = REAL( IBITS(advc_flags_m(k,j,i),15,1), KIND = wp )
     2962
     2963          k_ppp = k + 3 * ibit17
     2964          k_pp  = k + 2 * ( 1 - ibit15  )
     2965          k_mm  = k - 2 * ibit17
    27362966
    27372967          w_comp(k) = w(k,j-1,i) + w(k,j,i)
    27382968          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    2739                      ( 37.0_wp * ibit26 * adv_mom_5                           &
    2740                   +     7.0_wp * ibit25 * adv_mom_3                           &
    2741                   +              ibit24 * adv_mom_1                           &
     2969                     ( 37.0_wp * ibit17 * adv_mom_5                           &
     2970                  +     7.0_wp * ibit16 * adv_mom_3                           &
     2971                  +              ibit15 * adv_mom_1                           &
    27422972                     ) *                                                      &
    27432973                                ( v(k+1,j,i)   + v(k,j,i)    )                &
    2744               -      (  8.0_wp * ibit26 * adv_mom_5                           &
    2745                   +              ibit25 * adv_mom_3                           &
     2974              -      (  8.0_wp * ibit17 * adv_mom_5                           &
     2975                  +              ibit16 * adv_mom_3                           &
    27462976                     ) *                                                      &
    27472977                                ( v(k_pp,j,i)  + v(k-1,j,i)  )                &
    2748               +      (           ibit26 * adv_mom_5                           &
     2978              +      (           ibit17 * adv_mom_5                           &
    27492979                     ) *                                                      &
    27502980                                ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
     
    27522982
    27532983          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    2754                      ( 10.0_wp * ibit26 * adv_mom_5                           &
    2755                   +     3.0_wp * ibit25 * adv_mom_3                           &
    2756                   +              ibit24 * adv_mom_1                           &
     2984                     ( 10.0_wp * ibit17 * adv_mom_5                           &
     2985                  +     3.0_wp * ibit16 * adv_mom_3                           &
     2986                  +              ibit15 * adv_mom_1                           &
    27572987                     ) *                                                      &
    27582988                                ( v(k+1,j,i)   - v(k,j,i)    )                &
    2759               -      (  5.0_wp * ibit26 * adv_mom_5                           &
    2760                   +              ibit25 * adv_mom_3                           &
     2989              -      (  5.0_wp * ibit17 * adv_mom_5                           &
     2990                  +              ibit16 * adv_mom_3                           &
    27612991                     ) *                                                      &
    27622992                                ( v(k_pp,j,i)  - v(k-1,j,i)  )                &
    2763               +      (           ibit26 * adv_mom_5                           &
     2993              +      (           ibit17 * adv_mom_5                           &
    27642994                     ) *                                                      &
    27652995                                ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
     
    27692999       DO  k = nzb+3, nzt-2
    27703000
    2771           ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
    2772           ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
    2773           ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
     3001          ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp )
     3002          ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp )
     3003          ibit15 = REAL( IBITS(advc_flags_m(k,j,i),15,1), KIND = wp )
    27743004
    27753005          w_comp(k) = w(k,j-1,i) + w(k,j,i)
    27763006          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    2777                      ( 37.0_wp * ibit26 * adv_mom_5                           &
    2778                   +     7.0_wp * ibit25 * adv_mom_3                           &
    2779                   +              ibit24 * adv_mom_1                           &
    2780                      ) *                                                      &
    2781                                 ( v(k+1,j,i)   + v(k,j,i)    )                &
    2782               -      (  8.0_wp * ibit26 * adv_mom_5                           &
    2783                   +              ibit25 * adv_mom_3                           &
    2784                      ) *                                                      &
    2785                                 ( v(k+2,j,i)  + v(k-1,j,i)  )                 &
    2786               +      (           ibit26 * adv_mom_5                           &
     3007                     ( 37.0_wp * ibit17 * adv_mom_5                           &
     3008                  +     7.0_wp * ibit16 * adv_mom_3                           &
     3009                  +              ibit15 * adv_mom_1                           &
     3010                     ) *                                                      &
     3011                                ( v(k+1,j,i) + v(k,j,i)   )                   &
     3012              -      (  8.0_wp * ibit17 * adv_mom_5                           &
     3013                  +              ibit16 * adv_mom_3                           &
     3014                     ) *                                                      &
     3015                                ( v(k+2,j,i) + v(k-1,j,i) )                   &
     3016              +      (           ibit17 * adv_mom_5                           &
    27873017                     ) *                                                      &
    27883018                                ( v(k+3,j,i) + v(k-2,j,i) )                   &
     
    27903020
    27913021          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    2792                      ( 10.0_wp * ibit26 * adv_mom_5                           &
    2793                   +     3.0_wp * ibit25 * adv_mom_3                           &
    2794                   +              ibit24 * adv_mom_1                           &
    2795                      ) *                                                      &
    2796                                 ( v(k+1,j,i)   - v(k,j,i)    )                &
    2797               -      (  5.0_wp * ibit26 * adv_mom_5                           &
    2798                   +              ibit25 * adv_mom_3                           &
    2799                      ) *                                                      &
    2800                                 ( v(k+2,j,i)  - v(k-1,j,i)  )                 &
    2801               +      (           ibit26 * adv_mom_5                           &
     3022                     ( 10.0_wp * ibit17 * adv_mom_5                           &
     3023                  +     3.0_wp * ibit16 * adv_mom_3                           &
     3024                  +              ibit15 * adv_mom_1                           &
     3025                     ) *                                                      &
     3026                                ( v(k+1,j,i) - v(k,j,i)   )                   &
     3027              -      (  5.0_wp * ibit17 * adv_mom_5                           &
     3028                  +              ibit16 * adv_mom_3                           &
     3029                     ) *                                                      &
     3030                                ( v(k+2,j,i) - v(k-1,j,i) )                    &
     3031              +      (           ibit17 * adv_mom_5                           &
    28023032                     ) *                                                      &
    28033033                                ( v(k+3,j,i) - v(k-2,j,i) )                   &
     
    28093039!--       k index has to be modified near bottom and top, else array
    28103040!--       subscripts will be exceeded.
    2811           ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
    2812           ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
    2813           ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
    2814 
    2815           k_ppp = k + 3 * ibit26
    2816           k_pp  = k + 2 * ( 1 - ibit24  )
    2817           k_mm  = k - 2 * ibit26
     3041          ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp )
     3042          ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp )
     3043          ibit15 = REAL( IBITS(advc_flags_m(k,j,i),15,1), KIND = wp )
     3044
     3045          k_ppp = k + 3 * ibit17
     3046          k_pp  = k + 2 * ( 1 - ibit15  )
     3047          k_mm  = k - 2 * ibit17
    28183048
    28193049          w_comp(k) = w(k,j-1,i) + w(k,j,i)
    28203050          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
    2821                      ( 37.0_wp * ibit26 * adv_mom_5                           &
    2822                   +     7.0_wp * ibit25 * adv_mom_3                           &
    2823                   +              ibit24 * adv_mom_1                           &
     3051                     ( 37.0_wp * ibit17 * adv_mom_5                           &
     3052                  +     7.0_wp * ibit16 * adv_mom_3                           &
     3053                  +              ibit15 * adv_mom_1                           &
    28243054                     ) *                                                      &
    28253055                                ( v(k+1,j,i)   + v(k,j,i)    )                &
    2826               -      (  8.0_wp * ibit26 * adv_mom_5                           &
    2827                   +              ibit25 * adv_mom_3                           &
     3056              -      (  8.0_wp * ibit17 * adv_mom_5                           &
     3057                  +              ibit16 * adv_mom_3                           &
    28283058                     ) *                                                      &
    28293059                                ( v(k_pp,j,i)  + v(k-1,j,i)  )                &
    2830               +      (           ibit26 * adv_mom_5                           &
     3060              +      (           ibit17 * adv_mom_5                           &
    28313061                     ) *                                                      &
    28323062                                ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
     
    28343064
    28353065          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
    2836                      ( 10.0_wp * ibit26 * adv_mom_5                           &
    2837                   +     3.0_wp * ibit25 * adv_mom_3                           &
    2838                   +              ibit24 * adv_mom_1                           &
     3066                     ( 10.0_wp * ibit17 * adv_mom_5                           &
     3067                  +     3.0_wp * ibit16 * adv_mom_3                           &
     3068                  +              ibit15 * adv_mom_1                           &
    28393069                     ) *                                                      &
    28403070                                ( v(k+1,j,i)   - v(k,j,i)    )                &
    2841               -      (  5.0_wp * ibit26 * adv_mom_5                           &
    2842                   +              ibit25 * adv_mom_3                           &
     3071              -      (  5.0_wp * ibit17 * adv_mom_5                           &
     3072                  +              ibit16 * adv_mom_3                           &
    28433073                     ) *                                                      &
    28443074                                ( v(k_pp,j,i)  - v(k-1,j,i)  )                &
    2845               +      (           ibit26 * adv_mom_5                           &
     3075              +      (           ibit17 * adv_mom_5                           &
    28463076                     ) *                                                      &
    28473077                                ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
     
    28543084          diss_d    = diss_t(k-1)
    28553085         
    2856           ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp )
    2857           ibit19 = REAL( IBITS(advc_flags_1(k,j,i),19,1), KIND = wp )
    2858           ibit18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp )
     3086          ibit11 = REAL( IBITS(advc_flags_m(k,j,i),11,1), KIND = wp )
     3087          ibit10 = REAL( IBITS(advc_flags_m(k,j,i),10,1), KIND = wp )
     3088          ibit9  = REAL( IBITS(advc_flags_m(k,j,i),9,1), KIND = wp )
    28593089         
    2860           ibit23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp )
    2861           ibit22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp )
    2862           ibit21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp )
     3090          ibit14 = REAL( IBITS(advc_flags_m(k,j,i),14,1), KIND = wp )
     3091          ibit13 = REAL( IBITS(advc_flags_m(k,j,i),13,1), KIND = wp )
     3092          ibit12 = REAL( IBITS(advc_flags_m(k,j,i),12,1), KIND = wp )
    28633093         
    2864           ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
    2865           ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
    2866           ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
     3094          ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp )
     3095          ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp )
     3096          ibit15 = REAL( IBITS(advc_flags_m(k,j,i),15,1), KIND = wp )
    28673097!
    28683098!--       Calculate the divergence of the velocity field. A respective
     
    28703100!--       by a not sufficient reduction of divergences near topography.
    28713101          div = ( ( ( u_comp(k)     + gu )                                    &
    2872                                        * ( ibit18 + ibit19 + ibit20 )         &
     3102                                       * ( ibit9 + ibit10 + ibit11 )          &
    28733103                  - ( u(k,j-1,i) + u(k,j,i) )                                 &
    28743104                                       * (                                    &
    2875                         REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp )  &
    2876                       + REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp )  &
    2877                       + REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )  &
     3105                        REAL( IBITS(advc_flags_m(k,j,i-1),9,1), KIND = wp )  &
     3106                      + REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp )  &
     3107                      + REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp )  &
    28783108                                         )                                    &
    28793109                  ) * ddx                                                     &
    28803110               +  ( v_comp(k)                                                 &
    2881                                        * ( ibit21 + ibit22 + ibit23 )         &
     3111                                       * ( ibit12 + ibit13 + ibit14 )         &
    28823112                - ( v(k,j,i)     + v(k,j-1,i) )                               &
    28833113                                       * (                                    &
    2884                         REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )  &
    2885                       + REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp )  &
    2886                       + REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )  &
     3114                        REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp )  &
     3115                      + REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp )  &
     3116                      + REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp )  &
    28873117                                         )                                    &
    28883118                  ) * ddy                                                     &
    2889                +  ( w_comp(k)   * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )&
     3119               +  ( w_comp(k)   * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )&
    28903120                -   w_comp(k-1) * rho_air_zw(k-1)                             &
    28913121                                       * (                                    &
    2892                         REAL( IBITS(advc_flags_1(k-1,j,i),24,1), KIND = wp )  &
    2893                       + REAL( IBITS(advc_flags_1(k-1,j,i),25,1), KIND = wp )  &
    2894                       + REAL( IBITS(advc_flags_1(k-1,j,i),26,1), KIND = wp )  &
     3122                        REAL( IBITS(advc_flags_m(k-1,j,i),15,1), KIND = wp )  &
     3123                      + REAL( IBITS(advc_flags_m(k-1,j,i),16,1), KIND = wp )  &
     3124                      + REAL( IBITS(advc_flags_m(k-1,j,i),17,1), KIND = wp )  &
    28953125                                         )                                    &
    28963126                   ) * drho_air(k) * ddzw(k)                                  &
     
    30113241       INTEGER(iwp) ::  tn        !< number of OpenMP thread
    30123242       
    3013        REAL(wp)    ::  ibit27  !< flag indicating 1st-order scheme along x-direction
    3014        REAL(wp)    ::  ibit28  !< flag indicating 3rd-order scheme along x-direction
    3015        REAL(wp)    ::  ibit29  !< flag indicating 5th-order scheme along x-direction
    3016        REAL(wp)    ::  ibit30  !< flag indicating 1st-order scheme along y-direction
    3017        REAL(wp)    ::  ibit31  !< flag indicating 3rd-order scheme along y-direction
    3018        REAL(wp)    ::  ibit32  !< flag indicating 5th-order scheme along y-direction
    3019        REAL(wp)    ::  ibit33  !< flag indicating 1st-order scheme along z-direction
    3020        REAL(wp)    ::  ibit34  !< flag indicating 3rd-order scheme along z-direction
    3021        REAL(wp)    ::  ibit35  !< flag indicating 5th-order scheme along z-direction
     3243       REAL(wp)    ::  ibit18  !< flag indicating 1st-order scheme along x-direction
     3244       REAL(wp)    ::  ibit19  !< flag indicating 3rd-order scheme along x-direction
     3245       REAL(wp)    ::  ibit20  !< flag indicating 5th-order scheme along x-direction
     3246       REAL(wp)    ::  ibit21  !< flag indicating 1st-order scheme along y-direction
     3247       REAL(wp)    ::  ibit22  !< flag indicating 3rd-order scheme along y-direction
     3248       REAL(wp)    ::  ibit23  !< flag indicating 5th-order scheme along y-direction
     3249       REAL(wp)    ::  ibit24  !< flag indicating 1st-order scheme along z-direction
     3250       REAL(wp)    ::  ibit25  !< flag indicating 3rd-order scheme along z-direction
     3251       REAL(wp)    ::  ibit26  !< flag indicating 5th-order scheme along z-direction
    30223252       REAL(wp)    ::  diss_d  !< discretized artificial dissipation at top of the grid box
    30233253       REAL(wp)    ::  div     !< divergence on w-grid
     
    30563286
    30573287          DO  k = nzb+1, nzb_max_l
    3058              ibit32 = REAL( IBITS(advc_flags_2(k,j-1,i),0,1), KIND = wp )
    3059              ibit31 = REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )
    3060              ibit30 = REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp )
     3288             ibit23 = REAL( IBITS(advc_flags_m(k,j-1,i),23,1), KIND = wp )
     3289             ibit22 = REAL( IBITS(advc_flags_m(k,j-1,i),22,1), KIND = wp )
     3290             ibit21 = REAL( IBITS(advc_flags_m(k,j-1,i),21,1), KIND = wp )
    30613291
    30623292             v_comp(k)      = v(k+1,j,i) + v(k,j,i) - gv
    30633293             flux_s_w(k,tn) = v_comp(k) * (                                   &
    3064                             ( 37.0_wp * ibit32 * adv_mom_5                    &
    3065                          +     7.0_wp * ibit31 * adv_mom_3                    &
    3066                          +              ibit30 * adv_mom_1                    &
     3294                            ( 37.0_wp * ibit23 * adv_mom_5                    &
     3295                         +     7.0_wp * ibit22 * adv_mom_3                    &
     3296                         +              ibit21 * adv_mom_1                    &
    30673297                            ) *                                               &
    30683298                                        ( w(k,j,i)   + w(k,j-1,i) )           &
    3069                      -      (  8.0_wp * ibit32 * adv_mom_5                    &
    3070                          +              ibit31 * adv_mom_3                    &
     3299                     -      (  8.0_wp * ibit23 * adv_mom_5                    &
     3300                         +              ibit22 * adv_mom_3                    &
    30713301                            ) *                                               &
    30723302                                        ( w(k,j+1,i) + w(k,j-2,i) )           &
    3073                      +      (           ibit32 * adv_mom_5                    &
     3303                     +      (           ibit23 * adv_mom_5                    &
    30743304                            ) *                                               &
    30753305                                        ( w(k,j+2,i) + w(k,j-3,i) )           &
     
    30773307
    30783308             diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                          &
    3079                             ( 10.0_wp * ibit32 * adv_mom_5                    &
    3080                          +     3.0_wp * ibit31 * adv_mom_3                    &
    3081                          +              ibit30 * adv_mom_1                    &
     3309                            ( 10.0_wp * ibit23 * adv_mom_5                    &
     3310                         +     3.0_wp * ibit22 * adv_mom_3                    &
     3311                         +              ibit21 * adv_mom_1                    &
    30823312                            ) *                                               &
    30833313                                        ( w(k,j,i)   - w(k,j-1,i) )           &
    3084                      -      (  5.0_wp * ibit32 * adv_mom_5                    &
    3085                          +              ibit31 * adv_mom_3                    &
     3314                     -      (  5.0_wp * ibit23 * adv_mom_5                    &
     3315                         +              ibit22 * adv_mom_3                    &
    30863316                            ) *                                               &
    30873317                                        ( w(k,j+1,i) - w(k,j-2,i) )           &
    3088                      +      (           ibit32 * adv_mom_5                    &
     3318                     +      (           ibit23 * adv_mom_5                    &
    30893319                            ) *                                               &
    30903320                                        ( w(k,j+2,i) - w(k,j-3,i) )           &
     
    30973327             v_comp(k)      = v(k+1,j,i) + v(k,j,i) - gv
    30983328             flux_s_w(k,tn) = v_comp(k) * (                                   &
    3099                            37.0_wp * ( w(k,j,i) + w(k,j-1,i)  )              &
    3100                          -  8.0_wp * ( w(k,j+1,i) +w(k,j-2,i) )              &
     3329                           37.0_wp * ( w(k,j,i)   + w(k,j-1,i) )              &
     3330                         -  8.0_wp * ( w(k,j+1,i) + w(k,j-2,i) )              &
    31013331                         +           ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
    31023332             diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                          &
    3103                            10.0_wp * ( w(k,j,i) - w(k,j-1,i)  )              &
     3333                           10.0_wp * ( w(k,j,i)   - w(k,j-1,i) )              &
    31043334                         -  5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) )              &
    31053335                         +           ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
     
    31143344          DO  k = nzb+1, nzb_max_l
    31153345
    3116              ibit29 = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )
    3117              ibit28 = REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp )
    3118              ibit27 = REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp )
     3346             ibit20 = REAL( IBITS(advc_flags_m(k,j,i-1),20,1), KIND = wp )
     3347             ibit19 = REAL( IBITS(advc_flags_m(k,j,i-1),19,1), KIND = wp )
     3348             ibit18 = REAL( IBITS(advc_flags_m(k,j,i-1),18,1), KIND = wp )
    31193349
    31203350             u_comp(k)        = u(k+1,j,i) + u(k,j,i) - gu
    31213351             flux_l_w(k,j,tn) = u_comp(k) * (                                 &
    3122                              ( 37.0_wp * ibit29 * adv_mom_5                   &
    3123                           +     7.0_wp * ibit28 * adv_mom_3                   &
    3124                           +              ibit27 * adv_mom_1                   &
     3352                             ( 37.0_wp * ibit20 * adv_mom_5                   &
     3353                          +     7.0_wp * ibit19 * adv_mom_3                   &
     3354                          +              ibit18 * adv_mom_1                   &
    31253355                             ) *                                              &
    31263356                                        ( w(k,j,i)   + w(k,j,i-1) )           &
    3127                       -      (  8.0_wp * ibit29 * adv_mom_5                   &
    3128                           +              ibit28 * adv_mom_3                   &
     3357                      -      (  8.0_wp * ibit20 * adv_mom_5                   &
     3358                          +              ibit19 * adv_mom_3                   &
    31293359                             ) *                                              &
    31303360                                        ( w(k,j,i+1) + w(k,j,i-2) )           &
    3131                       +      (           ibit29 * adv_mom_5                   &
     3361                      +      (           ibit20 * adv_mom_5                   &
    31323362                             ) *                                              &
    31333363                                        ( w(k,j,i+2) + w(k,j,i-3) )           &
     
    31353365
    31363366               diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * (                      &
    3137                              ( 10.0_wp * ibit29 * adv_mom_5                   &
    3138                           +     3.0_wp * ibit28 * adv_mom_3                   &
    3139                           +              ibit27 * adv_mom_1                   &
     3367                             ( 10.0_wp * ibit20 * adv_mom_5                   &
     3368                          +     3.0_wp * ibit19 * adv_mom_3                   &
     3369                          +              ibit18 * adv_mom_1                   &
    31403370                             ) *                                              &
    31413371                                        ( w(k,j,i)   - w(k,j,i-1) )           &
    3142                       -      (  5.0_wp * ibit29 * adv_mom_5                   &
    3143                           +              ibit28 * adv_mom_3                   &
     3372                      -      (  5.0_wp * ibit20 * adv_mom_5                   &
     3373                          +              ibit19 * adv_mom_3                   &
    31443374                             ) *                                              &
    31453375                                        ( w(k,j,i+1) - w(k,j,i-2) )           &
    3146                       +      (           ibit29 * adv_mom_5                   &
     3376                      +      (           ibit20 * adv_mom_5                   &
    31473377                             ) *                                              &
    31483378                                        ( w(k,j,i+2) - w(k,j,i-3) )           &
     
    31553385             u_comp(k)        = u(k+1,j,i) + u(k,j,i) - gu
    31563386             flux_l_w(k,j,tn) = u_comp(k) * (                                 &
    3157                             37.0_wp * ( w(k,j,i) + w(k,j,i-1)  )             &
     3387                            37.0_wp * ( w(k,j,i)   + w(k,j,i-1) )             &
    31583388                          -  8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) )             &
    31593389                          +           ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
    31603390             diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * (                        &
    3161                             10.0_wp * ( w(k,j,i) - w(k,j,i-1)  )             &
     3391                            10.0_wp * ( w(k,j,i)   - w(k,j,i-1) )             &
    31623392                          -  5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) )             &
    31633393                          +           ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
     
    31713401       DO  k = nzb+1, nzb_max_l
    31723402
    3173           ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp )
    3174           ibit28 = REAL( IBITS(advc_flags_1(k,j,i),28,1), KIND = wp )
    3175           ibit27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp )
     3403          ibit20 = REAL( IBITS(advc_flags_m(k,j,i),20,1), KIND = wp )
     3404          ibit19 = REAL( IBITS(advc_flags_m(k,j,i),19,1), KIND = wp )
     3405          ibit18 = REAL( IBITS(advc_flags_m(k,j,i),18,1), KIND = wp )
    31763406
    31773407          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
    31783408          flux_r(k) = u_comp(k) * (                                           &
    3179                      ( 37.0_wp * ibit29 * adv_mom_5                           &
    3180                   +     7.0_wp * ibit28 * adv_mom_3                           &
    3181                   +              ibit27 * adv_mom_1                           &
     3409                     ( 37.0_wp * ibit20 * adv_mom_5                           &
     3410                  +     7.0_wp * ibit19 * adv_mom_3                           &
     3411                  +              ibit18 * adv_mom_1                           &
    31823412                     ) *                                                      &
    31833413                                    ( w(k,j,i+1) + w(k,j,i)   )               &
    3184               -      (  8.0_wp * ibit29 * adv_mom_5                           &
    3185                   +              ibit28 * adv_mom_3                           &
     3414              -      (  8.0_wp * ibit20 * adv_mom_5                           &
     3415                  +              ibit19 * adv_mom_3                           &
    31863416                     ) *                                                      &
    31873417                                    ( w(k,j,i+2) + w(k,j,i-1) )               &
    3188               +      (           ibit29 * adv_mom_5                           &
     3418              +      (           ibit20 * adv_mom_5                           &
    31893419                     ) *                                                      &
    31903420                                    ( w(k,j,i+3) + w(k,j,i-2) )               &
     
    31923422
    31933423          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
    3194                      ( 10.0_wp * ibit29 * adv_mom_5                           &
    3195                   +     3.0_wp * ibit28 * adv_mom_3                           &
    3196                   +              ibit27 * adv_mom_1                           &
     3424                     ( 10.0_wp * ibit20 * adv_mom_5                           &
     3425                  +     3.0_wp * ibit19 * adv_mom_3                           &
     3426                  +              ibit18 * adv_mom_1                           &
    31973427                     ) *                                                      &
    31983428                                    ( w(k,j,i+1) - w(k,j,i)   )               &
    3199               -      (  5.0_wp * ibit29 * adv_mom_5                           &
    3200                   +              ibit28 * adv_mom_3                           &
     3429              -      (  5.0_wp * ibit20 * adv_mom_5                           &
     3430                  +              ibit19 * adv_mom_3                           &
    32013431                     ) *                                                      &
    32023432                                    ( w(k,j,i+2) - w(k,j,i-1) )               &
    3203               +      (           ibit29 * adv_mom_5                           &
     3433              +      (           ibit20 * adv_mom_5                           &
    32043434                     ) *                                                      &
    32053435                                    ( w(k,j,i+3) - w(k,j,i-2) )               &
    32063436                                           )
    32073437
    3208           ibit32 = REAL( IBITS(advc_flags_2(k,j,i),0,1), KIND = wp )
    3209           ibit31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp )
    3210           ibit30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp )
     3438          ibit23 = REAL( IBITS(advc_flags_m(k,j,i),23,1), KIND = wp )
     3439          ibit22 = REAL( IBITS(advc_flags_m(k,j,i),22,1), KIND = wp )
     3440          ibit21 = REAL( IBITS(advc_flags_m(k,j,i),21,1), KIND = wp )
    32113441
    32123442          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
    32133443          flux_n(k) = v_comp(k) * (                                           &
    3214                      ( 37.0_wp * ibit32 * adv_mom_5                           &
    3215                   +     7.0_wp * ibit31 * adv_mom_3                           &
    3216                   +              ibit30 * adv_mom_1                           &
     3444                     ( 37.0_wp * ibit23 * adv_mom_5                           &
     3445                  +     7.0_wp * ibit22 * adv_mom_3                           &
     3446                  +              ibit21 * adv_mom_1                           &
    32173447                     ) *                                                      &
    32183448                                    ( w(k,j+1,i) + w(k,j,i)   )               &
    3219               -      (  8.0_wp * ibit32 * adv_mom_5                           &
    3220                   +              ibit31 * adv_mom_3                           &
     3449              -      (  8.0_wp * ibit23 * adv_mom_5                           &
     3450                  +              ibit22 * adv_mom_3                           &
    32213451                     ) *                                                      &
    32223452                                    ( w(k,j+2,i) + w(k,j-1,i) )               &
    3223               +      (           ibit32 * adv_mom_5                           &
     3453              +      (           ibit23 * adv_mom_5                           &
    32243454                     ) *                                                      &
    32253455                                    ( w(k,j+3,i) + w(k,j-2,i) )               &
     
    32273457
    32283458          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
    3229                      ( 10.0_wp * ibit32 * adv_mom_5                           &
    3230                   +     3.0_wp * ibit31 * adv_mom_3                           &
    3231                   +              ibit30 * adv_mom_1                           &
    3232                      ) *                                                      &
    3233                                     ( w(k,j+1,i) - w(k,j,i)  )                &
    3234               -      (  5.0_wp * ibit32 * adv_mom_5                           &
    3235                   +              ibit31 * adv_mom_3                           &
    3236                      ) *                                                      &
    3237                                    ( w(k,j+2,i) - w(k,j-1,i) )                &
    3238               +      (           ibit32 * adv_mom_5                           &
    3239                      ) *                                                      &
    3240                                    ( w(k,j+3,i) - w(k,j-2,i) )                &
     3459                     ( 10.0_wp * ibit23 * adv_mom_5                           &
     3460                  +     3.0_wp * ibit22 * adv_mom_3                           &
     3461                  +              ibit21 * adv_mom_1                           &
     3462                     ) *                                                      &
     3463                                    ( w(k,j+1,i) - w(k,j,i)   )               &
     3464              -      (  5.0_wp * ibit23 * adv_mom_5                           &
     3465                  +              ibit22 * adv_mom_3                           &
     3466                     ) *                                                      &
     3467                                   ( w(k,j+2,i)  - w(k,j-1,i) )               &
     3468              +      (           ibit23 * adv_mom_5                           &
     3469                     ) *                                                      &
     3470                                   ( w(k,j+3,i)  - w(k,j-2,i) )               &
    32413471                                           )
    32423472       ENDDO
     
    32753505!--    calculated explicetely for the tendency at
    32763506!--    the first w-level. For topography wall this is done implicitely by
    3277 !--    advc_flags_1.
     3507!--    advc_flags_m.
    32783508       k         = nzb + 1
    32793509       w_comp(k) = w(k,j,i) + w(k-1,j,i)
     
    32853515!--       k index has to be modified near bottom and top, else array
    32863516!--       subscripts will be exceeded.
    3287           ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
    3288           ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
    3289           ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
    3290 
    3291           k_ppp = k + 3 * ibit35
    3292           k_pp  = k + 2 * ( 1 - ibit33  )
    3293           k_mm  = k - 2 * ibit35
     3517          ibit26 = REAL( IBITS(advc_flags_m(k,j,i),26,1), KIND = wp )
     3518          ibit25 = REAL( IBITS(advc_flags_m(k,j,i),25,1), KIND = wp )
     3519          ibit24 = REAL( IBITS(advc_flags_m(k,j,i),24,1), KIND = wp )
     3520
     3521          k_ppp = k + 3 * ibit26
     3522          k_pp  = k + 2 * ( 1 - ibit24  )
     3523          k_mm  = k - 2 * ibit26
    32943524
    32953525          w_comp(k) = w(k+1,j,i) + w(k,j,i)
    32963526          flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
    3297                      ( 37.0_wp * ibit35 * adv_mom_5                           &
    3298                   +     7.0_wp * ibit34 * adv_mom_3                           &
    3299                   +              ibit33 * adv_mom_1                           &
    3300                      ) *                                                      &
    3301                                 ( w(k+1,j,i)  + w(k,j,i)     )                &
    3302               -      (  8.0_wp * ibit35 * adv_mom_5                           &
    3303                   +              ibit34 * adv_mom_3                           &
     3527                     ( 37.0_wp * ibit26 * adv_mom_5                           &
     3528                  +     7.0_wp * ibit25 * adv_mom_3                           &
     3529                  +              ibit24 * adv_mom_1                           &
     3530                     ) *                                                      &
     3531                                ( w(k+1,j,i)   + w(k,j,i)    )                &
     3532              -      (  8.0_wp * ibit26 * adv_mom_5                           &
     3533                  +              ibit25 * adv_mom_3                           &
    33043534                     ) *                                                      &
    33053535                                ( w(k_pp,j,i)  + w(k-1,j,i)  )                &
    3306               +      (           ibit35 * adv_mom_5                           &
     3536              +      (           ibit26 * adv_mom_5                           &
    33073537                     ) *                                                      &
    33083538                                ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
     
    33103540
    33113541          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
    3312                      ( 10.0_wp * ibit35 * adv_mom_5                           &
    3313                   +     3.0_wp * ibit34 * adv_mom_3                           &
    3314                   +              ibit33 * adv_mom_1                           &
     3542                     ( 10.0_wp * ibit26 * adv_mom_5                           &
     3543                  +     3.0_wp * ibit25 * adv_mom_3                           &
     3544                  +              ibit24 * adv_mom_1                           &
    33153545                     ) *                                                      &
    33163546                                ( w(k+1,j,i)   - w(k,j,i)    )                &
    3317               -      (  5.0_wp * ibit35 * adv_mom_5                           &
    3318                   +              ibit34 * adv_mom_3                           &
     3547              -      (  5.0_wp * ibit26 * adv_mom_5                           &
     3548                  +              ibit25 * adv_mom_3                           &
    33193549                     ) *                                                      &
    33203550                                ( w(k_pp,j,i)  - w(k-1,j,i)  )                &
    3321               +      (           ibit35 * adv_mom_5                           &
     3551              +      (           ibit26 * adv_mom_5                           &
    33223552                     ) *                                                      &
    33233553                                ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
     
    33273557       DO  k = nzb+3, nzt-2
    33283558       
    3329           ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
    3330           ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
    3331           ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
     3559          ibit26 = REAL( IBITS(advc_flags_m(k,j,i),26,1), KIND = wp )
     3560          ibit25 = REAL( IBITS(advc_flags_m(k,j,i),25,1), KIND = wp )
     3561          ibit24 = REAL( IBITS(advc_flags_m(k,j,i),24,1), KIND = wp )
    33323562
    33333563          w_comp(k) = w(k+1,j,i) + w(k,j,i)
    33343564          flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
    3335                      ( 37.0_wp * ibit35 * adv_mom_5                           &
    3336                   +     7.0_wp * ibit34 * adv_mom_3                           &
    3337                   +              ibit33 * adv_mom_1                           &
    3338                      ) *                                                      &
    3339                                 ( w(k+1,j,i)  + w(k,j,i)     )                &
    3340               -      (  8.0_wp * ibit35 * adv_mom_5                           &
    3341                   +              ibit34 * adv_mom_3                           &
    3342                      ) *                                                      &
    3343                                 ( w(k+2,j,i)  + w(k-1,j,i)  )                 &
    3344               +      (           ibit35 * adv_mom_5                           &
    3345                      ) *                                                      &
    3346                                 ( w(k+3,j,i) + w(k-2,j,i) )                   &
     3565                     ( 37.0_wp * ibit26 * adv_mom_5                           &
     3566                  +     7.0_wp * ibit25 * adv_mom_3                           &
     3567                  +              ibit24 * adv_mom_1                           &
     3568                     ) *                                                      &
     3569                                ( w(k+1,j,i)  + w(k,j,i)                   &
     3570              -      (  8.0_wp * ibit26 * adv_mom_5                           &
     3571                  +              ibit25 * adv_mom_3                           &
     3572                     ) *                                                      &
     3573                                ( w(k+2,j,i)  + w(k-1,j,i) )                  &
     3574              +      (           ibit26 * adv_mom_5                           &
     3575                     ) *                                                      &
     3576                                ( w(k+3,j,i)  + w(k-2,j,i) )                  &
    33473577                                                 )
    33483578
    33493579          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
    3350                      ( 10.0_wp * ibit35 * adv_mom_5                           &
    3351                   +     3.0_wp * ibit34 * adv_mom_3                           &
    3352                   +              ibit33 * adv_mom_1                           &
    3353                      ) *                                                      &
    3354                                 ( w(k+1,j,i)   - w(k,j,i)    )                &
    3355               -      (  5.0_wp * ibit35 * adv_mom_5                           &
    3356                   +              ibit34 * adv_mom_3                           &
    3357                      ) *                                                      &
    3358                                 ( w(k+2,j,i)  - w(k-1,j,i)  )                 &
    3359               +      (           ibit35 * adv_mom_5                           &
    3360                      ) *                                                      &
    3361                                 ( w(k+3,j,i) - w(k-2,j,i) )                   &
     3580                     ( 10.0_wp * ibit26 * adv_mom_5                           &
     3581                  +     3.0_wp * ibit25 * adv_mom_3                           &
     3582                  +              ibit24 * adv_mom_1                           &
     3583                     ) *                                                      &
     3584                                ( w(k+1,j,i) - w(k,j,i)    )                  &
     3585              -      (  5.0_wp * ibit26 * adv_mom_5                           &
     3586                  +              ibit25 * adv_mom_3                           &
     3587                     ) *                                                      &
     3588                                ( w(k+2,j,i) - w(k-1,j,i)  )                  &
     3589              +      (           ibit26 * adv_mom_5                           &
     3590                     ) *                                                      &
     3591                                ( w(k+3,j,i) - w(k-2,j,i)  )                  &
    33623592                                                          )
    33633593       ENDDO
     
    33673597!--       k index has to be modified near bottom and top, else array
    33683598!--       subscripts will be exceeded.
    3369           ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
    3370           ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
    3371           ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
    3372 
    3373           k_ppp = k + 3 * ibit35
    3374           k_pp  = k + 2 * ( 1 - ibit33  )
    3375           k_mm  = k - 2 * ibit35
     3599          ibit26 = REAL( IBITS(advc_flags_m(k,j,i),26,1), KIND = wp )
     3600          ibit25 = REAL( IBITS(advc_flags_m(k,j,i),25,1), KIND = wp )
     3601          ibit24 = REAL( IBITS(advc_flags_m(k,j,i),24,1), KIND = wp )
     3602
     3603          k_ppp = k + 3 * ibit26
     3604          k_pp  = k + 2 * ( 1 - ibit24  )
     3605          k_mm  = k - 2 * ibit26
    33763606
    33773607          w_comp(k) = w(k+1,j,i) + w(k,j,i)
    33783608          flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
    3379                      ( 37.0_wp * ibit35 * adv_mom_5                           &
    3380                   +     7.0_wp * ibit34 * adv_mom_3                           &
    3381                   +              ibit33 * adv_mom_1                           &
    3382                      ) *                                                      &
    3383                                 ( w(k+1,j,i)  + w(k,j,i)     )                &
    3384               -      (  8.0_wp * ibit35 * adv_mom_5                           &
    3385                   +              ibit34 * adv_mom_3                           &
     3609                     ( 37.0_wp * ibit26 * adv_mom_5                           &
     3610                  +     7.0_wp * ibit25 * adv_mom_3                           &
     3611                  +              ibit24 * adv_mom_1                           &
     3612                     ) *                                                      &
     3613                                ( w(k+1,j,i)   + w(k,j,i)    )                &
     3614              -      (  8.0_wp * ibit26 * adv_mom_5                           &
     3615                  +              ibit25 * adv_mom_3                           &
    33863616                     ) *                                                      &
    33873617                                ( w(k_pp,j,i)  + w(k-1,j,i)  )                &
    3388               +      (           ibit35 * adv_mom_5                           &
     3618              +      (           ibit26 * adv_mom_5                           &
    33893619                     ) *                                                      &
    33903620                                ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
     
    33923622
    33933623          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
    3394                      ( 10.0_wp * ibit35 * adv_mom_5                           &
    3395                   +     3.0_wp * ibit34 * adv_mom_3                           &
    3396                   +              ibit33 * adv_mom_1                           &
     3624                     ( 10.0_wp * ibit26 * adv_mom_5                           &
     3625                  +     3.0_wp * ibit25 * adv_mom_3                           &
     3626                  +              ibit24 * adv_mom_1                           &
    33973627                     ) *                                                      &
    33983628                                ( w(k+1,j,i)   - w(k,j,i)    )                &
    3399               -      (  5.0_wp * ibit35 * adv_mom_5                           &
    3400                   +              ibit34 * adv_mom_3                           &
     3629              -      (  5.0_wp * ibit26 * adv_mom_5                           &
     3630                  +              ibit25 * adv_mom_3                           &
    34013631                     ) *                                                      &
    34023632                                ( w(k_pp,j,i)  - w(k-1,j,i)  )                &
    3403               +      (           ibit35 * adv_mom_5                           &
     3633              +      (           ibit26 * adv_mom_5                           &
    34043634                     ) *                                                      &
    34053635                                ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
     
    34123642          diss_d    = diss_t(k-1)
    34133643         
    3414           ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp )
    3415           ibit28 = REAL( IBITS(advc_flags_1(k,j,i),28,1), KIND = wp )
    3416           ibit27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp )
     3644          ibit20 = REAL( IBITS(advc_flags_m(k,j,i),20,1), KIND = wp )
     3645          ibit19 = REAL( IBITS(advc_flags_m(k,j,i),19,1), KIND = wp )
     3646          ibit18 = REAL( IBITS(advc_flags_m(k,j,i),18,1), KIND = wp )
    34173647         
    3418           ibit32 = REAL( IBITS(advc_flags_2(k,j,i),0,1), KIND = wp )
    3419           ibit31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp )
    3420           ibit30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp )
     3648          ibit23 = REAL( IBITS(advc_flags_m(k,j,i),23,1), KIND = wp )
     3649          ibit22 = REAL( IBITS(advc_flags_m(k,j,i),22,1), KIND = wp )
     3650          ibit21 = REAL( IBITS(advc_flags_m(k,j,i),21,1), KIND = wp )
    34213651         
    3422           ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
    3423           ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
    3424           ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
     3652          ibit26 = REAL( IBITS(advc_flags_m(k,j,i),26,1), KIND = wp )
     3653          ibit25 = REAL( IBITS(advc_flags_m(k,j,i),25,1), KIND = wp )
     3654          ibit24 = REAL( IBITS(advc_flags_m(k,j,i),24,1), KIND = wp )
    34253655!
    34263656!--       Calculate the divergence of the velocity field. A respective
    34273657!--       correction is needed to overcome numerical instabilities introduced
    34283658!--       by a not sufficient reduction of divergences near topography.
    3429           div = ( ( ( u_comp(k) + gu ) * ( ibit27 + ibit28 + ibit29 )         &
     3659          div = ( ( ( u_comp(k) + gu ) * ( ibit18 + ibit19 + ibit20 )         &
    34303660                  - ( u(k+1,j,i) + u(k,j,i)   )                               &
    34313661                                    * (                                       &
    3432                      REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp )     &
    3433                    + REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp )     &
    3434                    + REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )     &
     3662                     REAL( IBITS(advc_flags_m(k,j,i-1),18,1), KIND = wp )     &
     3663                   + REAL( IBITS(advc_flags_m(k,j,i-1),19,1), KIND = wp )     &
     3664                   + REAL( IBITS(advc_flags_m(k,j,i-1),20,1), KIND = wp )     &
    34353665                                      )                                       &
    34363666                  ) * ddx                                                     &
    3437               +   ( ( v_comp(k) + gv ) * ( ibit30 + ibit31 + ibit32 )         &
     3667              +   ( ( v_comp(k) + gv ) * ( ibit21 + ibit22 + ibit23 )         &
    34383668                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
    34393669                                    * (                                       &
    3440                      REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp )     &
    3441                    + REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )     &
    3442                    + REAL( IBITS(advc_flags_2(k,j-1,i),0,1), KIND = wp )     &
     3670                     REAL( IBITS(advc_flags_m(k,j-1,i),21,1), KIND = wp )     &
     3671                   + REAL( IBITS(advc_flags_m(k,j-1,i),22,1), KIND = wp )     &
     3672                   + REAL( IBITS(advc_flags_m(k,j-1,i),23,1), KIND = wp )     &
    34433673                                      )                                       &
    34443674                  ) * ddy                                                     &
    34453675              +   ( w_comp(k)               * rho_air(k+1)                    &
    3446                                             * ( ibit33 + ibit34 + ibit35 )    &
     3676                                            * ( ibit24 + ibit25 + ibit26 )    &
    34473677                - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k)                      &
    34483678                                    * (                                       &
    3449                      REAL( IBITS(advc_flags_2(k-1,j,i),1,1), KIND = wp )      &
    3450                    + REAL( IBITS(advc_flags_2(k-1,j,i),2,1), KIND = wp )      &
    3451                    + REAL( IBITS(advc_flags_2(k-1,j,i),3,1), KIND = wp )      &
     3679                     REAL( IBITS(advc_flags_m(k-1,j,i),24,1), KIND = wp )     &
     3680                   + REAL( IBITS(advc_flags_m(k-1,j,i),25,1), KIND = wp )     &
     3681                   + REAL( IBITS(advc_flags_m(k-1,j,i),26,1), KIND = wp )     &
    34523682                                      )                                       &
    34533683                  ) * drho_air_zw(k) * ddzu(k+1)                              &
     
    35313761!> Scalar advection - Call for all grid points
    35323762!------------------------------------------------------------------------------!
    3533     SUBROUTINE advec_s_ws( sk, sk_char )
     3763    SUBROUTINE advec_s_ws( advc_flag, sk, sk_char,                             &
     3764                           non_cyclic_l, non_cyclic_n,                         &
     3765                           non_cyclic_r, non_cyclic_s )
    35343766
    35353767
    35363768       CHARACTER (LEN = *), INTENT(IN)    ::  sk_char !< string identifier, used for assign fluxes to the correct dimension in the analysis array
    3537        INTEGER(iwp) ::  sk_num !< integer identifier, used for assign fluxes to the correct dimension in the analysis array
     3769       INTEGER(iwp) ::  sk_num = 0 !< integer identifier, used for assign fluxes to the correct dimension in the analysis array
    35383770       
    3539        INTEGER(iwp) ::  i         !< grid index along x-direction
    3540        INTEGER(iwp) ::  j         !< grid index along y-direction
    3541        INTEGER(iwp) ::  k         !< grid index along z-direction
    3542        INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
    3543        INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
    3544        INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
    3545        INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms
    3546        INTEGER(iwp) ::  tn = 0    !< number of OpenMP thread
     3771       INTEGER(iwp) ::  i          !< grid index along x-direction
     3772       INTEGER(iwp) ::  j          !< grid index along y-direction
     3773       INTEGER(iwp) ::  k          !< grid index along z-direction
     3774       INTEGER(iwp) ::  k_mm       !< k-2 index in disretization, can be modified to avoid segmentation faults
     3775       INTEGER(iwp) ::  k_pp       !< k+2 index in disretization, can be modified to avoid segmentation faults
     3776       INTEGER(iwp) ::  k_ppp      !< k+3 index in disretization, can be modified to avoid segmentation faults
     3777       INTEGER(iwp) ::  nzb_max_l  !< index indicating upper bound for order degradation of horizontal advection terms
     3778       INTEGER(iwp) ::  tn = 0     !< number of OpenMP thread
    35473779       
     3780       INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::   &
     3781                                                  advc_flag !< flag array to control order of scalar advection
     3782                                                 
     3783       LOGICAL ::  non_cyclic_l    !< flag that indicates non-cyclic boundary on the left
     3784       LOGICAL ::  non_cyclic_n    !< flag that indicates non-cyclic boundary on the north
     3785       LOGICAL ::  non_cyclic_r    !< flag that indicates non-cyclic boundary on the right
     3786       LOGICAL ::  non_cyclic_s    !< flag that indicates non-cyclic boundary on the south 
    35483787!       
    35493788!--    sk is an array from parameter list. It should not be a pointer, because
     
    35523791!--    even worse, because the compiler cannot assume strided one in the
    35533792!--    caller side.
    3554        REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !<  advected scalar
     3793       REAL(wp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !<  advected scalar
    35553794
    35563795       REAL(wp) ::  ibit0  !< flag indicating 1st-order scheme along x-direction
     
    36123851!--    entire subdomain, in order to avoid unsymmetric loops which might be
    36133852!--    an issue for GPUs.
    3614        IF( bc_dirichlet_l  .OR.  bc_radiation_l  .OR.                          &
    3615            bc_dirichlet_r  .OR.  bc_radiation_r  .OR.                          &
    3616            bc_dirichlet_s  .OR.  bc_radiation_s  .OR.                          &
    3617            bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
     3853       IF( non_cyclic_l  .OR.  non_cyclic_r  .OR.                             &
     3854           non_cyclic_s  .OR.  non_cyclic_n )  THEN
    36183855          nzb_max_l = nzt
    36193856       ELSE
     
    36523889          DO  k = nzb+1, nzb_max_l
    36533890
    3654              ibit2 = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )
    3655              ibit1 = REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )
    3656              ibit0 = REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )
     3891             ibit2 = REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp )
     3892             ibit1 = REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp )
     3893             ibit0 = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )
    36573894
    36583895             u_comp                 = u(k,j,i) - u_gtrans + u_stokes_zu(k)
     
    37193956       !$ACC PRIVATE(flux_t, diss_t, flux_d, diss_d) &
    37203957       !$ACC PRIVATE(div, u_comp, u_comp_l, v_comp, v_comp_s) &
    3721        !$ACC PRESENT(advc_flags_1) &
     3958       !$ACC PRESENT(advc_flag) &
    37223959       !$ACC PRESENT(sk, u, v, w, u_stokes_zu, v_stokes_zu) &
    37233960       !$ACC PRESENT(drho_air, rho_air_zw, ddzw) &
     
    37363973          DO  k = nzb+1, nzb_max_l
    37373974
    3738              ibit5 = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )
    3739              ibit4 = REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )
    3740              ibit3 = REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )
     3975             ibit5 = REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp )
     3976             ibit4 = REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp )
     3977             ibit3 = REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp )
    37413978
    37423979             v_comp               = v(k,j,i) - v_gtrans + v_stokes_zu(k)
     
    37984035             DO  k = nzb+1, nzb_max_l
    37994036
    3800                 ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp )
    3801                 ibit1 = REAL( IBITS(advc_flags_1(k,j,i),1,1), KIND = wp )
    3802                 ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp )
     4037                ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp )
     4038                ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp )
     4039                ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp )
    38034040
    38044041                u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
     
    38364073!
    38374074!--             Recompute the left fluxes.
    3838                 ibit2_l = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )
    3839                 ibit1_l = REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )
    3840                 ibit0_l = REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )
     4075                ibit2_l = REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp )
     4076                ibit1_l = REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp )
     4077                ibit0_l = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )
    38414078
    38424079                u_comp_l = u(k,j,i) - u_gtrans + u_stokes_zu(k)
     
    38754112#endif
    38764113
    3877                 ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp )
    3878                 ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp )
    3879                 ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp )
     4114                ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp )
     4115                ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp )
     4116                ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp )
    38804117
    38814118                v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
     
    39134150!
    39144151!--             Recompute the south fluxes.
    3915                 ibit5_s = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )
    3916                 ibit4_s = REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )
    3917                 ibit3_s = REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )
     4152                ibit5_s = REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp )
     4153                ibit4_s = REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp )
     4154                ibit3_s = REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp )
    39184155
    39194156                v_comp_s = v(k,j,i) - v_gtrans + v_stokes_zu(k)
     
    39554192!--             k index has to be modified near bottom and top, else array
    39564193!--             subscripts will be exceeded.
    3957                 ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
    3958                 ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
    3959                 ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
     4194                ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     4195                ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     4196                ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
    39604197
    39614198                k_ppp = k + 3 * ibit8
     
    39984235                div   =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
    39994236                          - u(k,j,i)   * (                                     &
    4000                         REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )    &
    4001                       + REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )    &
    4002                       + REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )    &
     4237                        REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp )       &
     4238                      + REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp )       &
     4239                      + REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp )       &
    40034240                                         )                                     &
    40044241                          ) * ddx                                              &
    40054242                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
    40064243                          - v(k,j,i)   * (                                     &
    4007                         REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )    &
    4008                       + REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )    &
    4009                       + REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )    &
     4244                        REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp )       &
     4245                      + REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp )       &
     4246                      + REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp )       &
    40104247                                         )                                     &
    40114248                          ) * ddy                                              &
     
    40144251                          - w(k-1,j,i) * rho_air_zw(k-1) *                     &
    40154252                                         (                                     &
    4016                         REAL( IBITS(advc_flags_1(k-1,j,i),6,1), KIND = wp )    &
    4017                       + REAL( IBITS(advc_flags_1(k-1,j,i),7,1), KIND = wp )    &
    4018                       + REAL( IBITS(advc_flags_1(k-1,j,i),8,1), KIND = wp )    &
     4253                        REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp )       &
     4254                      + REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp )       &
     4255                      + REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp )       &
    40194256                                         )                                     &     
    40204257                          ) * drho_air(k) * ddzw(k)
     
    42024439!--             k index has to be modified near bottom and top, else array
    42034440!--             subscripts will be exceeded.
    4204                 ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
    4205                 ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
    4206                 ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
     4441                ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
     4442                ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     4443                ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp )
    42074444
    42084445                k_ppp = k + 3 * ibit8
     
    43894626       INTEGER(iwp) ::  tn = 0    !< number of OpenMP thread
    43904627       
    4391        REAL(wp)    ::  ibit9 !< flag indicating 1st-order scheme along x-direction
    4392        REAL(wp)    ::  ibit10 !< flag indicating 3rd-order scheme along x-direction
    4393        REAL(wp)    ::  ibit11 !< flag indicating 5th-order scheme along x-direction
     4628       REAL(wp)    ::  ibit0 !< flag indicating 1st-order scheme along x-direction
     4629       REAL(wp)    ::  ibit1 !< flag indicating 3rd-order scheme along x-direction
     4630       REAL(wp)    ::  ibit2 !< flag indicating 5th-order scheme along x-direction
    43944631#ifdef _OPENACC
    4395        REAL(wp)    ::  ibit9_l !< flag indicating 1st-order scheme along x-direction
    4396        REAL(wp)    ::  ibit10_l !< flag indicating 3rd-order scheme along x-direction
    4397        REAL(wp)    ::  ibit11_l !< flag indicating 5th-order scheme along x-direction
     4632       REAL(wp)    ::  ibit0_l !< flag indicating 1st-order scheme along x-direction
     4633       REAL(wp)    ::  ibit1_l !< flag indicating 3rd-order scheme along x-direction
     4634       REAL(wp)    ::  ibit2_l !< flag indicating 5th-order scheme along x-direction
    43984635#endif
    4399        REAL(wp)    ::  ibit12 !< flag indicating 1st-order scheme along y-direction
    4400        REAL(wp)    ::  ibit13 !< flag indicating 3rd-order scheme along y-direction
    4401        REAL(wp)    ::  ibit14 !< flag indicating 5th-order scheme along y-direction
     4636       REAL(wp)    ::  ibit3 !< flag indicating 1st-order scheme along y-direction
     4637       REAL(wp)    ::  ibit4 !< flag indicating 3rd-order scheme along y-direction
     4638       REAL(wp)    ::  ibit5 !< flag indicating 5th-order scheme along y-direction
    44024639#ifdef _OPENACC
    4403        REAL(wp)    ::  ibit12_s !< flag indicating 1st-order scheme along y-direction
    4404        REAL(wp)    ::  ibit13_s !< flag indicating 3rd-order scheme along y-direction
    4405        REAL(wp)    ::  ibit14_s !< flag indicating 5th-order scheme along y-direction
     4640       REAL(wp)    ::  ibit3_s !< flag indicating 1st-order scheme along y-direction
     4641       REAL(wp)    ::  ibit4_s !< flag indicating 3rd-order scheme along y-direction
     4642       REAL(wp)    ::  ibit5_s !< flag indicating 5th-order scheme along y-direction
    44064643#endif
    4407        REAL(wp)    ::  ibit15 !< flag indicating 1st-order scheme along z-direction
    4408        REAL(wp)    ::  ibit16 !< flag indicating 3rd-order scheme along z-direction
    4409        REAL(wp)    ::  ibit17 !< flag indicating 5th-order scheme along z-direction
     4644       REAL(wp)    ::  ibit6 !< flag indicating 1st-order scheme along z-direction
     4645       REAL(wp)    ::  ibit7 !< flag indicating 3rd-order scheme along z-direction
     4646       REAL(wp)    ::  ibit8 !< flag indicating 5th-order scheme along z-direction
    44104647       REAL(wp)    ::  diss_d !< artificial dissipation term at grid box bottom
    44114648       REAL(wp)    ::  div    !< diverence on u-grid
     
    44684705          DO  k = nzb+1, nzb_max_l
    44694706
    4470              ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )
    4471              ibit10 = REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )
    4472              ibit9  = REAL( IBITS(advc_flags_1(k,j,i-1),9,1), KIND = wp )
     4707             ibit2 = REAL( IBITS(advc_flags_m(k,j,i-1),2,1), KIND = wp )
     4708             ibit1 = REAL( IBITS(advc_flags_m(k,j,i-1),1,1), KIND = wp )
     4709             ibit0 = REAL( IBITS(advc_flags_m(k,j,i-1),0,1), KIND = wp )
    44734710
    44744711             u_comp                   = u(k,j,i) + u(k,j,i-1) - gu
    44754712             swap_flux_x_local_u(k,j) = u_comp * (                             &
    4476                                        ( 37.0_wp * ibit11 * adv_mom_5          &
    4477                                     +     7.0_wp * ibit10 * adv_mom_3          &
    4478                                     +              ibit9 * adv_mom_1          &
     4713                                       ( 37.0_wp * ibit2 * adv_mom_5           &
     4714                                    +     7.0_wp * ibit1 * adv_mom_3           &
     4715                                    +              ibit0 * adv_mom_1          &
    44794716                                       ) *                                     &
    44804717                                     ( u(k,j,i)   + u(k,j,i-1) )               &
    4481                                 -      (  8.0_wp * ibit11 * adv_mom_5          &
    4482                                     +              ibit10 * adv_mom_3          &
     4718                                -      (  8.0_wp * ibit2 * adv_mom_5           &
     4719                                    +              ibit1 * adv_mom_3           &