Index: /palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- /palm/trunk/SOURCE/advec_ws.f90	(revision 3550)
+++ /palm/trunk/SOURCE/advec_ws.f90	(revision 3551)
@@ -25,4 +25,12 @@
 ! -----------------
 ! $Id$
+! - Computation of vertical fluxes separated from computation of horizontal 
+!   fluxes. Loops are splitted in order to avoid indirect indexing and allow 
+!   for better vectorization. 
+! - Accelerate code by remove type-conversions of ibits
+! - replace pointer definition in scalar routine by simple explicit definition
+!   of the passed array
+! 
+! 3547 2018-11-21 13:21:24Z suehring
 ! variables documented
 ! 
@@ -193,5 +201,5 @@
 !
 ! 888 2012-04-20 15:03:46Z suehring
-! Number of IBITS() calls with identical arguments is reduced.
+! Number of REAL( IBITS() calls with identical arguments is reduced.
 !
 ! 862 2012-03-26 14:21:38Z suehring
@@ -1200,13 +1208,4 @@
        
        INTEGER(iwp) ::  i     !< grid index along x-direction
-       INTEGER(iwp) ::  ibit0 !< flag indicating 1st-order scheme along x-direction
-       INTEGER(iwp) ::  ibit1 !< flag indicating 3rd-order scheme along x-direction
-       INTEGER(iwp) ::  ibit2 !< flag indicating 5th-order scheme along x-direction
-       INTEGER(iwp) ::  ibit3 !< flag indicating 1st-order scheme along y-direction
-       INTEGER(iwp) ::  ibit4 !< flag indicating 3rd-order scheme along y-direction
-       INTEGER(iwp) ::  ibit5 !< flag indicating 5th-order scheme along y-direction
-       INTEGER(iwp) ::  ibit6 !< flag indicating 1st-order scheme along z-direction
-       INTEGER(iwp) ::  ibit7 !< flag indicating 3rd-order scheme along z-direction
-       INTEGER(iwp) ::  ibit8 !< flag indicating 5th-order scheme along z-direction
        INTEGER(iwp) ::  i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread
        INTEGER(iwp) ::  j     !< grid index along y-direction
@@ -1217,4 +1216,13 @@
        INTEGER(iwp) ::  tn    !< number of OpenMP thread
        
+       REAL(wp)     ::  ibit0  !< flag indicating 1st-order scheme along x-direction
+       REAL(wp)     ::  ibit1  !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp)     ::  ibit2  !< flag indicating 5th-order scheme along x-direction
+       REAL(wp)     ::  ibit3  !< flag indicating 1st-order scheme along y-direction
+       REAL(wp)     ::  ibit4  !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp)     ::  ibit5  !< flag indicating 5th-order scheme along y-direction
+       REAL(wp)     ::  ibit6  !< flag indicating 1st-order scheme along z-direction
+       REAL(wp)     ::  ibit7  !< flag indicating 3rd-order scheme along z-direction
+       REAL(wp)     ::  ibit8  !< flag indicating 5th-order scheme along z-direction
        REAL(wp)     ::  diss_d !< artificial dissipation term at grid box bottom
        REAL(wp)     ::  div    !< diverence on scalar grid
@@ -1222,10 +1230,12 @@
        REAL(wp)     ::  u_comp !< advection velocity along x-direction 
        REAL(wp)     ::  v_comp !< advection velocity along y-direction 
-       
-#if defined( __nopointer )
-       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !< advected scalar
-#else
-       REAL(wp), DIMENSION(:,:,:), POINTER    ::  sk     !< advected scalar
-#endif
+!       
+!--    sk is an array from parameter list. It should not be a pointer, because 
+!--    in that case the compiler can not assume a stride 1 and cannot perform 
+!--    a strided one vector load. Adding the CONTIGUOUS keyword makes things 
+!--    even worse, because the compiler cannot assume strided one in the 
+!--    caller side.
+       REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !<  advected scalar
+
        REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
        REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
@@ -1249,7 +1259,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit5 = IBITS(advc_flags_1(k,j-1,i),5,1)
-             ibit4 = IBITS(advc_flags_1(k,j-1,i),4,1)
-             ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
+             ibit5 = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )
+             ibit4 = REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )
+             ibit3 = REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )
 
              v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
@@ -1310,7 +1320,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit2 = IBITS(advc_flags_1(k,j,i-1),2,1)
-             ibit1 = IBITS(advc_flags_1(k,j,i-1),1,1)
-             ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
+             ibit2 = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )
+             ibit1 = REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )
+             ibit0 = REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )
 
              u_comp                     = u(k,j,i) - u_gtrans + u_stokes_zu(k)
@@ -1365,7 +1375,4 @@
            
        ENDIF
-
-       flux_t(0) = 0.0_wp
-       diss_t(0) = 0.0_wp
 !        
 !--    Now compute the fluxes and tendency terms for the horizontal and
@@ -1377,7 +1384,7 @@
 !--       flux at the end. 
 
-          ibit2 = IBITS(advc_flags_1(k,j,i),2,1)
-          ibit1 = IBITS(advc_flags_1(k,j,i),1,1)
-          ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
+          ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp )
+          ibit1 = REAL( IBITS(advc_flags_1(k,j,i),1,1), KIND = wp )
+          ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp )
 
           u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
@@ -1412,7 +1419,7 @@
                                        )
 
-          ibit5 = IBITS(advc_flags_1(k,j,i),5,1)
-          ibit4 = IBITS(advc_flags_1(k,j,i),4,1)
-          ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
+          ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp )
+          ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp )
+          ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp )
 
           v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
@@ -1446,91 +1453,4 @@
                              ( sk(k,j+3,i) - sk(k,j-2,i) )                    &
                                        )
-!
-!--       k index has to be modified near bottom and top, else array
-!--       subscripts will be exceeded.
-          ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
-          ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
-          ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
-
-          k_ppp = k + 3 * ibit8
-          k_pp  = k + 2 * ( 1 - ibit6  )
-          k_mm  = k - 2 * ibit8
-
-
-          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
-                     ( 37.0_wp * ibit8 * adv_sca_5                            &
-                  +     7.0_wp * ibit7 * adv_sca_3                            &
-                  +              ibit6 * adv_sca_1                            &
-                     ) *                                                      &
-                             ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
-              -      (  8.0_wp * ibit8 * adv_sca_5                            &
-                  +              ibit7 * adv_sca_3                            &
-                     ) *                                                      &
-                             ( sk(k_pp,j,i) + sk(k-1,j,i)  )                  &
-              +      (           ibit8 * adv_sca_5                            &
-                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                  &
-                                 )
-
-          diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
-                     ( 10.0_wp * ibit8 * adv_sca_5                            &
-                  +     3.0_wp * ibit7 * adv_sca_3                            &
-                  +              ibit6 * adv_sca_1                            &
-                     ) *                                                      &
-                             ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
-              -      (  5.0_wp * ibit8 * adv_sca_5                            &
-                  +              ibit7 * adv_sca_3                            &
-                     ) *                                                      &
-                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )                 &
-              +      (           ibit8 * adv_sca_5                            &
-                     ) *                                                      &
-                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
-                                         )
-       ENDDO
-
-       DO  k = nzb+1, nzb_max
-
-          flux_d    = flux_t(k-1)
-          diss_d    = diss_t(k-1)
-!
-!--       Calculate the divergence of the velocity field. A respective
-!--       correction is needed to overcome numerical instabilities caused
-!--       by a not sufficient reduction of divergences near topography. 
-          div         =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
-                          - u(k,j,i)   * ( IBITS(advc_flags_1(k,j,i-1),0,1)    &
-                                         + IBITS(advc_flags_1(k,j,i-1),1,1)    &
-                                         + IBITS(advc_flags_1(k,j,i-1),2,1)    &
-                                         )                                     &
-                          ) * ddx                                              &
-                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
-                          - v(k,j,i)   * ( IBITS(advc_flags_1(k,j-1,i),3,1)    &
-                                         + IBITS(advc_flags_1(k,j-1,i),4,1)    &
-                                         + IBITS(advc_flags_1(k,j-1,i),5,1)    &
-                                         )                                     &
-                          ) * ddy                                              &
-                        + ( w(k,j,i) * rho_air_zw(k) *                         &
-                                         ( ibit6 + ibit7 + ibit8 )             &
-                          - w(k-1,j,i) * rho_air_zw(k-1) *                     &
-                                         ( IBITS(advc_flags_1(k-1,j,i),6,1)    &
-                                         + IBITS(advc_flags_1(k-1,j,i),7,1)    &
-                                         + IBITS(advc_flags_1(k-1,j,i),8,1)    &
-                                         )                                     &      
-                          ) * drho_air(k) * ddzw(k)
-
-
-          tend(k,j,i) = tend(k,j,i) - (                                       &
-                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
-                          swap_diss_x_local(k,j,tn)            ) * ddx        &
-                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
-                          swap_diss_y_local(k,tn)              ) * ddy        &
-                      + ( ( flux_t(k) + diss_t(k) ) -                         &
-                          ( flux_d    + diss_d    )                           &
-                                                    ) * drho_air(k) * ddzw(k) &
-                                      ) + sk(k,j,i) * div
-
-          swap_flux_y_local(k,tn)   = flux_n(k)
-          swap_diss_y_local(k,tn)   = diss_n(k)
-          swap_flux_x_local(k,j,tn) = flux_r(k)
-          swap_diss_x_local(k,j,tn) = diss_r(k)
-
        ENDDO
 !
@@ -1559,47 +1479,137 @@
                     -  5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
                     +           ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
+
+       ENDDO
+!
+!--    Now, compute vertical fluxes. Split loop into a part treating the 
+!--    lowest 2 grid points with indirect indexing, a main loop without
+!--    indirect indexing, and a loop for the uppermost 2 grip points with
+!--    indirect indexing. This allows better vectorization for the main loop.
+!--    First, compute the flux at model surface, which need has to be 
+!--    calculated explicetely for the tendency at
+!--    the first w-level. For topography wall this is done implicitely by
+!--    advc_flags_1.
+       flux_t(nzb) = 0.0_wp
+       diss_t(nzb) = 0.0_wp
+       DO  k = nzb+1, nzb+2
+          ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
+          ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
+          ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
 !
 !--       k index has to be modified near bottom and top, else array
 !--       subscripts will be exceeded.
-          ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
-          ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
-          ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
-
           k_ppp = k + 3 * ibit8
           k_pp  = k + 2 * ( 1 - ibit6  )
           k_mm  = k - 2 * ibit8
 
+
           flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
-                    ( 37.0_wp * ibit8 * adv_sca_5                             &
-                 +     7.0_wp * ibit7 * adv_sca_3                             &
-                 +              ibit6 * adv_sca_1                             &
-                    ) *                                                       &
-                             ( sk(k+1,j,i)  + sk(k,j,i)   )                   &
-              -     (  8.0_wp * ibit8 * adv_sca_5                             &
+                     ( 37.0_wp * ibit8 * adv_sca_5                            &
+                  +     7.0_wp * ibit7 * adv_sca_3                            &
+                  +              ibit6 * adv_sca_1                            &
+                     ) *                                                      &
+                             ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
+              -      (  8.0_wp * ibit8 * adv_sca_5                            &
                   +              ibit7 * adv_sca_3                            &
-                    ) *                                                       &
-                             ( sk(k_pp,j,i) + sk(k-1,j,i) )                   &
-              +     (           ibit8 * adv_sca_5                             &
-                    ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                   &
+                     ) *                                                      &
+                             ( sk(k_pp,j,i) + sk(k-1,j,i)  )                  &
+              +      (           ibit8 * adv_sca_5                            &
+                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                  &
                                  )
 
           diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
-                    ( 10.0_wp * ibit8 * adv_sca_5                             &
-                 +     3.0_wp * ibit7 * adv_sca_3                             &
-                 +              ibit6 * adv_sca_1                             &
-                    ) *                                                       &
+                     ( 10.0_wp * ibit8 * adv_sca_5                            &
+                  +     3.0_wp * ibit7 * adv_sca_3                            &
+                  +              ibit6 * adv_sca_1                            &
+                     ) *                                                      &
                              ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
-              -     (  5.0_wp * ibit8 * adv_sca_5                             &
-                 +              ibit7 * adv_sca_3                             &
-                    ) *                                                       &
+              -      (  5.0_wp * ibit8 * adv_sca_5                            &
+                  +              ibit7 * adv_sca_3                            &
+                     ) *                                                      &
                              ( sk(k_pp,j,i)  - sk(k-1,j,i)  )                 &
-              +     (           ibit8 * adv_sca_5                             &
-                    ) *                                                       &
+              +      (           ibit8 * adv_sca_5                            &
+                     ) *                                                      &
                              ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
                                          )
-
-       ENDDO
-
-       DO  k = nzb_max+1, nzt
+       ENDDO 
+       
+       DO  k = nzb+3, nzt-2
+          ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
+          ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
+          ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
+
+          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
+                     ( 37.0_wp * ibit8 * adv_sca_5                            &
+                  +     7.0_wp * ibit7 * adv_sca_3                            &
+                  +              ibit6 * adv_sca_1                            &
+                     ) *                                                      &
+                             ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
+              -      (  8.0_wp * ibit8 * adv_sca_5                            &
+                  +              ibit7 * adv_sca_3                            &
+                     ) *                                                      &
+                             ( sk(k+2,j,i) + sk(k-1,j,i)  )                   &
+              +      (           ibit8 * adv_sca_5                            &
+                     ) *     ( sk(k+3,j,i)+ sk(k-2,j,i) )                     &
+                                                 )
+
+          diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
+                     ( 10.0_wp * ibit8 * adv_sca_5                            &
+                  +     3.0_wp * ibit7 * adv_sca_3                            &
+                  +              ibit6 * adv_sca_1                            &
+                     ) *                                                      &
+                             ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
+              -      (  5.0_wp * ibit8 * adv_sca_5                            &
+                  +              ibit7 * adv_sca_3                            &
+                     ) *                                                      &
+                             ( sk(k+2,j,i)  - sk(k-1,j,i)  )                  &
+              +      (           ibit8 * adv_sca_5                            &
+                     ) *                                                      &
+                             ( sk(k+3,j,i) - sk(k-2,j,i) )                    &
+                                                         )
+       ENDDO       
+       
+       DO  k = nzt-1, nzt
+          ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
+          ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
+          ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
+!
+!--       k index has to be modified near bottom and top, else array
+!--       subscripts will be exceeded.
+          k_ppp = k + 3 * ibit8
+          k_pp  = k + 2 * ( 1 - ibit6  )
+          k_mm  = k - 2 * ibit8
+
+
+          flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                            &
+                     ( 37.0_wp * ibit8 * adv_sca_5                            &
+                  +     7.0_wp * ibit7 * adv_sca_3                            &
+                  +              ibit6 * adv_sca_1                            &
+                     ) *                                                      &
+                             ( sk(k+1,j,i)  + sk(k,j,i)    )                  &
+              -      (  8.0_wp * ibit8 * adv_sca_5                            &
+                  +              ibit7 * adv_sca_3                            &
+                     ) *                                                      &
+                             ( sk(k_pp,j,i) + sk(k-1,j,i)  )                  &
+              +      (           ibit8 * adv_sca_5                            &
+                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                  &
+                                                 )
+
+          diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                    &
+                     ( 10.0_wp * ibit8 * adv_sca_5                            &
+                  +     3.0_wp * ibit7 * adv_sca_3                            &
+                  +              ibit6 * adv_sca_1                            &
+                     ) *                                                      &
+                             ( sk(k+1,j,i)   - sk(k,j,i)    )                 &
+              -      (  5.0_wp * ibit8 * adv_sca_5                            &
+                  +              ibit7 * adv_sca_3                            &
+                     ) *                                                      &
+                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )                 &
+              +      (           ibit8 * adv_sca_5                            &
+                     ) *                                                      &
+                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
+                                                         )
+       ENDDO 
+
+       DO  k = nzb+1, nzt
 
           flux_d    = flux_t(k-1)
@@ -1610,8 +1620,26 @@
 !--       correction is needed to overcome numerical instabilities introduced
 !--       by a not sufficient reduction of divergences near topography.
-          div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx                   &
-                        + ( v(k,j+1,i) - v(k,j,i)   ) * ddy                   &
-                        + ( w(k,j,i)   * rho_air_zw(k) -                      &
-                            w(k-1,j,i) * rho_air_zw(k-1)                      &
+          div         =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )            &
+                          - u(k,j,i)   * (                                    &
+                        REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )   &
+                      + REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )   &
+                      + REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )   &
+                                         )                                    &
+                          ) * ddx                                             &
+                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )            &
+                          - v(k,j,i)   * (                                    &
+                        REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )   &
+                      + REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )   &
+                      + REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )   &
+                                         )                                    &
+                          ) * ddy                                             &
+                        + ( w(k,j,i) * rho_air_zw(k) *                        &
+                                         ( ibit6 + ibit7 + ibit8 )            &
+                          - w(k-1,j,i) * rho_air_zw(k-1) *                    &
+                                         (                                    &
+                        REAL( IBITS(advc_flags_1(k-1,j,i),6,1), KIND = wp )   &
+                      + REAL( IBITS(advc_flags_1(k-1,j,i),7,1), KIND = wp )   &
+                      + REAL( IBITS(advc_flags_1(k-1,j,i),8,1), KIND = wp )   &
+                                         )                                    &      
                           ) * drho_air(k) * ddzw(k)
 
@@ -1777,13 +1805,4 @@
 
        INTEGER(iwp) ::  i      !< grid index along x-direction
-       INTEGER(iwp) ::  ibit9  !< flag indicating 1st-order scheme along x-direction
-       INTEGER(iwp) ::  ibit10 !< flag indicating 3rd-order scheme along x-direction
-       INTEGER(iwp) ::  ibit11 !< flag indicating 5th-order scheme along x-direction
-       INTEGER(iwp) ::  ibit12 !< flag indicating 1st-order scheme along y-direction
-       INTEGER(iwp) ::  ibit13 !< flag indicating 3rd-order scheme along y-direction
-       INTEGER(iwp) ::  ibit14 !< flag indicating 5th-order scheme along y-direction
-       INTEGER(iwp) ::  ibit15 !< flag indicating 1st-order scheme along z-direction
-       INTEGER(iwp) ::  ibit16 !< flag indicating 3rd-order scheme along z-direction
-       INTEGER(iwp) ::  ibit17 !< flag indicating 5th-order scheme along z-direction
        INTEGER(iwp) ::  i_omp  !< leftmost index on subdomain, or in case of OpenMP, on thread
        INTEGER(iwp) ::  j      !< grid index along y-direction
@@ -1794,4 +1813,13 @@
        INTEGER(iwp) ::  tn     !< number of OpenMP thread
        
+       REAL(wp)    ::  ibit9    !< flag indicating 1st-order scheme along x-direction
+       REAL(wp)    ::  ibit10   !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp)    ::  ibit11   !< flag indicating 5th-order scheme along x-direction
+       REAL(wp)    ::  ibit12   !< flag indicating 1st-order scheme along y-direction
+       REAL(wp)    ::  ibit13   !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp)    ::  ibit14   !< flag indicating 5th-order scheme along y-direction
+       REAL(wp)    ::  ibit15   !< flag indicating 1st-order scheme along z-direction
+       REAL(wp)    ::  ibit16   !< flag indicating 3rd-order scheme along z-direction
+       REAL(wp)    ::  ibit17   !< flag indicating 5th-order scheme along z-direction
        REAL(wp)    ::  diss_d   !< artificial dissipation term at grid box bottom
        REAL(wp)    ::  div      !< diverence on u-grid
@@ -1819,7 +1847,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit14 = IBITS(advc_flags_1(k,j-1,i),14,1)
-             ibit13 = IBITS(advc_flags_1(k,j-1,i),13,1)
-             ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)
+             ibit14 = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )
+             ibit13 = REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp )
+             ibit12 = REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp )
 
              v_comp(k)      = v(k,j,i) + v(k,j,i-1) - gv
@@ -1877,7 +1905,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit11 = IBITS(advc_flags_1(k,j,i-1),11,1)
-             ibit10 = IBITS(advc_flags_1(k,j,i-1),10,1)
-             ibit9  = IBITS(advc_flags_1(k,j,i-1),9,1)
+             ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )
+             ibit10 = REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )
+             ibit9  = REAL( IBITS(advc_flags_1(k,j,i-1),9,1),  KIND = wp )
 
              u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
@@ -1929,8 +1957,4 @@
           
        ENDIF
-
-       flux_t(0) = 0.0_wp
-       diss_t(0) = 0.0_wp
-       w_comp(0) = 0.0_wp
 !
 !--    Now compute the fluxes tendency terms for the horizontal and
@@ -1938,7 +1962,7 @@
        DO  k = nzb+1, nzb_max
 
-          ibit11 = IBITS(advc_flags_1(k,j,i),11,1)
-          ibit10 = IBITS(advc_flags_1(k,j,i),10,1)
-          ibit9  = IBITS(advc_flags_1(k,j,i),9,1)
+          ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp )
+          ibit10 = REAL( IBITS(advc_flags_1(k,j,i),10,1), KIND = wp )
+          ibit9  = REAL( IBITS(advc_flags_1(k,j,i),9,1),  KIND = wp )
 
           u_comp(k) = u(k,j,i+1) + u(k,j,i)
@@ -1973,7 +1997,7 @@
                                                 )
 
-          ibit14 = IBITS(advc_flags_1(k,j,i),14,1)
-          ibit13 = IBITS(advc_flags_1(k,j,i),13,1)
-          ibit12 = IBITS(advc_flags_1(k,j,i),12,1)
+          ibit14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp )
+          ibit13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp )
+          ibit12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp )
 
           v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
@@ -2007,10 +2031,48 @@
                                     ( u(k,j+3,i) - u(k,j-2,i) )               &
                                             )
+       ENDDO
+
+       DO  k = nzb_max+1, nzt
+
+          u_comp(k) = u(k,j,i+1) + u(k,j,i)
+          flux_r(k) = ( u_comp(k) - gu ) * (                                  &
+                         37.0_wp * ( u(k,j,i+1) + u(k,j,i)   )                &
+                       -  8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) )                &
+                       +           ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
+          diss_r(k) = - ABS( u_comp(k) - gu ) * (                             &
+                         10.0_wp * ( u(k,j,i+1) - u(k,j,i)   )                &
+                       -  5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) )                &
+                       +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
+
+          v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
+          flux_n(k) = v_comp(k) * (                                           &
+                         37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                &
+                       -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                &
+                       +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
+          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
+                         10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                &
+                       -  5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) )                &
+                       +           ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
+
+       ENDDO
+!
+!--    Now, compute vertical fluxes. Split loop into a part treating the 
+!--    lowest 2 grid points with indirect indexing, a main loop without
+!--    indirect indexing, and a loop for the uppermost 2 grip points with
+!--    indirect indexing. This allows better vectorization for the main loop.
+!--    First, compute the flux at model surface, which need has to be 
+!--    calculated explicetely for the tendency at
+!--    the first w-level. For topography wall this is done implicitely by
+!--    advc_flags_1.
+       flux_t(nzb) = 0.0_wp
+       diss_t(nzb) = 0.0_wp
+       w_comp(nzb) = 0.0_wp
+       DO  k = nzb+1, nzb+2
 !
 !--       k index has to be modified near bottom and top, else array
 !--       subscripts will be exceeded.
-          ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
-          ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
-          ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
+          ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
+          ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
+          ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
 
           k_ppp = k + 3 * ibit17
@@ -2049,6 +2111,88 @@
                                                            )
        ENDDO
-
-       DO  k = nzb+1, nzb_max
+       
+       DO  k = nzb+3, nzt-2
+
+          ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
+          ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
+          ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
+
+          w_comp(k) = w(k,j,i) + w(k,j,i-1)
+          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
+                     ( 37.0_wp * ibit17 * adv_mom_5                           &
+                  +     7.0_wp * ibit16 * adv_mom_3                           &
+                  +              ibit15 * adv_mom_1                           &
+                     ) *                                                      &
+                                ( u(k+1,j,i)  + u(k,j,i)     )                &
+              -      (  8.0_wp * ibit17 * adv_mom_5                           &
+                  +              ibit16 * adv_mom_3                           &
+                     ) *                                                      &
+                                ( u(k+2,j,i) + u(k-1,j,i)   )                 &
+              +      (           ibit17 * adv_mom_5                           &
+                     ) *                                                      &
+                                ( u(k+3,j,i) + u(k-2,j,i) )                   &
+                                                  )
+
+          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
+                     ( 10.0_wp * ibit17 * adv_mom_5                           &
+                  +     3.0_wp * ibit16 * adv_mom_3                           &
+                  +              ibit15 * adv_mom_1                           &
+                     ) *                                                      &
+                                ( u(k+1,j,i)   - u(k,j,i)    )                &
+              -      (  5.0_wp * ibit17 * adv_mom_5                           &
+                  +              ibit16 * adv_mom_3                           &
+                     ) *                                                      &
+                                ( u(k+2,j,i)  - u(k-1,j,i)  )                 &
+              +      (           ibit17 * adv_mom_5                           &
+                     ) *                                                      &
+                                ( u(k+3,j,i) - u(k-2,j,i) )                   &
+                                                           )
+       ENDDO
+       
+       DO  k = nzt-1, nzt
+!
+!--       k index has to be modified near bottom and top, else array
+!--       subscripts will be exceeded.
+          ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
+          ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
+          ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
+
+          k_ppp = k + 3 * ibit17
+          k_pp  = k + 2 * ( 1 - ibit15 )
+          k_mm  = k - 2 * ibit17
+
+          w_comp(k) = w(k,j,i) + w(k,j,i-1)
+          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
+                     ( 37.0_wp * ibit17 * adv_mom_5                           &
+                  +     7.0_wp * ibit16 * adv_mom_3                           &
+                  +              ibit15 * adv_mom_1                           &
+                     ) *                                                      &
+                                ( u(k+1,j,i)  + u(k,j,i)     )                &
+              -      (  8.0_wp * ibit17 * adv_mom_5                           &
+                  +              ibit16 * adv_mom_3                           &
+                     ) *                                                      &
+                                ( u(k_pp,j,i) + u(k-1,j,i)   )                &
+              +      (           ibit17 * adv_mom_5                           &
+                     ) *                                                      &
+                                ( u(k_ppp,j,i) + u(k_mm,j,i) )                &
+                                                  )
+
+          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
+                     ( 10.0_wp * ibit17 * adv_mom_5                           &
+                  +     3.0_wp * ibit16 * adv_mom_3                           &
+                  +              ibit15 * adv_mom_1                           &
+                     ) *                                                      &
+                                ( u(k+1,j,i)   - u(k,j,i)    )                &
+              -      (  5.0_wp * ibit17 * adv_mom_5                           &
+                  +              ibit16 * adv_mom_3                           &
+                     ) *                                                      &
+                                ( u(k_pp,j,i)  - u(k-1,j,i)  )                &
+              +      (           ibit17 * adv_mom_5                           &
+                     ) *                                                      &
+                                ( u(k_ppp,j,i) - u(k_mm,j,i) )                &
+                                                           )
+       ENDDO
+       
+       DO  k = nzb+1, nzt
 
           flux_d    = flux_t(k-1)
@@ -2057,30 +2201,32 @@
 !--       Calculate the divergence of the velocity field. A respective
 !--       correction is needed to overcome numerical instabilities introduced
-!--       by a not sufficient reduction of divergences near topography. 
+!--       by a not sufficient reduction of divergences near topography.
           div = ( ( u_comp(k)       * ( ibit9 + ibit10 + ibit11 )             &
                 - ( u(k,j,i)   + u(k,j,i-1)   )                               &
-                                    * ( IBITS(advc_flags_1(k,j,i-1),9,1)      &
-                                      + IBITS(advc_flags_1(k,j,i-1),10,1)     &
-                                      + IBITS(advc_flags_1(k,j,i-1),11,1)     &
+                                    * (                                       &
+                     REAL( IBITS(advc_flags_1(k,j,i-1),9,1),  KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )     &
                                       )                                       &
                   ) * ddx                                                     &
                +  ( ( v_comp(k) + gv ) * ( ibit12 + ibit13 + ibit14 )         &
                   - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
-                                    * ( IBITS(advc_flags_1(k,j-1,i),12,1)     &
-                                      + IBITS(advc_flags_1(k,j-1,i),13,1)     &
-                                      + IBITS(advc_flags_1(k,j-1,i),14,1)     &
+                                    * (                                       &
+                     REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )     &
                                       )                                       &
                   ) * ddy                                                     &
                +  ( w_comp(k)   * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )&
                 -   w_comp(k-1) * rho_air_zw(k-1)                             &
-                                    * ( IBITS(advc_flags_1(k-1,j,i),15,1)     &
-                                      + IBITS(advc_flags_1(k-1,j,i),16,1)     &
-                                      + IBITS(advc_flags_1(k-1,j,i),17,1)     &
+                                    * (                                       &
+                     REAL( IBITS(advc_flags_1(k-1,j,i),15,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k-1,j,i),16,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k-1,j,i),17,1), KIND = wp )     &
                                       )                                       &  
                   ) * drho_air(k) * ddzw(k)                                   &
                 ) * 0.5_wp
 
-
-          tend(k,j,i) = tend(k,j,i) - (                                       &
+           tend(k,j,i) = tend(k,j,i) - (                                      &
                             ( flux_r(k) + diss_r(k)                           &
                           -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx     &
@@ -2088,5 +2234,5 @@
                           -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
                           + ( ( flux_t(k) + diss_t(k) )                       &
-                          -   ( flux_d    + diss_d )                          &
+                          -   ( flux_d    + diss_d    )                       &
                                                     ) * drho_air(k) * ddzw(k) &
                                        ) + div * u(k,j,i)
@@ -2111,121 +2257,4 @@
            sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                         &
                 + ( flux_t(k)                                                  &
-                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                )     &
-                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )        )     &
-                  + diss_t(k)                                                  &
-                    *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)           )     &
-                    / ( ABS( w_comp(k) ) + 1.0E-20_wp                    )     &
-                  ) *   weight_substep(intermediate_timestep_count)
-       ENDDO
-
-       DO  k = nzb_max+1, nzt
-
-          u_comp(k) = u(k,j,i+1) + u(k,j,i)
-          flux_r(k) = ( u_comp(k) - gu ) * (                                  &
-                         37.0_wp * ( u(k,j,i+1) + u(k,j,i)   )                &
-                       -  8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) )                &
-                       +           ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
-          diss_r(k) = - ABS( u_comp(k) - gu ) * (                             &
-                         10.0_wp * ( u(k,j,i+1) - u(k,j,i)   )                &
-                       -  5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) )                &
-                       +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
-
-          v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
-          flux_n(k) = v_comp(k) * (                                           &
-                         37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                &
-                       -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                &
-                       +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
-          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
-                         10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                &
-                       -  5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) )                &
-                       +           ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
-!
-!--       k index has to be modified near bottom and top, else array
-!--       subscripts will be exceeded.
-          ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
-          ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
-          ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
-
-          k_ppp = k + 3 * ibit17
-          k_pp  = k + 2 * ( 1 - ibit15 )
-          k_mm  = k - 2 * ibit17
-
-          w_comp(k) = w(k,j,i) + w(k,j,i-1)
-          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
-                     ( 37.0_wp * ibit17 * adv_mom_5                           &
-                  +     7.0_wp * ibit16 * adv_mom_3                           &
-                  +              ibit15 * adv_mom_1                           &
-                     ) *                                                      &
-                                ( u(k+1,j,i)  + u(k,j,i)     )                &
-              -      (  8.0_wp * ibit17 * adv_mom_5                           &
-                  +              ibit16 * adv_mom_3                           &
-                     ) *                                                      &
-                                ( u(k_pp,j,i) + u(k-1,j,i)   )                &
-              +      (           ibit17 * adv_mom_5                           &
-                     ) *                                                      &
-                                ( u(k_ppp,j,i) + u(k_mm,j,i) )                &
-                                                  )
-
-          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
-                     ( 10.0_wp * ibit17 * adv_mom_5                           &
-                  +     3.0_wp * ibit16 * adv_mom_3                           &
-                  +              ibit15 * adv_mom_1                           &
-                     ) *                                                      &
-                                ( u(k+1,j,i)   - u(k,j,i)    )                &
-              -      (  5.0_wp * ibit17 * adv_mom_5                           &
-                  +              ibit16 * adv_mom_3                           &
-                     ) *                                                      &
-                                ( u(k_pp,j,i)  - u(k-1,j,i)  )                &
-              +      (           ibit17 * adv_mom_5                           &
-                     ) *                                                      &
-                                ( u(k_ppp,j,i) - u(k_mm,j,i) )                &
-                                                           )
-
-       ENDDO
-
-       DO  k = nzb_max+1, nzt
-
-          flux_d    = flux_t(k-1)
-          diss_d    = diss_t(k-1)
-!
-!--       Calculate the divergence of the velocity field. A respective
-!--       correction is needed to overcome numerical instabilities introduced
-!--       by a not sufficient reduction of divergences near topography.
-          div = ( ( u_comp(k)      - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx    &
-               +  ( v_comp(k) + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy    &
-               +  ( w_comp(k)   * rho_air_zw(k)   -                           &
-                    w_comp(k-1) * rho_air_zw(k-1)                             &
-                  ) * drho_air(k) * ddzw(k)                                   &
-                ) * 0.5_wp
-
-          tend(k,j,i) = tend(k,j,i) - (                                       &
-                            ( flux_r(k) + diss_r(k)                           &
-                          -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx     &
-                          + ( flux_n(k) + diss_n(k)                           &
-                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
-                          + ( ( flux_t(k) + diss_t(k) )                       &
-                          -   ( flux_d    + diss_d    )                       &
-                                                    ) * drho_air(k) * ddzw(k) &
-                                       ) + div * u(k,j,i)
-
-           flux_l_u(k,j,tn) = flux_r(k)
-           diss_l_u(k,j,tn) = diss_r(k)
-           flux_s_u(k,tn)   = flux_n(k)
-           diss_s_u(k,tn)   = diss_n(k)
-!
-!--        Statistical Evaluation of u'u'. The factor has to be applied for
-!--        right evaluation when gallilei_trans = .T. .
-           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                           &
-                + ( flux_r(k)                                                  &
-                    * ( u_comp(k) - 2.0_wp * hom(k,1,1,0)                   )  &
-                    / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )  &
-                  + diss_r(k)                                                  &
-                    *   ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0)              )  &
-                    / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp                  )  &
-                  ) *   weight_substep(intermediate_timestep_count)
-!
-!--        Statistical Evaluation of w'u'.
-           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                         &
-                + ( flux_t(k)                                                  &
                     * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
                     / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )  &
@@ -2272,13 +2301,4 @@
 
        INTEGER(iwp)  ::  i      !< grid index along x-direction
-       INTEGER(iwp)  ::  ibit18 !< flag indicating 1st-order scheme along x-direction
-       INTEGER(iwp)  ::  ibit19 !< flag indicating 3rd-order scheme along x-direction
-       INTEGER(iwp)  ::  ibit20 !< flag indicating 5th-order scheme along x-direction
-       INTEGER(iwp)  ::  ibit21 !< flag indicating 1st-order scheme along y-direction 
-       INTEGER(iwp)  ::  ibit22 !< flag indicating 3rd-order scheme along y-direction
-       INTEGER(iwp)  ::  ibit23 !< flag indicating 3rd-order scheme along y-direction
-       INTEGER(iwp)  ::  ibit24 !< flag indicating 1st-order scheme along z-direction
-       INTEGER(iwp)  ::  ibit25 !< flag indicating 3rd-order scheme along z-direction
-       INTEGER(iwp)  ::  ibit26 !< flag indicating 3rd-order scheme along z-direction
        INTEGER(iwp)  ::  i_omp  !< leftmost index on subdomain, or in case of OpenMP, on thread
        INTEGER(iwp)  ::  j      !< grid index along y-direction
@@ -2289,10 +2309,19 @@
        INTEGER(iwp)  ::  tn     !< number of OpenMP thread
        
-       REAL(wp)     ::  diss_d   !< artificial dissipation term at grid box bottom
-       REAL(wp)     ::  div      !< divergence on v-grid
-       REAL(wp)     ::  flux_d   !< 6th-order flux at grid box bottom
-       REAL(wp)     ::  gu       !< Galilei-transformation velocity along x
-       REAL(wp)     ::  gv       !< Galilei-transformation velocity along y
-       REAL(wp)     ::  v_comp_l !< advection velocity along y on leftmost grid point on subdomain
+       REAL(wp)      ::  ibit18   !< flag indicating 1st-order scheme along x-direction
+       REAL(wp)      ::  ibit19   !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp)      ::  ibit20   !< flag indicating 5th-order scheme along x-direction
+       REAL(wp)      ::  ibit21   !< flag indicating 1st-order scheme along y-direction 
+       REAL(wp)      ::  ibit22   !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp)      ::  ibit23   !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp)      ::  ibit24   !< flag indicating 1st-order scheme along z-direction
+       REAL(wp)      ::  ibit25   !< flag indicating 3rd-order scheme along z-direction
+       REAL(wp)      ::  ibit26   !< flag indicating 3rd-order scheme along z-direction
+       REAL(wp)      ::  diss_d   !< artificial dissipation term at grid box bottom
+       REAL(wp)      ::  div      !< divergence on v-grid
+       REAL(wp)      ::  flux_d   !< 6th-order flux at grid box bottom
+       REAL(wp)      ::  gu       !< Galilei-transformation velocity along x
+       REAL(wp)      ::  gv       !< Galilei-transformation velocity along y
+       REAL(wp)      ::  v_comp_l !< advection velocity along y on leftmost grid point on subdomain
        
        REAL(wp), DIMENSION(nzb:nzt+1)  ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
@@ -2315,7 +2344,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit20 = IBITS(advc_flags_1(k,j,i-1),20,1)
-             ibit19 = IBITS(advc_flags_1(k,j,i-1),19,1)
-             ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)
+             ibit20 = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )
+             ibit19 = REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp )
+             ibit18 = REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp )
 
              u_comp(k)        = u(k,j-1,i) + u(k,j,i) - gu
@@ -2373,7 +2402,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit23 = IBITS(advc_flags_1(k,j-1,i),23,1)
-             ibit22 = IBITS(advc_flags_1(k,j-1,i),22,1)
-             ibit21 = IBITS(advc_flags_1(k,j-1,i),21,1)
+             ibit23 = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )
+             ibit22 = REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp )
+             ibit21 = REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )
 
              v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
@@ -2425,8 +2454,4 @@
           
        ENDIF
-
-       flux_t(0) = 0.0_wp
-       diss_t(0) = 0.0_wp
-       w_comp(0) = 0.0_wp
 !
 !--    Now compute the fluxes and tendency terms for the horizontal and
@@ -2434,7 +2459,7 @@
        DO  k = nzb+1, nzb_max
 
-          ibit20 = IBITS(advc_flags_1(k,j,i),20,1)
-          ibit19 = IBITS(advc_flags_1(k,j,i),19,1)
-          ibit18 = IBITS(advc_flags_1(k,j,i),18,1)
+          ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp )
+          ibit19 = REAL( IBITS(advc_flags_1(k,j,i),19,1), KIND = wp )
+          ibit18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp )
  
           u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
@@ -2469,7 +2494,7 @@
                                            )
 
-          ibit23 = IBITS(advc_flags_1(k,j,i),23,1)
-          ibit22 = IBITS(advc_flags_1(k,j,i),22,1)
-          ibit21 = IBITS(advc_flags_1(k,j,i),21,1)
+          ibit23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp )
+          ibit22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp )
+          ibit21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp )
 
 
@@ -2504,10 +2529,51 @@
                                     ( v(k,j+3,i) - v(k,j-2,i) )               &
                                                 )
+       ENDDO
+
+       DO  k = nzb_max+1, nzt
+
+          u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
+          flux_r(k) = u_comp(k) * (                                           &
+                      37.0_wp * ( v(k,j,i+1) + v(k,j,i)   )                   &
+                    -  8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) )                   &
+                    +           ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
+
+          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
+                      10.0_wp * ( v(k,j,i+1) - v(k,j,i) )                     &
+                    -  5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) )                   &
+                    +           ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
+
+
+          v_comp(k) = v(k,j+1,i) + v(k,j,i)
+          flux_n(k) = ( v_comp(k) - gv ) * (                                  &
+                      37.0_wp * ( v(k,j+1,i) + v(k,j,i)   )                   &
+                    -  8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) )                   &
+                      +         ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
+
+          diss_n(k) = - ABS( v_comp(k) - gv ) * (                             &
+                      10.0_wp * ( v(k,j+1,i) - v(k,j,i)   )                   &
+                    -  5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) )                   &
+                    +           ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
+       ENDDO
+!
+!--    Now, compute vertical fluxes. Split loop into a part treating the 
+!--    lowest 2 grid points with indirect indexing, a main loop without
+!--    indirect indexing, and a loop for the uppermost 2 grip points with
+!--    indirect indexing. This allows better vectorization for the main loop.
+!--    First, compute the flux at model surface, which need has to be 
+!--    calculated explicetely for the tendency at
+!--    the first w-level. For topography wall this is done implicitely by
+!--    advc_flags_1.
+       flux_t(nzb) = 0.0_wp
+       diss_t(nzb) = 0.0_wp
+       w_comp(nzb) = 0.0_wp
+       
+       DO  k = nzb+1, nzb+2
 !
 !--       k index has to be modified near bottom and top, else array
 !--       subscripts will be exceeded.
-          ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
-          ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
-          ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
+          ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
+          ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
+          ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
 
           k_ppp = k + 3 * ibit26
@@ -2546,6 +2612,88 @@
                                                            )
        ENDDO
-
-       DO  k = nzb+1, nzb_max
+       
+       DO  k = nzb+3, nzt-2
+
+          ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
+          ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
+          ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
+
+          w_comp(k) = w(k,j-1,i) + w(k,j,i)
+          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
+                     ( 37.0_wp * ibit26 * adv_mom_5                           &
+                  +     7.0_wp * ibit25 * adv_mom_3                           &
+                  +              ibit24 * adv_mom_1                           &
+                     ) *                                                      &
+                                ( v(k+1,j,i)   + v(k,j,i)    )                &
+              -      (  8.0_wp * ibit26 * adv_mom_5                           &
+                  +              ibit25 * adv_mom_3                           &
+                     ) *                                                      &
+                                ( v(k+2,j,i)  + v(k-1,j,i)  )                 &
+              +      (           ibit26 * adv_mom_5                           &
+                     ) *                                                      &
+                                ( v(k+3,j,i) + v(k-2,j,i) )                   &
+                                                  )
+
+          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
+                     ( 10.0_wp * ibit26 * adv_mom_5                           &
+                  +     3.0_wp * ibit25 * adv_mom_3                           &
+                  +              ibit24 * adv_mom_1                           &
+                     ) *                                                      &
+                                ( v(k+1,j,i)   - v(k,j,i)    )                &
+              -      (  5.0_wp * ibit26 * adv_mom_5                           &
+                  +              ibit25 * adv_mom_3                           &
+                     ) *                                                      &
+                                ( v(k+2,j,i)  - v(k-1,j,i)  )                 &
+              +      (           ibit26 * adv_mom_5                           &
+                     ) *                                                      &
+                                ( v(k+3,j,i) - v(k-2,j,i) )                   &
+                                                           )
+       ENDDO
+       
+       DO  k = nzt-1, nzt
+!
+!--       k index has to be modified near bottom and top, else array
+!--       subscripts will be exceeded.
+          ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
+          ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
+          ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
+
+          k_ppp = k + 3 * ibit26
+          k_pp  = k + 2 * ( 1 - ibit24  )
+          k_mm  = k - 2 * ibit26
+
+          w_comp(k) = w(k,j-1,i) + w(k,j,i)
+          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
+                     ( 37.0_wp * ibit26 * adv_mom_5                           &
+                  +     7.0_wp * ibit25 * adv_mom_3                           &
+                  +              ibit24 * adv_mom_1                           &
+                     ) *                                                      &
+                                ( v(k+1,j,i)   + v(k,j,i)    )                &
+              -      (  8.0_wp * ibit26 * adv_mom_5                           &
+                  +              ibit25 * adv_mom_3                           &
+                     ) *                                                      &
+                                ( v(k_pp,j,i)  + v(k-1,j,i)  )                &
+              +      (           ibit26 * adv_mom_5                           &
+                     ) *                                                      &
+                                ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
+                                                  )
+
+          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
+                     ( 10.0_wp * ibit26 * adv_mom_5                           &
+                  +     3.0_wp * ibit25 * adv_mom_3                           &
+                  +              ibit24 * adv_mom_1                           &
+                     ) *                                                      &
+                                ( v(k+1,j,i)   - v(k,j,i)    )                &
+              -      (  5.0_wp * ibit26 * adv_mom_5                           &
+                  +              ibit25 * adv_mom_3                           &
+                     ) *                                                      &
+                                ( v(k_pp,j,i)  - v(k-1,j,i)  )                &
+              +      (           ibit26 * adv_mom_5                           &
+                     ) *                                                      &
+                                ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
+                                                           )
+       ENDDO
+       
+       DO  k = nzb+1, nzt
 
           flux_d    = flux_t(k-1)
@@ -2554,11 +2702,12 @@
 !--       Calculate the divergence of the velocity field. A respective
 !--       correction is needed to overcome numerical instabilities introduced
-!--       by a not sufficient reduction of divergences near topography. 
+!--       by a not sufficient reduction of divergences near topography.
           div = ( ( ( u_comp(k)     + gu )                                    &
                                        * ( ibit18 + ibit19 + ibit20 )         &
                   - ( u(k,j-1,i) + u(k,j,i) )                                 &
-                                       * ( IBITS(advc_flags_1(k,j,i-1),18,1)  &
-                                         + IBITS(advc_flags_1(k,j,i-1),19,1)  &
-                                         + IBITS(advc_flags_1(k,j,i-1),20,1)  &
+                                       * (                                    &
+                        REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp )  &
+                      + REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp )  &
+                      + REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )  &
                                          )                                    &
                   ) * ddx                                                     &
@@ -2566,18 +2715,19 @@
                                        * ( ibit21 + ibit22 + ibit23 )         &
                 - ( v(k,j,i)     + v(k,j-1,i) )                               &
-                                       * ( IBITS(advc_flags_1(k,j-1,i),21,1)  &
-                                         + IBITS(advc_flags_1(k,j-1,i),22,1)  &
-                                         + IBITS(advc_flags_1(k,j-1,i),23,1)  &
+                                       * (                                    &
+                        REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )  &
+                      + REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp )  &
+                      + REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )  &
                                          )                                    &
                   ) * ddy                                                     &
                +  ( w_comp(k)   * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )&
                 -   w_comp(k-1) * rho_air_zw(k-1)                             &
-                                       * ( IBITS(advc_flags_1(k-1,j,i),24,1)  &
-                                         + IBITS(advc_flags_1(k-1,j,i),25,1)  &
-                                         + IBITS(advc_flags_1(k-1,j,i),26,1)  &
+                                       * (                                    &
+                        REAL( IBITS(advc_flags_1(k-1,j,i),24,1), KIND = wp )  &
+                      + REAL( IBITS(advc_flags_1(k-1,j,i),25,1), KIND = wp )  &
+                      + REAL( IBITS(advc_flags_1(k-1,j,i),26,1), KIND = wp )  &
                                          )                                    &
                    ) * drho_air(k) * ddzw(k)                                  &
                 ) * 0.5_wp
-
 
           tend(k,j,i) = tend(k,j,i) - (                                       &
@@ -2595,5 +2745,4 @@
            flux_s_v(k,tn)   = flux_n(k)
            diss_s_v(k,tn)   = diss_n(k)
-
 !
 !--        Statistical Evaluation of v'v'. The factor has to be applied for
@@ -2619,124 +2768,4 @@
 
        ENDDO
-
-       DO  k = nzb_max+1, nzt
-
-          u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
-          flux_r(k) = u_comp(k) * (                                           &
-                      37.0_wp * ( v(k,j,i+1) + v(k,j,i)   )                   &
-                    -  8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) )                   &
-                    +           ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
-
-          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
-                      10.0_wp * ( v(k,j,i+1) - v(k,j,i) )                     &
-                    -  5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) )                   &
-                    +           ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
-
-
-          v_comp(k) = v(k,j+1,i) + v(k,j,i)
-          flux_n(k) = ( v_comp(k) - gv ) * (                                  &
-                      37.0_wp * ( v(k,j+1,i) + v(k,j,i)   )                   &
-                    -  8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) )                   &
-                      +         ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
-
-          diss_n(k) = - ABS( v_comp(k) - gv ) * (                             &
-                      10.0_wp * ( v(k,j+1,i) - v(k,j,i)   )                   &
-                    -  5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) )                   &
-                    +           ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
-!
-!--       k index has to be modified near bottom and top, else array
-!--       subscripts will be exceeded.
-          ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
-          ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
-          ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
-
-          k_ppp = k + 3 * ibit26
-          k_pp  = k + 2 * ( 1 - ibit24  )
-          k_mm  = k - 2 * ibit26
-
-          w_comp(k) = w(k,j-1,i) + w(k,j,i)
-          flux_t(k) = w_comp(k) * rho_air_zw(k) * (                           &
-                     ( 37.0_wp * ibit26 * adv_mom_5                           &
-                  +     7.0_wp * ibit25 * adv_mom_3                           &
-                  +              ibit24 * adv_mom_1                           &
-                     ) *                                                      &
-                                ( v(k+1,j,i)   + v(k,j,i)    )                &
-              -      (  8.0_wp * ibit26 * adv_mom_5                           &
-                  +              ibit25 * adv_mom_3                           &
-                     ) *                                                      &
-                                ( v(k_pp,j,i)  + v(k-1,j,i)  )                &
-              +      (           ibit26 * adv_mom_5                           &
-                     ) *                                                      &
-                                ( v(k_ppp,j,i) + v(k_mm,j,i) )                &
-                                                  )
-
-          diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * (                  &
-                     ( 10.0_wp * ibit26 * adv_mom_5                           &
-                  +     3.0_wp * ibit25 * adv_mom_3                           &
-                  +              ibit24 * adv_mom_1                           &
-                     ) *                                                      &
-                                ( v(k+1,j,i)   - v(k,j,i)    )                &
-              -      (  5.0_wp * ibit26 * adv_mom_5                           &
-                  +              ibit25 * adv_mom_3                           &
-                     ) *                                                      &
-                                ( v(k_pp,j,i)  - v(k-1,j,i)  )                &
-              +      (           ibit26 * adv_mom_5                           &
-                     ) *                                                      &
-                                ( v(k_ppp,j,i) - v(k_mm,j,i) )                &
-                                                           )
-       ENDDO
-
-       DO  k = nzb_max+1, nzt
-
-          flux_d    = flux_t(k-1)
-          diss_d    = diss_t(k-1)
-!
-!--       Calculate the divergence of the velocity field. A respective
-!--       correction is needed to overcome numerical instabilities introduced
-!--       by a not sufficient reduction of divergences near topography.
-          div = ( ( u_comp(k) + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx    &
-               +  ( v_comp(k)      - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy    &
-               +  ( w_comp(k)   * rho_air_zw(k)   -                           &
-                    w_comp(k-1) * rho_air_zw(k-1)                             &
-                  ) * drho_air(k) * ddzw(k)                                   &
-                ) * 0.5_wp
-
-          tend(k,j,i) = tend(k,j,i) - (                                       &
-                         ( flux_r(k) + diss_r(k)                              &
-                       -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx      &
-                       + ( flux_n(k) + diss_n(k)                              &
-                       -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy      &
-                       + ( ( flux_t(k) + diss_t(k) )                          &
-                       -   ( flux_d    + diss_d    )                          &
-                                                   ) * drho_air(k) * ddzw(k)  &
-                                      ) + v(k,j,i) * div
-
-           flux_l_v(k,j,tn) = flux_r(k)
-           diss_l_v(k,j,tn) = diss_r(k)
-           flux_s_v(k,tn)   = flux_n(k)
-           diss_s_v(k,tn)   = diss_n(k)
-!
-!--        Statistical Evaluation of v'v'. The factor has to be applied for
-!--        right evaluation when gallilei_trans = .T. .
-           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                           &
-                + ( flux_n(k)                                                  &
-                    * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)                   )  &
-                    / ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) )  &
-                  + diss_n(k)                                                  &
-                    *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0)              )  &
-                    / ( ABS( v_comp(k) - gv ) + 1.0E-20_wp                  )  &
-                  ) *   weight_substep(intermediate_timestep_count)
-!
-!--        Statistical Evaluation of w'u'.
-           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                         &
-                + ( flux_t(k)                                                  &
-                    * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                   )  &
-                    / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )           )  &
-                  + diss_t(k)                                                  &
-                    *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)              )  &
-                    / ( ABS( w_comp(k) ) + 1.0E-20_wp                       )  &
-                  ) *   weight_substep(intermediate_timestep_count)
-
-       ENDDO
        sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
 
@@ -2774,13 +2803,4 @@
 
        INTEGER(iwp) ::  i      !< grid index along x-direction
-       INTEGER(iwp) ::  ibit27 !< flag indicating 1st-order scheme along x-direction
-       INTEGER(iwp) ::  ibit28 !< flag indicating 3rd-order scheme along x-direction
-       INTEGER(iwp) ::  ibit29 !< flag indicating 5th-order scheme along x-direction
-       INTEGER(iwp) ::  ibit30 !< flag indicating 1st-order scheme along y-direction
-       INTEGER(iwp) ::  ibit31 !< flag indicating 3rd-order scheme along y-direction
-       INTEGER(iwp) ::  ibit32 !< flag indicating 5th-order scheme along y-direction
-       INTEGER(iwp) ::  ibit33 !< flag indicating 1st-order scheme along z-direction
-       INTEGER(iwp) ::  ibit34 !< flag indicating 3rd-order scheme along z-direction
-       INTEGER(iwp) ::  ibit35 !< flag indicating 5th-order scheme along z-direction
        INTEGER(iwp) ::  i_omp  !< leftmost index on subdomain, or in case of OpenMP, on thread
        INTEGER(iwp) ::  j      !< grid index along y-direction
@@ -2791,4 +2811,13 @@
        INTEGER(iwp) ::  tn     !< number of OpenMP thread
        
+       REAL(wp)    ::  ibit27  !< flag indicating 1st-order scheme along x-direction
+       REAL(wp)    ::  ibit28  !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp)    ::  ibit29  !< flag indicating 5th-order scheme along x-direction
+       REAL(wp)    ::  ibit30  !< flag indicating 1st-order scheme along y-direction
+       REAL(wp)    ::  ibit31  !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp)    ::  ibit32  !< flag indicating 5th-order scheme along y-direction
+       REAL(wp)    ::  ibit33  !< flag indicating 1st-order scheme along z-direction
+       REAL(wp)    ::  ibit34  !< flag indicating 3rd-order scheme along z-direction
+       REAL(wp)    ::  ibit35  !< flag indicating 5th-order scheme along z-direction
        REAL(wp)    ::  diss_d  !< discretized artificial dissipation at top of the grid box
        REAL(wp)    ::  div     !< divergence on w-grid
@@ -2815,7 +2844,7 @@
 
           DO  k = nzb+1, nzb_max
-             ibit32 = IBITS(advc_flags_2(k,j-1,i),0,1)
-             ibit31 = IBITS(advc_flags_1(k,j-1,i),31,1)
-             ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)
+             ibit32 = REAL( IBITS(advc_flags_2(k,j-1,i),0,1),  KIND = wp )
+             ibit31 = REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )
+             ibit30 = REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp )
 
              v_comp(k)      = v(k+1,j,i) + v(k,j,i) - gv
@@ -2859,5 +2888,5 @@
                          -  8.0_wp * ( w(k,j+1,i) +w(k,j-2,i)  )              &
                          +           ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
-             diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                         &
+             diss_s_w(k,tn) = - ABS( v_comp(k) ) * (                          &
                            10.0_wp * ( w(k,j,i) - w(k,j-1,i)   )              &
                          -  5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) )              &
@@ -2873,7 +2902,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit29 = IBITS(advc_flags_1(k,j,i-1),29,1)
-             ibit28 = IBITS(advc_flags_1(k,j,i-1),28,1)
-             ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)
+             ibit29 = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )
+             ibit28 = REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp )
+             ibit27 = REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp )
 
              u_comp(k)        = u(k+1,j,i) + u(k,j,i) - gu
@@ -2926,5 +2955,111 @@
        ENDIF
 !
-!--    The lower flux has to be calculated explicetely for the tendency at
+!--    Now compute the fluxes and tendency terms for the horizontal
+!--    and vertical parts.
+       DO  k = nzb+1, nzb_max
+
+          ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp )
+          ibit28 = REAL( IBITS(advc_flags_1(k,j,i),28,1), KIND = wp )
+          ibit27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp )
+
+          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
+          flux_r(k) = u_comp(k) * (                                           &
+                     ( 37.0_wp * ibit29 * adv_mom_5                           &
+                  +     7.0_wp * ibit28 * adv_mom_3                           &
+                  +              ibit27 * adv_mom_1                           &
+                     ) *                                                      &
+                                    ( w(k,j,i+1) + w(k,j,i)   )               &
+              -      (  8.0_wp * ibit29 * adv_mom_5                           &
+                  +              ibit28 * adv_mom_3                           &
+                     ) *                                                      &
+                                    ( w(k,j,i+2) + w(k,j,i-1) )               &
+              +      (           ibit29 * adv_mom_5                           &
+                     ) *                                                      &
+                                    ( w(k,j,i+3) + w(k,j,i-2) )               &
+                                  )
+
+          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
+                     ( 10.0_wp * ibit29 * adv_mom_5                           &
+                  +     3.0_wp * ibit28 * adv_mom_3                           &
+                  +              ibit27 * adv_mom_1                           &
+                     ) *                                                      &
+                                    ( w(k,j,i+1) - w(k,j,i)   )               &
+              -      (  5.0_wp * ibit29 * adv_mom_5                           &
+                  +              ibit28 * adv_mom_3                           &
+                     ) *                                                      &
+                                    ( w(k,j,i+2) - w(k,j,i-1) )               &
+              +      (           ibit29 * adv_mom_5                           &
+                     ) *                                                      &
+                                    ( w(k,j,i+3) - w(k,j,i-2) )               &
+                                           )
+
+          ibit32 = REAL( IBITS(advc_flags_2(k,j,i),0,1),  KIND = wp )
+          ibit31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp )
+          ibit30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp )
+
+          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
+          flux_n(k) = v_comp(k) * (                                           &
+                     ( 37.0_wp * ibit32 * adv_mom_5                           &
+                  +     7.0_wp * ibit31 * adv_mom_3                           &
+                  +              ibit30 * adv_mom_1                           &
+                     ) *                                                      &
+                                    ( w(k,j+1,i) + w(k,j,i)   )               &
+              -      (  8.0_wp * ibit32 * adv_mom_5                           &
+                  +              ibit31 * adv_mom_3                           &
+                     ) *                                                      &
+                                    ( w(k,j+2,i) + w(k,j-1,i) )               &
+              +      (           ibit32 * adv_mom_5                           &
+                     ) *                                                      &
+                                    ( w(k,j+3,i) + w(k,j-2,i) )               &
+                                  )
+
+          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
+                     ( 10.0_wp * ibit32 * adv_mom_5                           &
+                  +     3.0_wp * ibit31 * adv_mom_3                           &
+                  +              ibit30 * adv_mom_1                           &
+                     ) *                                                      &
+                                    ( w(k,j+1,i) - w(k,j,i)  )                &
+              -      (  5.0_wp * ibit32 * adv_mom_5                           &
+                  +              ibit31 * adv_mom_3                           &
+                     ) *                                                      &
+                                   ( w(k,j+2,i) - w(k,j-1,i) )                &
+              +      (           ibit32 * adv_mom_5                           &
+                     ) *                                                      &
+                                   ( w(k,j+3,i) - w(k,j-2,i) )                &
+                                           )
+       ENDDO
+
+       DO  k = nzb_max+1, nzt
+
+          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
+          flux_r(k) = u_comp(k) * (                                           &
+                      37.0_wp * ( w(k,j,i+1) + w(k,j,i)   )                   &
+                    -  8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) )                   &
+                    +           ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
+
+          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
+                      10.0_wp * ( w(k,j,i+1) - w(k,j,i)   )                   &
+                    -  5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) )                   &
+                    +           ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
+
+          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
+          flux_n(k) = v_comp(k) * (                                           &
+                      37.0_wp * ( w(k,j+1,i) + w(k,j,i)   )                   &
+                    -  8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) )                   &
+                    +           ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
+
+          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
+                      10.0_wp * ( w(k,j+1,i) - w(k,j,i)   )                   &
+                    -  5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) )                   &
+                    +           ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
+       ENDDO
+
+!
+!--    Now, compute vertical fluxes. Split loop into a part treating the 
+!--    lowest 2 grid points with indirect indexing, a main loop without
+!--    indirect indexing, and a loop for the uppermost 2 grip points with
+!--    indirect indexing. This allows better vectorization for the main loop.
+!--    First, compute the flux at model surface, which need has to be 
+!--    calculated explicetely for the tendency at
 !--    the first w-level. For topography wall this is done implicitely by
 !--    advc_flags_1.
@@ -2933,84 +3068,12 @@
        flux_t(0) = w_comp(k)       * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
        diss_t(0) = -ABS(w_comp(k)) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
-!
-!--    Now compute the fluxes and tendency terms for the horizontal
-!--    and vertical parts.
-       DO  k = nzb+1, nzb_max
-
-          ibit29 = IBITS(advc_flags_1(k,j,i),29,1)
-          ibit28 = IBITS(advc_flags_1(k,j,i),28,1)
-          ibit27 = IBITS(advc_flags_1(k,j,i),27,1)
-
-          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
-          flux_r(k) = u_comp(k) * (                                           &
-                     ( 37.0_wp * ibit29 * adv_mom_5                           &
-                  +     7.0_wp * ibit28 * adv_mom_3                           &
-                  +              ibit27 * adv_mom_1                           &
-                     ) *                                                      &
-                                    ( w(k,j,i+1) + w(k,j,i)   )               &
-              -      (  8.0_wp * ibit29 * adv_mom_5                           &
-                  +              ibit28 * adv_mom_3                           &
-                     ) *                                                      &
-                                    ( w(k,j,i+2) + w(k,j,i-1) )               &
-              +      (           ibit29 * adv_mom_5                           &
-                     ) *                                                      &
-                                    ( w(k,j,i+3) + w(k,j,i-2) )               &
-                                  )
-
-          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
-                     ( 10.0_wp * ibit29 * adv_mom_5                           &
-                  +     3.0_wp * ibit28 * adv_mom_3                           &
-                  +              ibit27 * adv_mom_1                           &
-                     ) *                                                      &
-                                    ( w(k,j,i+1) - w(k,j,i)   )               &
-              -      (  5.0_wp * ibit29 * adv_mom_5                           &
-                  +              ibit28 * adv_mom_3                           &
-                     ) *                                                      &
-                                    ( w(k,j,i+2) - w(k,j,i-1) )               &
-              +      (           ibit29 * adv_mom_5                           &
-                     ) *                                                      &
-                                    ( w(k,j,i+3) - w(k,j,i-2) )               &
-                                           )
-
-          ibit32 = IBITS(advc_flags_2(k,j,i),0,1)
-          ibit31 = IBITS(advc_flags_1(k,j,i),31,1)
-          ibit30 = IBITS(advc_flags_1(k,j,i),30,1)
-
-          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
-          flux_n(k) = v_comp(k) * (                                           &
-                     ( 37.0_wp * ibit32 * adv_mom_5                           &
-                  +     7.0_wp * ibit31 * adv_mom_3                           &
-                  +              ibit30 * adv_mom_1                           &
-                     ) *                                                      &
-                                    ( w(k,j+1,i) + w(k,j,i)   )               &
-              -      (  8.0_wp * ibit32 * adv_mom_5                           &
-                  +              ibit31 * adv_mom_3                           &
-                     ) *                                                      &
-                                    ( w(k,j+2,i) + w(k,j-1,i) )               &
-              +      (           ibit32 * adv_mom_5                           &
-                     ) *                                                      &
-                                    ( w(k,j+3,i) + w(k,j-2,i) )               &
-                                  )
-
-          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
-                     ( 10.0_wp * ibit32 * adv_mom_5                           &
-                  +     3.0_wp * ibit31 * adv_mom_3                           &
-                  +              ibit30 * adv_mom_1                           &
-                     ) *                                                      &
-                                    ( w(k,j+1,i) - w(k,j,i)  )                &
-              -      (  5.0_wp * ibit32 * adv_mom_5                           &
-                  +              ibit31 * adv_mom_3                           &
-                     ) *                                                      &
-                                   ( w(k,j+2,i) - w(k,j-1,i) )                &
-              +      (           ibit32 * adv_mom_5                           &
-                     ) *                                                      &
-                                   ( w(k,j+3,i) - w(k,j-2,i) )                &
-                                           )
+       
+       DO  k = nzb+1, nzb+2
 !
 !--       k index has to be modified near bottom and top, else array
 !--       subscripts will be exceeded.
-          ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
-          ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
-          ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
+          ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
+          ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
+          ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
 
           k_ppp = k + 3 * ibit35
@@ -3049,6 +3112,88 @@
                                                           )
        ENDDO
-
-       DO  k = nzb+1, nzb_max
+       
+       DO  k = nzb+3, nzt-2
+       
+          ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
+          ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
+          ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
+
+          w_comp(k) = w(k+1,j,i) + w(k,j,i)
+          flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
+                     ( 37.0_wp * ibit35 * adv_mom_5                           &
+                  +     7.0_wp * ibit34 * adv_mom_3                           &
+                  +              ibit33 * adv_mom_1                           &
+                     ) *                                                      &
+                                ( w(k+1,j,i)  + w(k,j,i)     )                &
+              -      (  8.0_wp * ibit35 * adv_mom_5                           &
+                  +              ibit34 * adv_mom_3                           &
+                     ) *                                                      &
+                                ( w(k+2,j,i)  + w(k-1,j,i)  )                 &
+              +      (           ibit35 * adv_mom_5                           &
+                     ) *                                                      &
+                                ( w(k+3,j,i) + w(k-2,j,i) )                   &
+                                                 )
+
+          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
+                     ( 10.0_wp * ibit35 * adv_mom_5                           &
+                  +     3.0_wp * ibit34 * adv_mom_3                           &
+                  +              ibit33 * adv_mom_1                           &
+                     ) *                                                      &
+                                ( w(k+1,j,i)   - w(k,j,i)    )                &
+              -      (  5.0_wp * ibit35 * adv_mom_5                           &
+                  +              ibit34 * adv_mom_3                           &
+                     ) *                                                      &
+                                ( w(k+2,j,i)  - w(k-1,j,i)  )                 &
+              +      (           ibit35 * adv_mom_5                           &
+                     ) *                                                      &
+                                ( w(k+3,j,i) - w(k-2,j,i) )                   &
+                                                          )
+       ENDDO
+       
+       DO  k = nzt-1, nzt
+!
+!--       k index has to be modified near bottom and top, else array
+!--       subscripts will be exceeded.
+          ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
+          ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
+          ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
+
+          k_ppp = k + 3 * ibit35
+          k_pp  = k + 2 * ( 1 - ibit33  )
+          k_mm  = k - 2 * ibit35
+
+          w_comp(k) = w(k+1,j,i) + w(k,j,i)
+          flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
+                     ( 37.0_wp * ibit35 * adv_mom_5                           &
+                  +     7.0_wp * ibit34 * adv_mom_3                           &
+                  +              ibit33 * adv_mom_1                           &
+                     ) *                                                      &
+                                ( w(k+1,j,i)  + w(k,j,i)     )                &
+              -      (  8.0_wp * ibit35 * adv_mom_5                           &
+                  +              ibit34 * adv_mom_3                           &
+                     ) *                                                      &
+                                ( w(k_pp,j,i)  + w(k-1,j,i)  )                &
+              +      (           ibit35 * adv_mom_5                           &
+                     ) *                                                      &
+                                ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
+                                                 )
+
+          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
+                     ( 10.0_wp * ibit35 * adv_mom_5                           &
+                  +     3.0_wp * ibit34 * adv_mom_3                           &
+                  +              ibit33 * adv_mom_1                           &
+                     ) *                                                      &
+                                ( w(k+1,j,i)   - w(k,j,i)    )                &
+              -      (  5.0_wp * ibit35 * adv_mom_5                           &
+                  +              ibit34 * adv_mom_3                           &
+                     ) *                                                      &
+                                ( w(k_pp,j,i)  - w(k-1,j,i)  )                &
+              +      (           ibit35 * adv_mom_5                           &
+                     ) *                                                      &
+                                ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
+                                                          )
+       ENDDO
+       
+       DO  k = nzb+1, nzt
 
           flux_d    = flux_t(k-1)
@@ -3057,17 +3202,19 @@
 !--       Calculate the divergence of the velocity field. A respective
 !--       correction is needed to overcome numerical instabilities introduced
-!--       by a not sufficient reduction of divergences near topography. 
+!--       by a not sufficient reduction of divergences near topography.
           div = ( ( ( u_comp(k) + gu ) * ( ibit27 + ibit28 + ibit29 )         &
                   - ( u(k+1,j,i) + u(k,j,i)   )                               & 
-                                    * ( IBITS(advc_flags_1(k,j,i-1),27,1)     &
-                                      + IBITS(advc_flags_1(k,j,i-1),28,1)     &
-                                      + IBITS(advc_flags_1(k,j,i-1),29,1)     &
+                                    * (                                       &
+                     REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )     &
                                       )                                       &
                   ) * ddx                                                     &
               +   ( ( v_comp(k) + gv ) * ( ibit30 + ibit31 + ibit32 )         & 
                   - ( v(k+1,j,i) + v(k,j,i)   )                               &
-                                    * ( IBITS(advc_flags_1(k,j-1,i),30,1)     &
-                                      + IBITS(advc_flags_1(k,j-1,i),31,1)     &
-                                      + IBITS(advc_flags_2(k,j-1,i),0,1)      &
+                                    * (                                       &
+                     REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_2(k,j-1,i),0,1),  KIND = wp )     &
                                       )                                       &
                   ) * ddy                                                     &
@@ -3075,117 +3222,9 @@
                                             * ( ibit33 + ibit34 + ibit35 )    &
                 - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k)                      &
-                                    * ( IBITS(advc_flags_2(k-1,j,i),1,1)      &
-                                      + IBITS(advc_flags_2(k-1,j,i),2,1)      &
-                                      + IBITS(advc_flags_2(k-1,j,i),3,1)      &
+                                    * (                                       &
+                     REAL( IBITS(advc_flags_2(k-1,j,i),1,1), KIND = wp )      &
+                   + REAL( IBITS(advc_flags_2(k-1,j,i),2,1), KIND = wp )      &
+                   + REAL( IBITS(advc_flags_2(k-1,j,i),3,1), KIND = wp )      &
                                       )                                       & 
-                  ) * drho_air_zw(k) * ddzu(k+1)                              &
-                ) * 0.5_wp
-
-
-          tend(k,j,i) = tend(k,j,i) - (                                       &
-                      ( flux_r(k) + diss_r(k)                                 &
-                    -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx         &
-                    + ( flux_n(k) + diss_n(k)                                 &
-                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy         &
-                    + ( ( flux_t(k) + diss_t(k) )                             &
-                    -   ( flux_d    + diss_d    )                             &
-                                              ) * drho_air_zw(k) * ddzu(k+1)  &
-                                      ) + div * w(k,j,i)
-
-          flux_l_w(k,j,tn) = flux_r(k)
-          diss_l_w(k,j,tn) = diss_r(k)
-          flux_s_w(k,tn)   = flux_n(k)
-          diss_s_w(k,tn)   = diss_n(k)
-!
-!--       Statistical Evaluation of w'w'.
-          sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                          &
-                      + ( flux_t(k)                                           &
-                       * ( w_comp(k) - 2.0_wp * hom(k,1,3,0)                ) &
-                       / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) )        ) &
-                        + diss_t(k)                                           &
-                       *   ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0)           ) &
-                       / ( ABS( w_comp(k) ) + 1.0E-20_wp                    ) &
-                        ) *   weight_substep(intermediate_timestep_count)
-
-       ENDDO
-
-       DO  k = nzb_max+1, nzt
-
-          u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
-          flux_r(k) = u_comp(k) * (                                           &
-                      37.0_wp * ( w(k,j,i+1) + w(k,j,i)   )                   &
-                    -  8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) )                   &
-                    +           ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
-
-          diss_r(k) = - ABS( u_comp(k) ) * (                                  &
-                      10.0_wp * ( w(k,j,i+1) - w(k,j,i)   )                   &
-                    -  5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) )                   &
-                    +           ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
-
-          v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
-          flux_n(k) = v_comp(k) * (                                           &
-                      37.0_wp * ( w(k,j+1,i) + w(k,j,i)   )                   &
-                    -  8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) )                   &
-                    +           ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
-
-          diss_n(k) = - ABS( v_comp(k) ) * (                                  &
-                      10.0_wp * ( w(k,j+1,i) - w(k,j,i)   )                   &
-                    -  5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) )                   &
-                    +           ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
-!
-!--       k index has to be modified near bottom and top, else array
-!--       subscripts will be exceeded.
-          ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
-          ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
-          ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
-
-          k_ppp = k + 3 * ibit35
-          k_pp  = k + 2 * ( 1 - ibit33  )
-          k_mm  = k - 2 * ibit35
-
-          w_comp(k) = w(k+1,j,i) + w(k,j,i)
-          flux_t(k) = w_comp(k) * rho_air(k+1) * (                            &
-                     ( 37.0_wp * ibit35 * adv_mom_5                           &
-                  +     7.0_wp * ibit34 * adv_mom_3                           &
-                  +              ibit33 * adv_mom_1                           &
-                     ) *                                                      &
-                                ( w(k+1,j,i)  + w(k,j,i)     )                &
-              -      (  8.0_wp * ibit35 * adv_mom_5                           &
-                  +              ibit34 * adv_mom_3                           &
-                     ) *                                                      &
-                                ( w(k_pp,j,i)  + w(k-1,j,i)  )                &
-              +      (           ibit35 * adv_mom_5                           &
-                     ) *                                                      &
-                                ( w(k_ppp,j,i) + w(k_mm,j,i) )                &
-                                                 )
-
-          diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * (                   &
-                     ( 10.0_wp * ibit35 * adv_mom_5                           &
-                  +     3.0_wp * ibit34 * adv_mom_3                           &
-                  +              ibit33 * adv_mom_1                           &
-                     ) *                                                      &
-                                ( w(k+1,j,i)   - w(k,j,i)    )                &
-              -      (  5.0_wp * ibit35 * adv_mom_5                           &
-                  +              ibit34 * adv_mom_3                           &
-                     ) *                                                      &
-                                ( w(k_pp,j,i)  - w(k-1,j,i)  )                &
-              +      (           ibit35 * adv_mom_5                           &
-                     ) *                                                      &
-                                ( w(k_ppp,j,i) - w(k_mm,j,i) )                &
-                                                          )
-       ENDDO
-
-       DO  k = nzb_max+1, nzt
-
-          flux_d    = flux_t(k-1)
-          diss_d    = diss_t(k-1)
-!
-!--       Calculate the divergence of the velocity field. A respective
-!--       correction is needed to overcome numerical instabilities introduced
-!--       by a not sufficient reduction of divergences near topography.
-          div = ( ( u_comp(k) + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx      &
-              +   ( v_comp(k) + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy      &
-              +   ( w_comp(k)               * rho_air(k+1) -                  &
-                  ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k)                      &
                   ) * drho_air_zw(k) * ddzu(k+1)                              &
                 ) * 0.5_wp
@@ -3260,13 +3299,4 @@
        
        INTEGER(iwp) ::  i      !< grid index along x-direction
-       INTEGER(iwp) ::  ibit0  !< flag indicating 1st-order scheme along x-direction
-       INTEGER(iwp) ::  ibit1  !< flag indicating 3rd-order scheme along x-direction
-       INTEGER(iwp) ::  ibit2  !< flag indicating 5th-order scheme along x-direction
-       INTEGER(iwp) ::  ibit3  !< flag indicating 1st-order scheme along y-direction
-       INTEGER(iwp) ::  ibit4  !< flag indicating 3rd-order scheme along y-direction
-       INTEGER(iwp) ::  ibit5  !< flag indicating 5th-order scheme along y-direction
-       INTEGER(iwp) ::  ibit6  !< flag indicating 1st-order scheme along z-direction
-       INTEGER(iwp) ::  ibit7  !< flag indicating 3rd-order scheme along z-direction
-       INTEGER(iwp) ::  ibit8  !< flag indicating 5th-order scheme along z-direction
        INTEGER(iwp) ::  j      !< grid index along y-direction
        INTEGER(iwp) ::  k      !< grid index along z-direction
@@ -3276,10 +3306,21 @@
        INTEGER(iwp) ::  tn = 0 !< number of OpenMP thread
        
-#if defined( __nopointer )
-       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !< advected scalar
-#else
-       REAL(wp), DIMENSION(:,:,:), POINTER ::  sk !< advected scalar
-#endif
-
+!       
+!--    sk is an array from parameter list. It should not be a pointer, because 
+!--    in that case the compiler can not assume a stride 1 and cannot perform 
+!--    a strided one vector load. Adding the CONTIGUOUS keyword makes things 
+!--    even worse, because the compiler cannot assume strided one in the 
+!--    caller side.
+       REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !<  advected scalar
+
+       REAL(wp) ::  ibit0  !< flag indicating 1st-order scheme along x-direction
+       REAL(wp) ::  ibit1  !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp) ::  ibit2  !< flag indicating 5th-order scheme along x-direction
+       REAL(wp) ::  ibit3  !< flag indicating 1st-order scheme along y-direction
+       REAL(wp) ::  ibit4  !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp) ::  ibit5  !< flag indicating 5th-order scheme along y-direction
+       REAL(wp) ::  ibit6  !< flag indicating 1st-order scheme along z-direction
+       REAL(wp) ::  ibit7  !< flag indicating 3rd-order scheme along z-direction
+       REAL(wp) ::  ibit8  !< flag indicating 5th-order scheme along z-direction
        REAL(wp) ::  diss_d !< artificial dissipation term at grid box bottom 
        REAL(wp) ::  div    !< diverence on scalar grid
@@ -3309,7 +3350,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit2 = IBITS(advc_flags_1(k,j,i-1),2,1)
-             ibit1 = IBITS(advc_flags_1(k,j,i-1),1,1)
-             ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
+             ibit2 = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )
+             ibit1 = REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )
+             ibit0 = REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )
 
              u_comp                 = u(k,j,i) - u_gtrans + u_stokes_zu(k)
@@ -3370,7 +3411,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit5 = IBITS(advc_flags_1(k,j-1,i),5,1)
-             ibit4 = IBITS(advc_flags_1(k,j-1,i),4,1)
-             ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
+             ibit5 = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )
+             ibit4 = REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )
+             ibit3 = REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )
 
              v_comp               = v(k,j,i) - v_gtrans + v_stokes_zu(k)
@@ -3433,7 +3474,7 @@
              DO  k = nzb+1, nzb_max
 
-                ibit2 = IBITS(advc_flags_1(k,j,i),2,1)
-                ibit1 = IBITS(advc_flags_1(k,j,i),1,1)
-                ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
+                ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp )
+                ibit1 = REAL( IBITS(advc_flags_1(k,j,i),1,1), KIND = wp )
+                ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp )
 
                 u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
@@ -3468,7 +3509,7 @@
                                              )
 
-                ibit5 = IBITS(advc_flags_1(k,j,i),5,1)
-                ibit4 = IBITS(advc_flags_1(k,j,i),4,1)
-                ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
+                ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp )
+                ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp )
+                ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp )
 
                 v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
@@ -3505,7 +3546,7 @@
 !--             k index has to be modified near bottom and top, else array
 !--             subscripts will be exceeded.
-                ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
-                ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
-                ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
+                ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
+                ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
+                ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
 
                 k_ppp = k + 3 * ibit8
@@ -3547,13 +3588,15 @@
 !--             by a not sufficient reduction of divergences near topography.
                 div   =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
-                          - u(k,j,i)   * ( IBITS(advc_flags_1(k,j,i-1),0,1)    &
-                                         + IBITS(advc_flags_1(k,j,i-1),1,1)    &
-                                         + IBITS(advc_flags_1(k,j,i-1),2,1)    &
+                          - u(k,j,i)   * (                                     &
+                        REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )    &
+                      + REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )    &
+                      + REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )    &
                                          )                                     &
                           ) * ddx                                              &
                         + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
-                          - v(k,j,i)   * ( IBITS(advc_flags_1(k,j-1,i),3,1)    &
-                                         + IBITS(advc_flags_1(k,j-1,i),4,1)    &
-                                         + IBITS(advc_flags_1(k,j-1,i),5,1)    &
+                          - v(k,j,i)   * (                                     &
+                        REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )    &
+                      + REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )    &
+                      + REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )    &
                                          )                                     &
                           ) * ddy                                              &
@@ -3561,7 +3604,8 @@
                                          ( ibit6 + ibit7 + ibit8 )             &
                           - w(k-1,j,i) * rho_air_zw(k-1) *                     &
-                                         ( IBITS(advc_flags_1(k-1,j,i),6,1)    &
-                                         + IBITS(advc_flags_1(k-1,j,i),7,1)    &
-                                         + IBITS(advc_flags_1(k-1,j,i),8,1)    &
+                                         (                                     &
+                        REAL( IBITS(advc_flags_1(k-1,j,i),6,1), KIND = wp )    &
+                      + REAL( IBITS(advc_flags_1(k-1,j,i),7,1), KIND = wp )    &
+                      + REAL( IBITS(advc_flags_1(k-1,j,i),8,1), KIND = wp )    &
                                          )                                     &      
                           ) * drho_air(k) * ddzw(k)
@@ -3611,7 +3655,7 @@
 !--             k index has to be modified near bottom and top, else array
 !--             subscripts will be exceeded.
-                ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
-                ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
-                ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
+                ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp )
+                ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp )
+                ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp )
 
                 k_ppp = k + 3 * ibit8
@@ -3817,13 +3861,4 @@
 
        INTEGER(iwp) ::  i      !< grid index along x-direction
-       INTEGER(iwp) ::  ibit9  !< flag indicating 1st-order scheme along x-direction
-       INTEGER(iwp) ::  ibit10 !< flag indicating 3rd-order scheme along x-direction
-       INTEGER(iwp) ::  ibit11 !< flag indicating 5th-order scheme along x-direction
-       INTEGER(iwp) ::  ibit12 !< flag indicating 1st-order scheme along y-direction
-       INTEGER(iwp) ::  ibit13 !< flag indicating 3rd-order scheme along y-direction
-       INTEGER(iwp) ::  ibit14 !< flag indicating 5th-order scheme along y-direction
-       INTEGER(iwp) ::  ibit15 !< flag indicating 1st-order scheme along z-direction
-       INTEGER(iwp) ::  ibit16 !< flag indicating 3rd-order scheme along z-direction
-       INTEGER(iwp) ::  ibit17 !< flag indicating 5th-order scheme along z-direction
        INTEGER(iwp) ::  j      !< grid index along y-direction
        INTEGER(iwp) ::  k      !< grid index along z-direction
@@ -3833,4 +3868,13 @@
        INTEGER(iwp) ::  tn = 0 !< number of OpenMP thread
        
+       REAL(wp)    ::  ibit9  !< flag indicating 1st-order scheme along x-direction
+       REAL(wp)    ::  ibit10 !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp)    ::  ibit11 !< flag indicating 5th-order scheme along x-direction
+       REAL(wp)    ::  ibit12 !< flag indicating 1st-order scheme along y-direction
+       REAL(wp)    ::  ibit13 !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp)    ::  ibit14 !< flag indicating 5th-order scheme along y-direction
+       REAL(wp)    ::  ibit15 !< flag indicating 1st-order scheme along z-direction
+       REAL(wp)    ::  ibit16 !< flag indicating 3rd-order scheme along z-direction
+       REAL(wp)    ::  ibit17 !< flag indicating 5th-order scheme along z-direction
        REAL(wp)    ::  diss_d !< artificial dissipation term at grid box bottom
        REAL(wp)    ::  div    !< diverence on u-grid
@@ -3864,20 +3908,20 @@
           DO  k = nzb+1, nzb_max
 
-             ibit11 = IBITS(advc_flags_1(k,j,i-1),11,1)
-             ibit10 = IBITS(advc_flags_1(k,j,i-1),10,1)
-             ibit9  = IBITS(advc_flags_1(k,j,i-1),9,1)
+             ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )
+             ibit10 = REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )
+             ibit9  = REAL( IBITS(advc_flags_1(k,j,i-1),9,1),  KIND = wp )
 
              u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
              swap_flux_x_local_u(k,j) = u_comp(k) * (                          &
-                                       ( 37.0_wp * ibit11 * adv_mom_5             &
-                                    +     7.0_wp * ibit10 * adv_mom_3             &
-                                    +              ibit9  * adv_mom_1             &
+                                       ( 37.0_wp * ibit11 * adv_mom_5          &
+                                    +     7.0_wp * ibit10 * adv_mom_3          &
+                                    +              ibit9  * adv_mom_1          &
                                        ) *                                     &
                                      ( u(k,j,i)   + u(k,j,i-1) )               &
-                                -      (  8.0_wp * ibit11 * adv_mom_5             &
-                                    +              ibit10 * adv_mom_3             &
+                                -      (  8.0_wp * ibit11 * adv_mom_5          &
+                                    +              ibit10 * adv_mom_3          &
                                        ) *                                     &
                                      ( u(k,j,i+1) + u(k,j,i-2) )               &
-                                +      (           ibit11 * adv_mom_5             &
+                                +      (           ibit11 * adv_mom_5          &
                                        ) *                                     &
                                      ( u(k,j,i+2) + u(k,j,i-3) )               &
@@ -3885,14 +3929,14 @@
 
               swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * (                &
-                                       ( 10.0_wp * ibit11 * adv_mom_5             &
-                                    +     3.0_wp * ibit10 * adv_mom_3             &
-                                    +              ibit9  * adv_mom_1             &
+                                       ( 10.0_wp * ibit11 * adv_mom_5          &
+                                    +     3.0_wp * ibit10 * adv_mom_3          &
+                                    +              ibit9  * adv_mom_1          &
                                        ) *                                     &
                                      ( u(k,j,i)   - u(k,j,i-1) )               &
-                                -      (  5.0_wp * ibit11 * adv_mom_5             &
-                                    +              ibit10 * adv_mom_3             &
+                                -      (  5.0_wp * ibit11 * adv_mom_5          &
+                                    +              ibit10 * adv_mom_3          &
                                        ) *                                     &
                                      ( u(k,j,i+1) - u(k,j,i-2) )               &
-                                +      (           ibit11 * adv_mom_5             &
+                                +      (           ibit11 * adv_mom_5          &
                                        ) *                                     &
                                      ( u(k,j,i+2) - u(k,j,i-3) )               &
@@ -3905,10 +3949,10 @@
              u_comp(k)         = u(k,j,i) + u(k,j,i-1) - gu
              swap_flux_x_local_u(k,j) = u_comp(k) * (                          &
-                             37.0_wp * ( u(k,j,i) + u(k,j,i-1)   )                &
-                           -  8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) )                &
+                             37.0_wp * ( u(k,j,i) + u(k,j,i-1)   )             &
+                           -  8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) )             &
                            +           ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
              swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * (                   &
-                             10.0_wp * ( u(k,j,i) - u(k,j,i-1)   )                &
-                           -  5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) )                &
+                             10.0_wp * ( u(k,j,i) - u(k,j,i-1)   )             &
+                           -  5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) )             &
                            +           ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
 
@@ -3922,20 +3966,20 @@
           DO  k = nzb+1, nzb_max
 
-             ibit14 = IBITS(advc_flags_1(k,j-1,i),14,1)
-             ibit13 = IBITS(advc_flags_1(k,j-1,i),13,1)
-             ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)
+             ibit14 = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )
+             ibit13 = REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp )
+             ibit12 = REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp )
 
              v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
              swap_flux_y_local_u(k) = v_comp * (                              &
-                                   ( 37.0_wp * ibit14 * adv_mom_5                &
-                                +     7.0_wp * ibit13 * adv_mom_3                &
-                                +              ibit12 * adv_mom_1                &
+                                   ( 37.0_wp * ibit14 * adv_mom_5             &
+                                +     7.0_wp * ibit13 * adv_mom_3             &
+                                +              ibit12 * adv_mom_1             &
                                    ) *                                        &
                                      ( u(k,j,i)   + u(k,j-1,i) )              &
-                            -      (  8.0_wp * ibit14 * adv_mom_5                &
-                            +                  ibit13 * adv_mom_3                    &
+                            -      (  8.0_wp * ibit14 * adv_mom_5             &
+                            +                  ibit13 * adv_mom_3             &
                                    ) *                                        &
                                      ( u(k,j+1,i) + u(k,j-2,i) )              &
-                        +      (               ibit14 * adv_mom_5                    &
+                        +      (               ibit14 * adv_mom_5             &
                                ) *                                            &
                                      ( u(k,j+2,i) + u(k,j-3,i) )              &
@@ -3943,14 +3987,14 @@
 
              swap_diss_y_local_u(k) = - ABS ( v_comp ) * (                    &
-                                   ( 10.0_wp * ibit14 * adv_mom_5                &
-                                +     3.0_wp * ibit13 * adv_mom_3                &
-                                +              ibit12 * adv_mom_1                &
+                                   ( 10.0_wp * ibit14 * adv_mom_5             &
+                                +     3.0_wp * ibit13 * adv_mom_3             &
+                                +              ibit12 * adv_mom_1             &
                                    ) *                                        &
                                      ( u(k,j,i)   - u(k,j-1,i) )              &
-                            -      (  5.0_wp * ibit14 * adv_mom_5                &
-                                +              ibit13 * adv_mom_3                &
+                            -      (  5.0_wp * ibit14 * adv_mom_5             &
+                                +              ibit13 * adv_mom_3             &
                                    ) *                                        &
                                      ( u(k,j+1,i) - u(k,j-2,i) )              &
-                            +      (           ibit14 * adv_mom_5                &
+                            +      (           ibit14 * adv_mom_5             &
                                    ) *                                        &
                                      ( u(k,j+2,i) - u(k,j-3,i) )              &
@@ -3983,20 +4027,20 @@
              DO  k = nzb+1, nzb_max
 
-                ibit11 = IBITS(advc_flags_1(k,j,i),11,1)
-                ibit10 = IBITS(advc_flags_1(k,j,i),10,1)
-                ibit9  = IBITS(advc_flags_1(k,j,i),9,1)
+                ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp )
+                ibit10 = REAL( IBITS(advc_flags_1(k,j,i),10,1), KIND = wp )
+                ibit9  = REAL( IBITS(advc_flags_1(k,j,i),9,1),  KIND = wp )
 
                 u_comp(k) = u(k,j,i+1) + u(k,j,i)
                 flux_r(k) = ( u_comp(k) - gu ) * (                           &
-                          ( 37.0_wp * ibit11 * adv_mom_5                        &
-                       +     7.0_wp * ibit10 * adv_mom_3                        &
-                       +              ibit9  * adv_mom_1                        &
+                          ( 37.0_wp * ibit11 * adv_mom_5                     &
+                       +     7.0_wp * ibit10 * adv_mom_3                     &
+                       +              ibit9  * adv_mom_1                     &
                           ) *                                                &
                                  ( u(k,j,i+1) + u(k,j,i)   )                 &
-                   -      (  8.0_wp * ibit11 * adv_mom_5                        &
-                       +              ibit10 * adv_mom_3                        &
+                   -      (  8.0_wp * ibit11 * adv_mom_5                     &
+                       +              ibit10 * adv_mom_3                     &
                           ) *                                                &
                                  ( u(k,j,i+2) + u(k,j,i-1) )                 &
-                   +      (           ibit11 * adv_mom_5                        &
+                   +      (           ibit11 * adv_mom_5                     &
                           ) *                                                &
                                  ( u(k,j,i+3) + u(k,j,i-2) )                 &
@@ -4004,34 +4048,34 @@
 
                 diss_r(k) = - ABS( u_comp(k) - gu ) * (                      &
-                          ( 10.0_wp * ibit11 * adv_mom_5                        &
-                       +     3.0_wp * ibit10 * adv_mom_3                        & 
-                       +              ibit9  * adv_mom_1                        &
+                          ( 10.0_wp * ibit11 * adv_mom_5                     &
+                       +     3.0_wp * ibit10 * adv_mom_3                     & 
+                       +              ibit9  * adv_mom_1                     &
                           ) *                                                &
                                  ( u(k,j,i+1) - u(k,j,i)  )                  &
-                   -      (  5.0_wp * ibit11 * adv_mom_5                        &
-                       +              ibit10 * adv_mom_3                        &
+                   -      (  5.0_wp * ibit11 * adv_mom_5                     &
+                       +              ibit10 * adv_mom_3                     &
                           ) *                                                &
                                  ( u(k,j,i+2) - u(k,j,i-1) )                 &
-                   +      (           ibit11 * adv_mom_5                        &
+                   +      (           ibit11 * adv_mom_5                     &
                           ) *                                                &
                                  ( u(k,j,i+3) - u(k,j,i-2) )                 &
                                                      )
 
-                ibit14 = IBITS(advc_flags_1(k,j,i),14,1)
-                ibit13 = IBITS(advc_flags_1(k,j,i),13,1)
-                ibit12 = IBITS(advc_flags_1(k,j,i),12,1)
+                ibit14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp )
+                ibit13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp )
+                ibit12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp )
 
                 v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
                 flux_n(k) = v_comp * (                                       &
-                          ( 37.0_wp * ibit14 * adv_mom_5                        &
-                       +     7.0_wp * ibit13 * adv_mom_3                        &
-                       +              ibit12 * adv_mom_1                        &
+                          ( 37.0_wp * ibit14 * adv_mom_5                     &
+                       +     7.0_wp * ibit13 * adv_mom_3                     &
+                       +              ibit12 * adv_mom_1                     &
                           ) *                                                &
                                  ( u(k,j+1,i) + u(k,j,i)   )                 &
-                   -      (  8.0_wp * ibit14 * adv_mom_5                        &
-                       +              ibit13 * adv_mom_3                        &
+                   -      (  8.0_wp * ibit14 * adv_mom_5                     &
+                       +              ibit13 * adv_mom_3                     &
                           ) *                                                &
                                  ( u(k,j+2,i) + u(k,j-1,i) )                 &
-                   +      (           ibit14 * adv_mom_5                        & 
+                   +      (           ibit14 * adv_mom_5                     & 
                           ) *                                                &
                                  ( u(k,j+3,i) + u(k,j-2,i) )                 &
@@ -4039,14 +4083,14 @@
 
                 diss_n(k) = - ABS ( v_comp ) * (                             &
-                          ( 10.0_wp * ibit14 * adv_mom_5                        &
-                       +     3.0_wp * ibit13 * adv_mom_3                        &
-                       +              ibit12 * adv_mom_1                        &
+                          ( 10.0_wp * ibit14 * adv_mom_5                     &
+                       +     3.0_wp * ibit13 * adv_mom_3                     &
+                       +              ibit12 * adv_mom_1                     &
                           ) *                                                &
                                  ( u(k,j+1,i) - u(k,j,i)  )                  &
-                   -      (  5.0_wp * ibit14 * adv_mom_5                        &
-                       +              ibit13 * adv_mom_3                        &
+                   -      (  5.0_wp * ibit14 * adv_mom_5                     &
+                       +              ibit13 * adv_mom_3                     &
                           ) *                                                &
                                  ( u(k,j+2,i) - u(k,j-1,i) )                 &
-                   +      (           ibit14 * adv_mom_5                        &
+                   +      (           ibit14 * adv_mom_5                     &
                           ) *                                                &
                                  ( u(k,j+3,i) - u(k,j-2,i) )                 &
@@ -4055,7 +4099,7 @@
 !--             k index has to be modified near bottom and top, else array
 !--             subscripts will be exceeded.
-                ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
-                ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
-                ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
+                ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
+                ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
+                ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
 
                 k_ppp = k + 3 * ibit17
@@ -4065,14 +4109,14 @@
                 w_comp    = w(k,j,i) + w(k,j,i-1)
                 flux_t(k) = w_comp * rho_air_zw(k) * (                       &
-                          ( 37.0_wp * ibit17 * adv_mom_5                        &
-                       +     7.0_wp * ibit16 * adv_mom_3                        &
-                       +              ibit15 * adv_mom_1                        & 
+                          ( 37.0_wp * ibit17 * adv_mom_5                     &
+                       +     7.0_wp * ibit16 * adv_mom_3                     &
+                       +              ibit15 * adv_mom_1                     & 
                           ) *                                                &
                              ( u(k+1,j,i)  + u(k,j,i)     )                  &
-                   -      (  8.0_wp * ibit17 * adv_mom_5                        &
-                       +              ibit16 * adv_mom_3                        &
+                   -      (  8.0_wp * ibit17 * adv_mom_5                     &
+                       +              ibit16 * adv_mom_3                     &
                           ) *                                                &
                              ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
-                   +      (           ibit17 * adv_mom_5                        &
+                   +      (           ibit17 * adv_mom_5                     &
                           ) *                                                &
                              ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
@@ -4080,14 +4124,14 @@
 
                 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
-                          ( 10.0_wp * ibit17 * adv_mom_5                        &
-                       +     3.0_wp * ibit16 * adv_mom_3                        &
-                       +              ibit15 * adv_mom_1                        &
+                          ( 10.0_wp * ibit17 * adv_mom_5                     &
+                       +     3.0_wp * ibit16 * adv_mom_3                     &
+                       +              ibit15 * adv_mom_1                     &
                           ) *                                                &
                              ( u(k+1,j,i)   - u(k,j,i)    )                  &
-                   -      (  5.0_wp * ibit17 * adv_mom_5                        &
-                       +              ibit16 * adv_mom_3                        &
+                   -      (  5.0_wp * ibit17 * adv_mom_5                     &
+                       +              ibit16 * adv_mom_3                     &
                           ) *                                                &
                              ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
-                   +      (           ibit17 * adv_mom_5                        &
+                   +      (           ibit17 * adv_mom_5                     &
                            ) *                                               &
                              ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
@@ -4099,21 +4143,24 @@
                 div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 )             &
                 - ( u(k,j,i)   + u(k,j,i-1)   )                               &
-                                    * ( IBITS(advc_flags_1(k,j,i-1),9,1)      &
-                                      + IBITS(advc_flags_1(k,j,i-1),10,1)     &
-                                      + IBITS(advc_flags_1(k,j,i-1),11,1)     &
+                                    * (                                       &
+                    REAL( IBITS(advc_flags_1(k,j,i-1),9,1),  KIND = wp )      &
+                  + REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )      &
+                  + REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )      &
                                       )                                       &
                   ) * ddx                                                     &
                +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
                   - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
-                                    * ( IBITS(advc_flags_1(k,j-1,i),12,1)     &
-                                      + IBITS(advc_flags_1(k,j-1,i),13,1)     &
-                                      + IBITS(advc_flags_1(k,j-1,i),14,1)     &
+                                    * (                                       &
+                     REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )     &
                                       )                                       &
                   ) * ddy                                                     &
                +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
                 - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
-                                    * ( IBITS(advc_flags_1(k-1,j,i),15,1)     &
-                                      + IBITS(advc_flags_1(k-1,j,i),16,1)     &
-                                      + IBITS(advc_flags_1(k-1,j,i),17,1)     &
+                                    * (                                       &
+                     REAL( IBITS(advc_flags_1(k-1,j,i),15,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k-1,j,i),16,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k-1,j,i),17,1), KIND = wp )     &
                                       )                                       &  
                   ) * drho_air(k) * ddzw(k)                                   &
@@ -4186,7 +4233,7 @@
 !--             k index has to be modified near bottom and top, else array
 !--             subscripts will be exceeded.
-                ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
-                ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
-                ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
+                ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )
+                ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )
+                ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )
 
                 k_ppp = k + 3 * ibit17
@@ -4309,13 +4356,4 @@
 
        INTEGER(iwp) ::  i      !< grid index along x-direction
-       INTEGER(iwp) ::  ibit18 !< flag indicating 1st-order scheme along x-direction
-       INTEGER(iwp) ::  ibit19 !< flag indicating 3rd-order scheme along x-direction
-       INTEGER(iwp) ::  ibit20 !< flag indicating 5th-order scheme along x-direction
-       INTEGER(iwp) ::  ibit21 !< flag indicating 1st-order scheme along y-direction
-       INTEGER(iwp) ::  ibit22 !< flag indicating 3rd-order scheme along y-direction
-       INTEGER(iwp) ::  ibit23 !< flag indicating 5th-order scheme along y-direction
-       INTEGER(iwp) ::  ibit24 !< flag indicating 1st-order scheme along z-direction
-       INTEGER(iwp) ::  ibit25 !< flag indicating 3rd-order scheme along z-direction
-       INTEGER(iwp) ::  ibit26 !< flag indicating 5th-order scheme along z-direction
        INTEGER(iwp) ::  j      !< grid index along y-direction
        INTEGER(iwp) ::  k      !< grid index along z-direction
@@ -4325,4 +4363,13 @@
        INTEGER(iwp) ::  tn = 0 !< number of OpenMP thread
        
+       REAL(wp)    ::  ibit18 !< flag indicating 1st-order scheme along x-direction
+       REAL(wp)    ::  ibit19 !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp)    ::  ibit20 !< flag indicating 5th-order scheme along x-direction
+       REAL(wp)    ::  ibit21 !< flag indicating 1st-order scheme along y-direction
+       REAL(wp)    ::  ibit22 !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp)    ::  ibit23 !< flag indicating 5th-order scheme along y-direction
+       REAL(wp)    ::  ibit24 !< flag indicating 1st-order scheme along z-direction
+       REAL(wp)    ::  ibit25 !< flag indicating 3rd-order scheme along z-direction
+       REAL(wp)    ::  ibit26 !< flag indicating 5th-order scheme along z-direction
        REAL(wp)    ::  diss_d !< artificial dissipation term at grid box bottom
        REAL(wp)    ::  div    !< diverence on v-grid
@@ -4355,7 +4402,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit20 = IBITS(advc_flags_1(k,j,i-1),20,1)
-             ibit19 = IBITS(advc_flags_1(k,j,i-1),19,1)
-             ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)
+             ibit20 = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )
+             ibit19 = REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp )
+             ibit18 = REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp )
 
              u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
@@ -4413,7 +4460,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit23 = IBITS(advc_flags_1(k,j-1,i),23,1)
-             ibit22 = IBITS(advc_flags_1(k,j-1,i),22,1)
-             ibit21 = IBITS(advc_flags_1(k,j-1,i),21,1)
+             ibit23 = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )
+             ibit22 = REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp )
+             ibit21 = REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )
 
              v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
@@ -4473,7 +4520,7 @@
              DO  k = nzb+1, nzb_max
 
-                ibit20 = IBITS(advc_flags_1(k,j,i),20,1)
-                ibit19 = IBITS(advc_flags_1(k,j,i),19,1)
-                ibit18 = IBITS(advc_flags_1(k,j,i),18,1)
+                ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp )
+                ibit19 = REAL( IBITS(advc_flags_1(k,j,i),19,1), KIND = wp )
+                ibit18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp )
 
                 u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
@@ -4508,7 +4555,7 @@
                                               )
 
-                ibit23 = IBITS(advc_flags_1(k,j,i),23,1)
-                ibit22 = IBITS(advc_flags_1(k,j,i),22,1)
-                ibit21 = IBITS(advc_flags_1(k,j,i),21,1)
+                ibit23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp )
+                ibit22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp )
+                ibit21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp )
 
                 v_comp(k) = v(k,j+1,i) + v(k,j,i)
@@ -4545,7 +4592,7 @@
 !--             k index has to be modified near bottom and top, else array
 !--             subscripts will be exceeded.
-                ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
-                ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
-                ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
+                ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
+                ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
+                ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
 
                 k_ppp = k + 3 * ibit26
@@ -4590,7 +4637,8 @@
                                        * ( ibit18 + ibit19 + ibit20 )         &
                 - ( u(k,j-1,i)   + u(k,j,i) )                                 &
-                                       * ( IBITS(advc_flags_1(k,j,i-1),18,1)  &
-                                         + IBITS(advc_flags_1(k,j,i-1),19,1)  &
-                                         + IBITS(advc_flags_1(k,j,i-1),20,1)  &
+                                       * (                                    &
+                        REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp )  &
+                      + REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp )  &
+                      + REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )  &
                                          )                                    &
                   ) * ddx                                                     &
@@ -4598,7 +4646,8 @@
                                        * ( ibit21 + ibit22 + ibit23 )         &
                 - ( v(k,j,i)     + v(k,j-1,i) )                               &
-                                       * ( IBITS(advc_flags_1(k,j-1,i),21,1)  &
-                                         + IBITS(advc_flags_1(k,j-1,i),22,1)  &
-                                         + IBITS(advc_flags_1(k,j-1,i),23,1)  &
+                                       * (                                    &
+                        REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )  &
+                      + REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp )  &
+                      + REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )  &
                                          )                                    &
                   ) * ddy                                                     &
@@ -4606,7 +4655,8 @@
                                        * ( ibit24 + ibit25 + ibit26 )         &
                 - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
-                                       * ( IBITS(advc_flags_1(k-1,j,i),24,1)  &
-                                         + IBITS(advc_flags_1(k-1,j,i),25,1)  &
-                                         + IBITS(advc_flags_1(k-1,j,i),26,1)  &
+                                       * (                                    &
+                        REAL( IBITS(advc_flags_1(k-1,j,i),24,1), KIND = wp )  &
+                      + REAL( IBITS(advc_flags_1(k-1,j,i),25,1), KIND = wp )  &
+                      + REAL( IBITS(advc_flags_1(k-1,j,i),26,1), KIND = wp )  &
                                          )                                    &
                    ) * drho_air(k) * ddzw(k)                                  &
@@ -4684,7 +4734,7 @@
 !--             k index has to be modified near bottom and top, else array
 !--             subscripts will be exceeded.
-                ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
-                ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
-                ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
+                ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )
+                ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )
+                ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )
 
                 k_ppp = k + 3 * ibit26
@@ -4811,13 +4861,4 @@
 
        INTEGER(iwp) ::  i      !< grid index along x-direction
-       INTEGER(iwp) ::  ibit27 !< flag indicating 1st-order scheme along x-direction
-       INTEGER(iwp) ::  ibit28 !< flag indicating 3rd-order scheme along x-direction
-       INTEGER(iwp) ::  ibit29 !< flag indicating 5th-order scheme along x-direction
-       INTEGER(iwp) ::  ibit30 !< flag indicating 1st-order scheme along y-direction
-       INTEGER(iwp) ::  ibit31 !< flag indicating 3rd-order scheme along y-direction
-       INTEGER(iwp) ::  ibit32 !< flag indicating 5th-order scheme along y-direction
-       INTEGER(iwp) ::  ibit33 !< flag indicating 1st-order scheme along z-direction
-       INTEGER(iwp) ::  ibit34 !< flag indicating 3rd-order scheme along z-direction
-       INTEGER(iwp) ::  ibit35 !< flag indicating 5th-order scheme along z-direction
        INTEGER(iwp) ::  j      !< grid index along y-direction
        INTEGER(iwp) ::  k      !< grid index along z-direction
@@ -4827,4 +4868,13 @@
        INTEGER(iwp) ::  tn = 0 !< number of OpenMP thread
        
+       REAL(wp)    ::  ibit27 !< flag indicating 1st-order scheme along x-direction
+       REAL(wp)    ::  ibit28 !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp)    ::  ibit29 !< flag indicating 5th-order scheme along x-direction
+       REAL(wp)    ::  ibit30 !< flag indicating 1st-order scheme along y-direction
+       REAL(wp)    ::  ibit31 !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp)    ::  ibit32 !< flag indicating 5th-order scheme along y-direction
+       REAL(wp)    ::  ibit33 !< flag indicating 1st-order scheme along z-direction
+       REAL(wp)    ::  ibit34 !< flag indicating 3rd-order scheme along z-direction
+       REAL(wp)    ::  ibit35 !< flag indicating 5th-order scheme along z-direction
        REAL(wp)    ::  diss_d !< artificial dissipation term at grid box bottom
        REAL(wp)    ::  div    !< divergence on w-grid
@@ -4857,7 +4907,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit29 = IBITS(advc_flags_1(k,j,i-1),29,1)
-             ibit28 = IBITS(advc_flags_1(k,j,i-1),28,1)
-             ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)
+             ibit29 = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )
+             ibit28 = REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp )
+             ibit27 = REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp )
 
              u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
@@ -4915,7 +4965,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit32 = IBITS(advc_flags_2(k,j-1,i),0,1)
-             ibit31 = IBITS(advc_flags_1(k,j-1,i),31,1)
-             ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)
+             ibit32 = REAL( IBITS(advc_flags_2(k,j-1,i),0,1),  KIND = wp )
+             ibit31 = REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )
+             ibit30 = REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp )
 
              v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
@@ -4981,7 +5031,7 @@
              DO  k = nzb+1, nzb_max
 
-                ibit29 = IBITS(advc_flags_1(k,j,i),29,1)
-                ibit28 = IBITS(advc_flags_1(k,j,i),28,1)
-                ibit27 = IBITS(advc_flags_1(k,j,i),27,1)
+                ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp )
+                ibit28 = REAL( IBITS(advc_flags_1(k,j,i),28,1), KIND = wp )
+                ibit27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp )
 
                 u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
@@ -5016,7 +5066,7 @@
                                               )
 
-                ibit32 = IBITS(advc_flags_2(k,j,i),0,1)
-                ibit31 = IBITS(advc_flags_1(k,j,i),31,1)
-                ibit30 = IBITS(advc_flags_1(k,j,i),30,1)
+                ibit32 = REAL( IBITS(advc_flags_2(k,j,i),0,1),  KIND = wp )
+                ibit31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp )
+                ibit30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp )
 
                 v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
@@ -5053,7 +5103,7 @@
 !--             k index has to be modified near bottom and top, else array
 !--             subscripts will be exceeded.
-                ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
-                ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
-                ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
+                ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
+                ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
+                ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
 
                 k_ppp = k + 3 * ibit35
@@ -5097,21 +5147,24 @@
                 div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 )      &
                   - ( u(k+1,j,i) + u(k,j,i)   )                               & 
-                                    * ( IBITS(advc_flags_1(k,j,i-1),27,1)     &
-                                      + IBITS(advc_flags_1(k,j,i-1),28,1)     &
-                                      + IBITS(advc_flags_1(k,j,i-1),29,1)     &
+                                    * (                                       &
+                     REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )     &
                                       )                                       &
                   ) * ddx                                                     &
               +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
                   - ( v(k+1,j,i) + v(k,j,i)   )                               &
-                                    * ( IBITS(advc_flags_1(k,j-1,i),30,1)     &
-                                      + IBITS(advc_flags_1(k,j-1,i),31,1)     &
-                                      + IBITS(advc_flags_2(k,j-1,i),0,1)      &
+                                    * (                                       &
+                     REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )     &
+                   + REAL( IBITS(advc_flags_2(k,j-1,i),0,1),  KIND = wp )     &
                                       )                                       &
                   ) * ddy                                                     &
               +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
                 - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
-                                    * ( IBITS(advc_flags_2(k-1,j,i),1,1)      &
-                                      + IBITS(advc_flags_2(k-1,j,i),2,1)      &
-                                      + IBITS(advc_flags_2(k-1,j,i),3,1)      &
+                                    * (                                       &
+                     REAL( IBITS(advc_flags_2(k-1,j,i),1,1), KIND = wp )      &
+                   + REAL( IBITS(advc_flags_2(k-1,j,i),2,1), KIND = wp )      &
+                   + REAL( IBITS(advc_flags_2(k-1,j,i),3,1), KIND = wp )      &
                                       )                                       & 
                   ) * drho_air_zw(k) * ddzu(k+1)                              &
@@ -5176,7 +5229,7 @@
 !--             k index has to be modified near bottom and top, else array
 !--             subscripts will be exceeded.
-                ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
-                ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
-                ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
+                ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )
+                ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )
+                ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )
 
                 k_ppp = k + 3 * ibit35
