Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 2101)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 2118)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! OpenACC version of subroutines removed
 ! 
 ! Former revisions:
@@ -181,7 +181,6 @@
 
     PRIVATE
-    PUBLIC   advec_s_ws, advec_s_ws_acc, advec_u_ws, advec_u_ws_acc,          &
-             advec_v_ws, advec_v_ws_acc, advec_w_ws, advec_w_ws_acc,          &
-             ws_init, ws_init_flags, ws_statistics
+    PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, ws_init,          &
+             ws_init_flags, ws_statistics
 
     INTERFACE ws_init
@@ -207,8 +206,4 @@
     END INTERFACE advec_u_ws
 
-    INTERFACE advec_u_ws_acc
-       MODULE PROCEDURE advec_u_ws_acc
-    END INTERFACE advec_u_ws_acc
-
     INTERFACE advec_v_ws
        MODULE PROCEDURE advec_v_ws
@@ -216,16 +211,8 @@
     END INTERFACE advec_v_ws
 
-    INTERFACE advec_v_ws_acc
-       MODULE PROCEDURE advec_v_ws_acc
-    END INTERFACE advec_v_ws_acc
-
     INTERFACE advec_w_ws
        MODULE PROCEDURE advec_w_ws
        MODULE PROCEDURE advec_w_ws_ij
     END INTERFACE advec_w_ws
-
-    INTERFACE advec_w_ws_acc
-       MODULE PROCEDURE advec_w_ws_acc
-    END INTERFACE advec_w_ws_acc
 
  CONTAINS
@@ -4029,514 +4016,4 @@
 
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Scalar advection - Call for all grid points - accelerator version
-!------------------------------------------------------------------------------!
-    SUBROUTINE advec_s_ws_acc ( sk, sk_char )
-
-       USE arrays_3d,                                                         &
-           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
-
-       USE constants,                                                         &
-           ONLY:  adv_sca_1, adv_sca_3, adv_sca_5
-
-       USE control_parameters,                                                &
-           ONLY:  intermediate_timestep_count, monotonic_adjustment, u_gtrans,&
-                  v_gtrans
-
-       USE grid_variables,                                                    &
-           ONLY:  ddx, ddy
-
-       USE indices,                                                           &
-           ONLY:  i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg,  &
-                  nzb, nzb_max, nzt, wall_flags_0
-
-       USE kinds
-       
-!        USE statistics,                                                       &
-!            ONLY:  sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,          &
-!                   sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep
-
-       IMPLICIT NONE
-
-       CHARACTER (LEN = *), INTENT(IN)    :: sk_char !<
-
-       INTEGER(iwp) ::  i      !<
-       INTEGER(iwp) ::  ibit0  !<
-       INTEGER(iwp) ::  ibit1  !<
-       INTEGER(iwp) ::  ibit2  !<
-       INTEGER(iwp) ::  ibit3  !<
-       INTEGER(iwp) ::  ibit4  !<
-       INTEGER(iwp) ::  ibit5  !<
-       INTEGER(iwp) ::  ibit6  !<
-       INTEGER(iwp) ::  ibit7  !<
-       INTEGER(iwp) ::  ibit8  !<
-       INTEGER(iwp) ::  j      !<
-       INTEGER(iwp) ::  k      !<
-       INTEGER(iwp) ::  k_mm   !<
-       INTEGER(iwp) ::  k_mmm  !<
-       INTEGER(iwp) ::  k_pp   !<
-       INTEGER(iwp) ::  k_ppp  !<
-       INTEGER(iwp) ::  tn = 0 !<
-
-       REAL(wp)    ::  diss_d !<
-       REAL(wp)    ::  diss_l !<
-       REAL(wp)    ::  diss_n !<
-       REAL(wp)    ::  diss_r !<
-       REAL(wp)    ::  diss_s !<
-       REAL(wp)    ::  diss_t !<
-       REAL(wp)    ::  div    !<
-       REAL(wp)    ::  flux_d !<
-       REAL(wp)    ::  flux_l !<
-       REAL(wp)    ::  flux_n !<
-       REAL(wp)    ::  flux_r !<
-       REAL(wp)    ::  flux_s !<
-       REAL(wp)    ::  flux_t !<
-       REAL(wp)    ::  fd_1   !<
-       REAL(wp)    ::  fl_1   !<
-       REAL(wp)    ::  fn_1   !<
-       REAL(wp)    ::  fr_1   !<
-       REAL(wp)    ::  fs_1   !<
-       REAL(wp)    ::  ft_1   !<
-       REAL(wp)    ::  phi_d  !<
-       REAL(wp)    ::  phi_l  !<
-       REAL(wp)    ::  phi_n  !<
-       REAL(wp)    ::  phi_r  !<
-       REAL(wp)    ::  phi_s  !<
-       REAL(wp)    ::  phi_t  !<
-       REAL(wp)    ::  rd     !<
-       REAL(wp)    ::  rl     !<
-       REAL(wp)    ::  rn     !<
-       REAL(wp)    ::  rr     !<
-       REAL(wp)    ::  rs     !<
-       REAL(wp)    ::  rt     !<
-       REAL(wp)    ::  u_comp !<
-       REAL(wp)    ::  v_comp !<
-
-       REAL(wp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  ::  sk !<
-
-!
-!--    Computation of fluxes and tendency terms
-       !$acc kernels present( ddzw, sk, tend, u, v, w, wall_flags_0 )
-       DO  i = i_left, i_right
-          DO  j = j_south, j_north
-             DO  k = nzb+1, nzt
-
-                ibit2 = IBITS(wall_flags_0(k,j,i-1),2,1)
-                ibit1 = IBITS(wall_flags_0(k,j,i-1),1,1)
-                ibit0 = IBITS(wall_flags_0(k,j,i-1),0,1)
-
-                u_comp              = u(k,j,i) - u_gtrans
-                flux_l              = u_comp * (                              &
-                                               ( 37.0_wp * ibit2 * adv_sca_5  &
-                                            +     7.0_wp * ibit1 * adv_sca_3  &
-                                            +              ibit0 * adv_sca_1  &
-                                               ) *                            &
-                                         ( sk(k,j,i)   + sk(k,j,i-1)    )     &
-                                        -      (  8.0_wp * ibit2 * adv_sca_5  &
-                                            +              ibit1 * adv_sca_3  &
-                                               ) *                            &
-                                         ( sk(k,j,i+1) + sk(k,j,i-2)    )     &
-                                        +      (           ibit2 * adv_sca_5  &
-                                               ) *                            &
-                                         ( sk(k,j,i+2) + sk(k,j,i-3)    )     &
-                                               )
-
-                diss_l              = -ABS( u_comp ) * (                      &
-                                               ( 10.0_wp * ibit2 * adv_sca_5  &
-                                            +     3.0_wp * ibit1 * adv_sca_3  &
-                                            +              ibit0 * adv_sca_1  &
-                                               ) *                            &
-                                         ( sk(k,j,i)   - sk(k,j,i-1)    )     &
-                                        -      (  5.0_wp * ibit2 * adv_sca_5  &
-                                            +              ibit1 * adv_sca_3  &
-                                               ) *                            &
-                                         ( sk(k,j,i+1) - sk(k,j,i-2)    )     &
-                                        +      (           ibit2 * adv_sca_5  &
-                                               ) *                            &
-                                         ( sk(k,j,i+2) - sk(k,j,i-3)    )     &
-                                                        )
-
-                ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
-                ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
-                ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
-
-                u_comp    = u(k,j,i+1) - u_gtrans
-                flux_r    = u_comp * (                                        &
-                          ( 37.0_wp * ibit2 * adv_sca_5                       &
-                      +      7.0_wp * ibit1 * adv_sca_3                       &
-                      +               ibit0 * adv_sca_1                       &
-                          ) *                                                 &
-                             ( sk(k,j,i+1) + sk(k,j,i)   )                    &
-                   -      (  8.0_wp * ibit2 * adv_sca_5                       &
-                       +              ibit1 * adv_sca_3                       &
-                          ) *                                                 &
-                             ( sk(k,j,i+2) + sk(k,j,i-1) )                    &
-                   +      (           ibit2 * adv_sca_5                       &
-                          ) *                                                 &
-                             ( sk(k,j,i+3) + sk(k,j,i-2) )                    &
-                                     )
-
-                diss_r    = -ABS( u_comp ) * (                                &
-                          ( 10.0_wp * ibit2 * adv_sca_5                       &
-                       +     3.0_wp * ibit1 * adv_sca_3                       &
-                       +              ibit0 * adv_sca_1                       &
-                          ) *                                                 &
-                             ( sk(k,j,i+1) - sk(k,j,i)   )                    &
-                   -      (  5.0_wp * ibit2 * adv_sca_5                       &
-                       +              ibit1 * adv_sca_3                       &
-                          ) *                                                 &
-                             ( sk(k,j,i+2) - sk(k,j,i-1) )                    &
-                   +      (           ibit2 * adv_sca_5                       &
-                          ) *                                                 &
-                             ( sk(k,j,i+3) - sk(k,j,i-2) )                    &
-                                             )
-
-                ibit5 = IBITS(wall_flags_0(k,j-1,i),5,1)
-                ibit4 = IBITS(wall_flags_0(k,j-1,i),4,1)
-                ibit3 = IBITS(wall_flags_0(k,j-1,i),3,1)
-
-                v_comp    = v(k,j,i) - v_gtrans
-                flux_s    = v_comp * (                                        &
-                          ( 37.0_wp * ibit5 * adv_sca_5                       &
-                       +     7.0_wp * ibit4 * adv_sca_3                       &
-                       +              ibit3 * adv_sca_1                       &
-                          ) *                                                 &
-                             ( sk(k,j,i)  + sk(k,j-1,i)     )                 &
-                    -     (  8.0_wp * ibit5 * adv_sca_5                       &
-                       +              ibit4 * adv_sca_3                       &
-                          ) *                                                 &
-                             ( sk(k,j+1,i) + sk(k,j-2,i)    )                 &
-                   +      (           ibit5 * adv_sca_5                       &
-                          ) *                                                 &
-                             ( sk(k,j+2,i) + sk(k,j-3,i)    )                 &
-                                     )
-
-                diss_s    = -ABS( v_comp ) * (                                &
-                          ( 10.0_wp * ibit5 * adv_sca_5                       &
-                       +     3.0_wp * ibit4 * adv_sca_3                       &
-                       +              ibit3 * adv_sca_1                       &
-                          ) *                                                 &
-                             ( sk(k,j,i)   - sk(k,j-1,i)  )                   &
-                   -      (  5.0_wp * ibit5 * adv_sca_5                       &
-                       +              ibit4 * adv_sca_3                       &
-                          ) *                                                 &
-                             ( sk(k,j+1,i) - sk(k,j-2,i)  )                   &
-                   +      (           ibit5 * adv_sca_5                       &
-                          ) *                                                 &
-                             ( sk(k,j+2,i) - sk(k,j-3,i)  )                   &
-                                             )
-
-                ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
-                ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
-                ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
-
-                v_comp    = v(k,j+1,i) - v_gtrans
-                flux_n    = v_comp * (                                        &
-                          ( 37.0_wp * ibit5 * adv_sca_5                       &
-                       +     7.0_wp * ibit4 * adv_sca_3                       &
-                       +              ibit3 * adv_sca_1                       &
-                          ) *                                                 &
-                             ( sk(k,j+1,i) + sk(k,j,i)   )                    &
-                   -      (  8.0_wp * ibit5 * adv_sca_5                       &
-                       +              ibit4 * adv_sca_3                       &
-                          ) *                                                 &
-                             ( sk(k,j+2,i) + sk(k,j-1,i) )                    &
-                   +      (           ibit5 * adv_sca_5                       &
-                          ) *                                                 &
-                             ( sk(k,j+3,i) + sk(k,j-2,i) )                    &
-                                     )
-
-                diss_n    = -ABS( v_comp ) * (                                &
-                          ( 10.0_wp * ibit5 * adv_sca_5                       &
-                       +     3.0_wp * ibit4 * adv_sca_3                       &
-                       +              ibit3 * adv_sca_1                       &
-                          ) *                                                 &
-                             ( sk(k,j+1,i) - sk(k,j,i)    )                   &
-                   -      (  5.0_wp * ibit5 * adv_sca_5                       &
-                       +              ibit4 * adv_sca_3                       &
-                          ) *                                                 &
-                             ( sk(k,j+2,i) - sk(k,j-1,i)  )                   &
-                   +      (           ibit5 * adv_sca_5                       &
-                          ) *                                                 &
-                             ( sk(k,j+3,i) - sk(k,j-2,i)  )                   &
-                                             )
-
-!
-!--             indizes k_m, k_mm, ... should be known at these point
-                ibit8 = IBITS(wall_flags_0(k-1,j,i),8,1)
-                ibit7 = IBITS(wall_flags_0(k-1,j,i),7,1)
-                ibit6 = IBITS(wall_flags_0(k-1,j,i),6,1)
-
-                k_pp  = k + 2 * ibit8
-                k_mm  = k - 2 * ( ibit7 + ibit8 )
-                k_mmm = k - 3 * ibit8
-
-                flux_d    = w(k-1,j,i) * (                                    &
-                           ( 37.0_wp * ibit8 * adv_sca_5                      &
-                        +     7.0_wp * ibit7 * adv_sca_3                      &
-                        +              ibit6 * adv_sca_1                      &
-                           ) *                                                &
-                                   ( sk(k,j,i)    + sk(k-1,j,i)  )            &
-                    -      (  8.0_wp * ibit8 * adv_sca_5                      &
-                          +            ibit7 * adv_sca_3                      &
-                           ) *                                                &
-                                   ( sk(k+1,j,i) + sk(k_mm,j,i)  )            &
-                    +      (           ibit8 * adv_sca_5                      &
-                           ) *     ( sk(k_pp,j,i)+ sk(k_mmm,j,i) )            &
-                                         )
-
-                diss_d    = -ABS( w(k-1,j,i) ) * (                            &
-                           ( 10.0_wp * ibit8 * adv_sca_5                      &
-                        +     3.0_wp * ibit7 * adv_sca_3                      &
-                        +              ibit6 * adv_sca_1                      &
-                           ) *                                                &
-                                   ( sk(k,j,i)    - sk(k-1,j,i)   )           &
-                    -      (  5.0_wp * ibit8 * adv_sca_5                      &
-                        +              ibit7 * adv_sca_3                      &
-                           ) *                                                &
-                                   ( sk(k+1,j,i)  - sk(k_mm,j,i)  )           &
-                    +      (           ibit8 * adv_sca_5                      &
-                           ) *                                                &
-                                   ( sk(k_pp,j,i) - sk(k_mmm,j,i) )           &
-                                                 )
-
-                ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
-                ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
-                ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
-
-                k_ppp = k + 3 * ibit8
-                k_pp  = k + 2 * ( 1 - ibit6  )
-                k_mm  = k - 2 * ibit8
-
-                flux_t    = 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    = -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) )           &
-                                         )
-!
-!--             Apply monotonic adjustment.
-                IF ( monotonic_adjustment )  THEN
-!
-!--                At first, calculate first order fluxes.
-                   u_comp = u(k,j,i) - u_gtrans
-                   fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )         &
-                         -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )         &
-                             ) * adv_sca_1 
-
-                   u_comp = u(k,j,i+1) - u_gtrans
-                   fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )         &
-                         -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )         &
-                             ) * adv_sca_1 
-
-                   v_comp = v(k,j,i) - v_gtrans
-                   fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )         &
-                         -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )         &
-                             ) * adv_sca_1 
-
-                   v_comp = v(k,j+1,i) - v_gtrans
-                   fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )         &
-                         -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )         &
-                             ) * adv_sca_1 
-
-                   fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
-                        -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )      &
-                            ) * adv_sca_1 * rho_air_zw(k)
-
-                   ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
-                        -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )        &
-                            ) * adv_sca_1 * rho_air_zw(k)
-!
-!--                Calculate ratio of upwind gradients. Note, Min/Max is just
-!--                to avoid if statements.
-                   rl     = ( MAX( 0.0_wp, u(k,j,i) - u_gtrans ) *            & 
-                               ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /&
-                                    ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )  &
-                                  ) +                                         & 
-                              MIN( 0.0_wp, u(k,j,i) - u_gtrans ) *            &
-                               ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /&
-                                    ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )  &
-                                  )                                           &
-                            ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
-
-                   rr     = ( MAX( 0.0_wp, u(k,j,i+1) - u_gtrans ) *          & 
-                               ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /&
-                                    ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )  &
-                                  ) +                                         & 
-                              MIN( 0.0_wp, u(k,j,i+1) - u_gtrans ) *          &
-                               ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /&
-                                    ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )  &
-                                  )                                           &
-                            ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
-
-                   rs     = ( MAX( 0.0_wp, v(k,j,i) - v_gtrans ) *            & 
-                               ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /&
-                                    ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )  &
-                                  ) +                                         & 
-                              MIN( 0.0_wp, v(k,j,i) - v_gtrans ) *            &
-                               ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /&
-                                    ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )  &
-                                  )                                           &
-                            ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
-
-                   rn     = ( MAX( 0.0_wp, v(k,j+1,i) - v_gtrans ) *          & 
-                               ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /&
-                                    ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )  &
-                                  ) +                                         & 
-                              MIN( 0.0_wp, v(k,j+1,i) - v_gtrans ) *          &
-                               ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /&
-                                    ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )  &
-                                  )                                           &
-                            ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
-!   
-!--                Reuse k_mm and compute k_mmm for the vertical gradient ratios. 
-!--                Note, for vertical advection below the third grid point above 
-!--                surface ( or below the model top) rd and rt are set to 0, i.e. 
-!--                use of first order scheme is enforced.
-                   k_mmm  = k - 3 * ibit8
-
-                   rd     = ( MAX( 0.0_wp, w(k-1,j,i) ) *                     & 
-                            ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) / &
-                                 ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )   &
-                               ) +                                            & 
-                              MIN( 0.0_wp, w(k-1,j,i) ) *                     &
-                            ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /     &
-                                 ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )  &
-                               )                                              &
-                            ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp ) 
- 
-                   rt     = ( MAX( 0.0_wp, w(k,j,i) ) *                       & 
-                            ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /   &
-                                 ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )     &
-                               ) +                                            & 
-                              MIN( 0.0_wp, w(k,j,i) ) *                       &
-                            ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /   &
-                                 ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )     &
-                               )                                              &
-                            ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
-!
-!--                Calculate empirical limiter function (van Albada2 limiter).
-                   phi_l = MIN( 1.0_wp, ( 2.0_wp * ABS( rl ) ) /              &
-                                        ( rl**2 + 1.0_wp ) ) 
-                   phi_r = MIN( 1.0_wp, ( 2.0_wp * ABS( rr ) ) /              &
-                                        ( rr**2 + 1.0_wp ) ) 
-                   phi_s = MIN( 1.0_wp, ( 2.0_wp * ABS( rs ) ) /              &
-                                        ( rs**2 + 1.0_wp ) ) 
-                   phi_n = MIN( 1.0_wp, ( 2.0_wp * ABS( rn ) ) /              &
-                                        ( rn**2 + 1.0_wp ) ) 
-                   phi_d = MIN( 1.0_wp, ( 2.0_wp * ABS( rd ) ) /              &
-                                        ( rd**2 + 1.0_wp ) ) 
-                   phi_t = MIN( 1.0_wp, ( 2.0_wp * ABS( rt ) ) /              &
-                                        ( rt**2 + 1.0_wp ) ) 
-!
-!--                Calculate the resulting monotone flux. 
-                   flux_l = fl_1 - phi_l * ( fl_1 - flux_l )
-                   flux_r = fr_1 - phi_r * ( fr_1 - flux_r )
-                   flux_s = fs_1 - phi_s * ( fs_1 - flux_s )
-                   flux_n = fn_1 - phi_n * ( fn_1 - flux_n )
-                   flux_d = fd_1 - phi_d * ( fd_1 - flux_d )
-                   flux_t = ft_1 - phi_t * ( ft_1 - flux_t )
-!
-!--                Moreover, modify dissipation flux according to the limiter. 
-                   diss_l = diss_l * phi_l
-                   diss_r = diss_r * phi_r
-                   diss_s = diss_s * phi_s
-                   diss_n = diss_n * phi_n
-                   diss_d = diss_d * phi_d
-                   diss_t = diss_t * phi_t
-
-                ENDIF
-!
-!--             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(wall_flags_0(k,j,i-1),0,1)    &
-                                         + IBITS(wall_flags_0(k,j,i-1),1,1)    &
-                                         + IBITS(wall_flags_0(k,j,i-1),2,1)    &
-                                         )                                     &
-                          ) * rho_air(k) * ddx                                 &
-                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
-                          - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
-                                         + IBITS(wall_flags_0(k,j-1,i),4,1)    &
-                                         + IBITS(wall_flags_0(k,j-1,i),5,1)    &
-                                         )                                     &
-                          ) * rho_air(k) * ddy                                 &
-                        + ( w(k,j,i) * rho_air_zw(k) *                         &
-                                         ( ibit6 + ibit7 + ibit8 )             &
-                          - w(k-1,j,i) * rho_air_zw(k-1) *                     &
-                                         ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
-                                         + IBITS(wall_flags_0(k-1,j,i),7,1)    &
-                                         + IBITS(wall_flags_0(k-1,j,i),8,1)    &
-                                         )                                     &      
-                          ) * ddzw(k)
-
-
-                tend(k,j,i) = - (                                             &
-                               ( flux_r + diss_r - flux_l - diss_l ) * ddx    &
-                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy    &
-                             + ( ( flux_t + diss_t ) -                        &
-                                 ( flux_d + diss_d )                          &
-                                                    ) * drho_air(k) * ddzw(k) &
-                                ) + div * sk(k,j,i)
-
-!++
-!--             Evaluation of statistics
-!                SELECT CASE ( sk_char )
-!
-!                   CASE ( 'pt' )
-!                      sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)         &
-!                       + ( flux_t + diss_t )                                &
-!                       *   weight_substep(intermediate_timestep_count)
-!                   CASE ( 'sa' )
-!                      sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)         &
-!                       + ( flux_t + diss_t )                                &
-!                       *   weight_substep(intermediate_timestep_count)
-!                   CASE ( 'q' )
-!                      sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)           &
-!                      + ( flux_t + diss_t )                                 &
-!                      *   weight_substep(intermediate_timestep_count)
-!                   CASE ( 'qr' )
-!                      sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn)         &
-!                      + ( flux_t + diss_t )                                 &
-!                      *   weight_substep(intermediate_timestep_count)
-!                   CASE ( 'nr' )
-!                      sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn)         &
-!                      + ( flux_t + diss_t )                                 &
-!                      *   weight_substep(intermediate_timestep_count)
-!
-!                END SELECT
-
-             ENDDO
-         ENDDO
-      ENDDO
-      !$acc end kernels
-
-    END SUBROUTINE advec_s_ws_acc
-
 
 !------------------------------------------------------------------------------!
@@ -5039,367 +4516,4 @@
 ! Description:
 ! ------------
-!> Advection of u - Call for all grid points - accelerator version
-!------------------------------------------------------------------------------!
-    SUBROUTINE advec_u_ws_acc
-
-       USE arrays_3d,                                                          &
-           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
-
-       USE constants,                                                          &
-           ONLY:  adv_mom_1, adv_mom_3, adv_mom_5
-
-       USE control_parameters,                                                 &
-           ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
-
-       USE grid_variables,                                                     &
-           ONLY:  ddx, ddy
-
-       USE indices,                                                            &
-           ONLY:  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb,  &
-                  nzb_max, nzt, wall_flags_0
-           
-       USE kinds
-       
-!        USE statistics,                                                       &
-!            ONLY:  hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i      !<
-       INTEGER(iwp) ::  ibit9  !<
-       INTEGER(iwp) ::  ibit10 !<
-       INTEGER(iwp) ::  ibit11 !<
-       INTEGER(iwp) ::  ibit12 !<
-       INTEGER(iwp) ::  ibit13 !<
-       INTEGER(iwp) ::  ibit14 !<
-       INTEGER(iwp) ::  ibit15 !<
-       INTEGER(iwp) ::  ibit16 !<
-       INTEGER(iwp) ::  ibit17 !<
-       INTEGER(iwp) ::  j      !<
-       INTEGER(iwp) ::  k      !<
-       INTEGER(iwp) ::  k_mmm  !<
-       INTEGER(iwp) ::  k_mm   !<
-       INTEGER(iwp) ::  k_pp   !<
-       INTEGER(iwp) ::  k_ppp  !<
-       INTEGER(iwp) ::  tn = 0 !<
-
-       REAL(wp)    ::  diss_d   !<
-       REAL(wp)    ::  diss_l   !<
-       REAL(wp)    ::  diss_n   !<
-       REAL(wp)    ::  diss_r   !<
-       REAL(wp)    ::  diss_s   !<
-       REAL(wp)    ::  diss_t   !<
-       REAL(wp)    ::  div      !<
-       REAL(wp)    ::  flux_d   !<
-       REAL(wp)    ::  flux_l   !<
-       REAL(wp)    ::  flux_n   !<
-       REAL(wp)    ::  flux_r   !<
-       REAL(wp)    ::  flux_s   !<
-       REAL(wp)    ::  flux_t   !<
-       REAL(wp)    ::  gu       !<
-       REAL(wp)    ::  gv       !<
-       REAL(wp)    ::  u_comp   !<
-       REAL(wp)    ::  u_comp_l !<
-       REAL(wp)    ::  v_comp   !<
-       REAL(wp)    ::  v_comp_s !<
-       REAL(wp)    ::  w_comp   !<
-
-
-       gu = 2.0_wp * u_gtrans
-       gv = 2.0_wp * v_gtrans
-
-!
-!--    Computation of fluxes and tendency terms
-       !$acc  kernels present( ddzw, tend, u, v, w, wall_flags_0 )
-       DO i = i_left, i_right
-          DO  j = j_south, j_north
-             DO  k = nzb+1, nzt
-
-                ibit11 = IBITS(wall_flags_0(k,j,i-1),11,1)
-                ibit10 = IBITS(wall_flags_0(k,j,i-1),10,1)
-                ibit9  = IBITS(wall_flags_0(k,j,i-1),9,1)
-
-                u_comp_l           = u(k,j,i) + u(k,j,i-1) - gu
-                flux_l             = u_comp_l * (                          &
-                                    ( 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             &
-                                    ) *                                     &
-                                  ( u(k,j,i+1) + u(k,j,i-2) )               &
-                             +      (           ibit11 * adv_mom_5             &
-                                    ) *                                     &
-                                  ( u(k,j,i+2) + u(k,j,i-3) )               &
-                                                )
-
-                diss_l             = - ABS( u_comp_l ) * (                &
-                                   ( 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             &
-                                   ) *                                     &
-                                 ( u(k,j,i+1) - u(k,j,i-2) )               &
-                            +      (           ibit11 * adv_mom_5             &
-                                   ) *                                     &
-                                 ( u(k,j,i+2) - u(k,j,i-3) )               &
-                                                         )
-
-                ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
-                ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
-                ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
-
-                u_comp    = u(k,j,i+1) + u(k,j,i)
-                flux_r    = ( u_comp   - gu ) * (                           &
-                          ( 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                        &
-                          ) *                                                &
-                                 ( u(k,j,i+2) + u(k,j,i-1) )                 &
-                   +      (           ibit11 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( u(k,j,i+3) + u(k,j,i-2) )                 &
-                                                 )
-
-                diss_r    = - ABS( u_comp    - gu ) * (                      &
-                          ( 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                        &
-                          ) *                                                &
-                                 ( u(k,j,i+2) - u(k,j,i-1) )                 &
-                   +      (           ibit11 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( u(k,j,i+3) - u(k,j,i-2) )                 &
-                                                     )
-
-                ibit14 = IBITS(wall_flags_0(k,j-1,i),14,1)
-                ibit13 = IBITS(wall_flags_0(k,j-1,i),13,1)
-                ibit12 = IBITS(wall_flags_0(k,j-1,i),12,1)
-
-                v_comp_s                 = v(k,j,i) + v(k,j,i-1) - gv
-                flux_s                   = v_comp_s * (                       &
-                                   ( 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                    &
-                                   ) *                                        &
-                                     ( u(k,j+1,i) + u(k,j-2,i) )              &
-                        +      (               ibit14 * adv_mom_5                    &
-                               ) *                                            &
-                                     ( u(k,j+2,i) + u(k,j-3,i) )              &
-                                               )
-
-                diss_s                  = - ABS ( v_comp_s ) * (              &
-                                   ( 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                &
-                                   ) *                                        &
-                                     ( u(k,j+1,i) - u(k,j-2,i) )              &
-                            +      (           ibit14 * adv_mom_5                &
-                                   ) *                                        &
-                                     ( u(k,j+2,i) - u(k,j-3,i) )              &
-                                                         )
-
-
-                ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
-                ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
-                ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
-
-                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
-                flux_n    = v_comp * (                                       &
-                          ( 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                        &
-                          ) *                                                &
-                                 ( u(k,j+2,i) + u(k,j-1,i) )                 &
-                   +      (           ibit14 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( u(k,j+3,i) + u(k,j-2,i) )                 &
-                                                 )
-
-                diss_n    = - ABS ( v_comp ) * (                             &
-                          ( 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                        &
-                          ) *                                                &
-                                 ( u(k,j+2,i) - u(k,j-1,i) )                 &
-                   +      (           ibit14 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( u(k,j+3,i) - u(k,j-2,i) )                 &
-                                                      )
-
-                ibit17 = IBITS(wall_flags_0(k-1,j,i),17,1)
-                ibit16 = IBITS(wall_flags_0(k-1,j,i),16,1)
-                ibit15 = IBITS(wall_flags_0(k-1,j,i),15,1)
-
-                k_pp  = k + 2 * ibit17
-                k_mm  = k - 2 * ( ibit16 + ibit17 )
-                k_mmm = k - 3 * ibit17
-
-                w_comp    = w(k-1,j,i) + w(k-1,j,i-1)
-                flux_d    = w_comp  * (                                      &
-                          ( 37.0_wp * ibit17 * adv_mom_5                        &
-                       +     7.0_wp * ibit16 * adv_mom_3                        &
-                       +              ibit15 * adv_mom_1                        &
-                          ) *                                                &
-                             ( u(k,j,i)    + u(k-1,j,i)   )                  &
-                   -      (  8.0_wp * ibit17 * adv_mom_5                        &
-                       +              ibit16 * adv_mom_3                        &
-                          ) *                                                &
-                             ( u(k+1,j,i) + u(k_mm,j,i)   )                  &
-                   +      (           ibit17 * adv_mom_5                        &
-                          ) *                                                 &
-                             ( u(k_pp,j,i) + u(k_mmm,j,i) )                  &
-                                      )
-
-                diss_d    = - ABS( w_comp ) * (                              &
-                          ( 10.0_wp * ibit17 * adv_mom_5                        &
-                       +     3.0_wp * ibit16 * adv_mom_3                        &
-                       +              ibit15 * adv_mom_1                        &
-                          ) *                                                &
-                             ( u(k,j,i)     - u(k-1,j,i)  )                  &
-                   -      (  5.0_wp * ibit17 * adv_mom_5                        &
-                       +              ibit16 * adv_mom_3                        &
-                          ) *                                                &
-                             ( u(k+1,j,i)  - u(k_mm,j,i)  )                  &
-                   +      (           ibit17 * adv_mom_5                        &
-                           ) *                                               &
-                             ( u(k_pp,j,i) - u(k_mmm,j,i) )                  &
-                                              )
-!
-!--             k index has to be modified near bottom and top, else array
-!--             subscripts will be exceeded.
-                ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
-                ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
-                ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
-
-                k_ppp = k + 3 * ibit17
-                k_pp  = k + 2 * ( 1 - ibit15  )
-                k_mm  = k - 2 * ibit17
-
-                w_comp    = w(k,j,i) + w(k,j,i-1)
-                flux_t    = w_comp * 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    = - 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                        &
-                          ) *                                                &
-                             ( 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) )                  &
-                                              )
-!
-!--             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_comp    * ( ibit9 + ibit10 + ibit11 )             &
-                - ( u(k,j,i)   + u(k,j,i-1)   )                               &
-                                    * ( IBITS(wall_flags_0(k,j,i-1),9,1)      &
-                                      + IBITS(wall_flags_0(k,j,i-1),10,1)     &
-                                      + IBITS(wall_flags_0(k,j,i-1),11,1)     &
-                                      )                                       &
-                  ) * rho_air(k) * ddx                                        &
-               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
-                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
-                                    * ( IBITS(wall_flags_0(k,j-1,i),12,1)     &
-                                      + IBITS(wall_flags_0(k,j-1,i),13,1)     &
-                                      + IBITS(wall_flags_0(k,j-1,i),14,1)     &
-                                      )                                       &
-                  ) * rho_air(k) * 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(wall_flags_0(k-1,j,i),15,1)     &
-                                      + IBITS(wall_flags_0(k-1,j,i),16,1)     &
-                                      + IBITS(wall_flags_0(k-1,j,i),17,1)     &
-                                      )                                       &
-                  ) * ddzw(k)   &
-                ) * 0.5_wp
-
-
-                tend(k,j,i) = - (                                              &
-                               ( flux_r + diss_r - flux_l - diss_l ) * ddx     &
-                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy     &
-                             + ( ( flux_t + diss_t ) -                         &
-                                 ( flux_d + diss_d )                           &
-                                                     ) * drho_air(k) * ddzw(k) &
-                                ) + div * u(k,j,i)
-
-!++
-!--             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    *                                 &
-!                                ( u_comp    - 2.0_wp * hom(k,1,1,0) )            &
-!                              / ( u_comp    - gu + 1.0E-20_wp   )             &
-!                              +   diss_r    *                                 &
-!                                  ABS( u_comp    - 2.0_wp * hom(k,1,1,0) )       &
-!                              / ( ABS( u_comp    - 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    + diss_t    )                     &
-!                              *   weight_substep(intermediate_timestep_count)
-             ENDDO
-          ENDDO
-       ENDDO
-       !$acc end kernels
-
-!++
-!       sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
-
-    END SUBROUTINE advec_u_ws_acc
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
 !> Advection of v - Call for all grid points
 !------------------------------------------------------------------------------!
@@ -5906,370 +5020,4 @@
     
     
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Advection of v - Call for all grid points - accelerator version
-!------------------------------------------------------------------------------!
-    SUBROUTINE advec_v_ws_acc
-
-       USE arrays_3d,                                                          &
-           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
-
-       USE constants,                                                          &
-           ONLY:  adv_mom_1, adv_mom_3, adv_mom_5
-
-       USE control_parameters,                                                 &
-           ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
-
-       USE grid_variables,                                                     &
-           ONLY:  ddx, ddy
-
-       USE indices,                                                            &
-           ONLY:  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb,  &
-                  nzb_max, nzt, wall_flags_0
-           
-       USE kinds
-       
-!        USE statistics,                                                       &
-!            ONLY:  hom, sums_vs2_ws_l, sums_wsvs_ws_l, weight_substep
-
-       IMPLICIT NONE
-
-
-       INTEGER(iwp) ::  i      !<
-       INTEGER(iwp) ::  ibit18 !<
-       INTEGER(iwp) ::  ibit19 !<
-       INTEGER(iwp) ::  ibit20 !<
-       INTEGER(iwp) ::  ibit21 !<
-       INTEGER(iwp) ::  ibit22 !<
-       INTEGER(iwp) ::  ibit23 !<
-       INTEGER(iwp) ::  ibit24 !<
-       INTEGER(iwp) ::  ibit25 !<
-       INTEGER(iwp) ::  ibit26 !<
-       INTEGER(iwp) ::  j      !<
-       INTEGER(iwp) ::  k      !<
-       INTEGER(iwp) ::  k_mm   !<
-       INTEGER(iwp) ::  k_mmm  !<
-       INTEGER(iwp) ::  k_pp   !<
-       INTEGER(iwp) ::  k_ppp  !<
-       INTEGER(iwp) ::  tn = 0 !<
-
-       REAL(wp)    ::  diss_d   !<
-       REAL(wp)    ::  diss_l   !<
-       REAL(wp)    ::  diss_n   !<
-       REAL(wp)    ::  diss_r   !<
-       REAL(wp)    ::  diss_s   !<
-       REAL(wp)    ::  diss_t   !<
-       REAL(wp)    ::  div      !<
-       REAL(wp)    ::  flux_d   !<
-       REAL(wp)    ::  flux_l   !<
-       REAL(wp)    ::  flux_n   !<
-       REAL(wp)    ::  flux_r   !<
-       REAL(wp)    ::  flux_s   !<
-       REAL(wp)    ::  flux_t   !<
-       REAL(wp)    ::  gu       !<
-       REAL(wp)    ::  gv       !<
-       REAL(wp)    ::  u_comp   !<
-       REAL(wp)    ::  u_comp_l !<
-       REAL(wp)    ::  v_comp   !<
-       REAL(wp)    ::  v_comp_s !<
-       REAL(wp)    ::  w_comp   !<
-
-       gu = 2.0_wp * u_gtrans
-       gv = 2.0_wp * v_gtrans
-
-!
-!--    Computation of fluxes and tendency terms
-       !$acc kernels present( ddzw, tend, u, v, w, wall_flags_0 )
-       DO  i = i_left, i_right
-          DO  j = j_south, j_north
-             DO  k = nzb+1, nzt
-
-                ibit20 = IBITS(wall_flags_0(k,j,i-1),20,1)
-                ibit19 = IBITS(wall_flags_0(k,j,i-1),19,1)
-                ibit18 = IBITS(wall_flags_0(k,j,i-1),18,1)
-
-                u_comp_l                 = u(k,j-1,i) + u(k,j,i) - gu
-                flux_l                   = u_comp_l * (                          &
-                                      ( 37.0_wp * ibit20 * adv_mom_5              &
-                                   +     7.0_wp * ibit19 * adv_mom_3              &
-                                   +              ibit18 * adv_mom_1              &
-                                      ) *                                      &
-                                     ( v(k,j,i)   + v(k,j,i-1) )               &
-                               -      (  8.0_wp * ibit20 * adv_mom_5              &
-                                   +              ibit19 * adv_mom_3              &
-                                      ) *                                      &
-                                     ( v(k,j,i+1) + v(k,j,i-2) )               &
-                               +      (           ibit20 * adv_mom_5              &
-                                      ) *                                      &
-                                     ( v(k,j,i+2) + v(k,j,i-3) )               &
-                                                 )
-
-                diss_l                   = - ABS( u_comp_l ) * (                 &
-                                      ( 10.0_wp * ibit20 * adv_mom_5              &
-                                   +     3.0_wp * ibit19 * adv_mom_3              &
-                                   +              ibit18 * adv_mom_1              &
-                                      ) *                                      &
-                                     ( v(k,j,i)   - v(k,j,i-1) )               &
-                               -      (  5.0_wp * ibit20 * adv_mom_5              &
-                                   +              ibit19 * adv_mom_3              &
-                                      ) *                                      &
-                                     ( v(k,j,i+1) - v(k,j,i-2) )               &
-                               +      (           ibit20 * adv_mom_5              &
-                                      ) *                                      &
-                                     ( v(k,j,i+2) - v(k,j,i-3) )               &
-                                                           )
-
-                ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
-                ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
-                ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
-
-                u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
-                flux_r    = u_comp * (                                       &
-                          ( 37.0_wp * ibit20 * adv_mom_5                        &
-                       +     7.0_wp * ibit19 * adv_mom_3                        &
-                       +              ibit18 * adv_mom_1                        &
-                          ) *                                                &
-                                 ( v(k,j,i+1) + v(k,j,i)   )                 &
-                   -      (  8.0_wp * ibit20 * adv_mom_5                        &
-                       +              ibit19 * adv_mom_3                        &
-                          ) *                                                &
-                                 ( v(k,j,i+2) + v(k,j,i-1) )                 &
-                   +      (           ibit20 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( v(k,j,i+3) + v(k,j,i-2) )                 &
-                                     )
-
-                diss_r    = - ABS( u_comp ) * (                              &
-                          ( 10.0_wp * ibit20 * adv_mom_5                        &
-                       +     3.0_wp * ibit19 * adv_mom_3                        &
-                       +              ibit18 * adv_mom_1                        &
-                          ) *                                                &
-                                 ( v(k,j,i+1) - v(k,j,i)  )                  &
-                   -      (  5.0_wp * ibit20 * adv_mom_5                        &
-                       +              ibit19 * adv_mom_3                        &
-                          ) *                                                &
-                                 ( v(k,j,i+2) - v(k,j,i-1) )                 &
-                   +      (           ibit20 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( v(k,j,i+3) - v(k,j,i-2) )                 &
-                                              )
-
-                ibit23 = IBITS(wall_flags_0(k,j-1,i),23,1)
-                ibit22 = IBITS(wall_flags_0(k,j-1,i),22,1)
-                ibit21 = IBITS(wall_flags_0(k,j-1,i),21,1)
-
-
-                v_comp_s              = v(k,j,i) + v(k,j-1,i) - gv
-                flux_s                = v_comp_s    * (                       &
-                                   ( 37.0_wp * ibit23 * adv_mom_5                &
-                                +     7.0_wp * ibit22 * adv_mom_3                &
-                                +              ibit21 * adv_mom_1                &
-                                   ) *                                        &
-                                     ( v(k,j,i)   + v(k,j-1,i) )              &
-                            -      (  8.0_wp * ibit23 * adv_mom_5                &
-                                +              ibit22 * adv_mom_3                &
-                                   ) *                                        &
-                                     ( v(k,j+1,i) + v(k,j-2,i) )              &
-                            +      (           ibit23 * adv_mom_5                &
-                                   ) *                                        &
-                                     ( v(k,j+2,i) + v(k,j-3,i) )              &
-                                                 )
-
-                diss_s                = - ABS( v_comp_s ) * (                 &
-                                   ( 10.0_wp * ibit23 * adv_mom_5                &
-                                +     3.0_wp * ibit22 * adv_mom_3                &
-                                +              ibit21 * adv_mom_1                &
-                                   ) *                                        &
-                                     ( v(k,j,i)   - v(k,j-1,i) )              &
-                            -      (  5.0_wp * ibit23 * adv_mom_5                &
-                                +              ibit22 * adv_mom_3                &
-                                   ) *                                        &
-                                     ( v(k,j+1,i) - v(k,j-2,i) )              &
-                            +      (           ibit23 * adv_mom_5                &
-                                   ) *                                        &
-                                     ( v(k,j+2,i) - v(k,j-3,i) )              &
-                                                          )
-
-                ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
-                ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
-                ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
-
-                v_comp = v(k,j+1,i) + v(k,j,i)
-                flux_n = ( v_comp - gv ) * (                                 &
-                          ( 37.0_wp * ibit23 * adv_mom_5                        &
-                       +     7.0_wp * ibit22 * adv_mom_3                        &
-                       +              ibit21 * adv_mom_1                        &
-                          ) *                                                &
-                                 ( v(k,j+1,i) + v(k,j,i)   )                 &
-                   -      (  8.0_wp * ibit23 * adv_mom_5                        &
-                       +              ibit22 * adv_mom_3                        &
-                          ) *                                                &
-                                 ( v(k,j+2,i) + v(k,j-1,i) )                 &
-                   +      (           ibit23 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( v(k,j+3,i) + v(k,j-2,i) )                 &
-                                     )
-
-                diss_n = - ABS( v_comp - gv ) * (                         &
-                          ( 10.0_wp * ibit23 * adv_mom_5                        &
-                       +     3.0_wp * ibit22 * adv_mom_3                        &
-                       +              ibit21 * adv_mom_1                        &
-                          ) *                                                &
-                                 ( v(k,j+1,i) - v(k,j,i)  )                  &
-                   -      (  5.0_wp * ibit23 * adv_mom_5                        &
-                       +              ibit22 * adv_mom_3                        &
-                          ) *                                                &
-                                 ( v(k,j+2,i) - v(k,j-1,i) )                 &
-                   +      (           ibit23 * adv_mom_5                        &
-                          ) *                                                &
-                                 ( v(k,j+3,i) - v(k,j-2,i) )                 &
-                                                     )
-
-                ibit26 = IBITS(wall_flags_0(k-1,j,i),26,1)
-                ibit25 = IBITS(wall_flags_0(k-1,j,i),25,1)
-                ibit24 = IBITS(wall_flags_0(k-1,j,i),24,1)
-
-                k_pp  = k + 2 * ibit26
-                k_mm  = k - 2 * ( ibit25 + ibit26 )
-                k_mmm = k - 3 * ibit26
-
-                w_comp    = w(k-1,j-1,i) + w(k-1,j,i)
-                flux_d    = w_comp  * (                                      &
-                          ( 37.0_wp * ibit26 * adv_mom_5                        &
-                       +     7.0_wp * ibit25 * adv_mom_3                        &
-                       +              ibit24 * adv_mom_1                        &
-                          ) *                                                &
-                             ( v(k,j,i)     + v(k-1,j,i)  )                  &
-                   -      (  8.0_wp * ibit26 * adv_mom_5                        &
-                       +              ibit25 * adv_mom_3                        &
-                          ) *                                                &
-                             ( v(k+1,j,i)  + v(k_mm,j,i)  )                  &
-                   +      (           ibit26 * adv_mom_5                        &
-                          ) *                                                &
-                             ( v(k_pp,j,i) + v(k_mmm,j,i) )                  &
-                                      )
-
-                diss_d    = - ABS( w_comp ) * (                              &
-                          ( 10.0_wp * ibit26 * adv_mom_5                        &
-                       +     3.0_wp * ibit25 * adv_mom_3                        &
-                       +              ibit24 * adv_mom_1                        &
-                          ) *                                                &
-                             ( v(k,j,i)     - v(k-1,j,i)  )                  &
-                   -      (  5.0_wp * ibit26 * adv_mom_5                        &
-                       +              ibit25 * adv_mom_3                        &
-                          ) *                                                &
-                             ( v(k+1,j,i)  - v(k_mm,j,i)  )                  &
-                   +      (           ibit26 * adv_mom_5                        &
-                          ) *                                                &
-                             ( v(k_pp,j,i) - v(k_mmm,j,i) )                  &
-                                               )
-!
-!--             k index has to be modified near bottom and top, else array
-!--             subscripts will be exceeded.
-                ibit26 = IBITS(wall_flags_0(k,j,i),26,1)
-                ibit25 = IBITS(wall_flags_0(k,j,i),25,1)
-                ibit24 = IBITS(wall_flags_0(k,j,i),24,1)
-
-                k_ppp = k + 3 * ibit26
-                k_pp  = k + 2 * ( 1 - ibit24  )
-                k_mm  = k - 2 * ibit26
-
-                w_comp    = w(k,j-1,i) + w(k,j,i)
-                flux_t    = w_comp * 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    = - ABS( w_comp ) * 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) )                  &
-                                               )
-!
-!--             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_comp     + gu )                                 &
-                                       * ( ibit18 + ibit19 + ibit20 )         &
-                - ( u(k,j-1,i)   + u(k,j,i) )                                 &
-                                       * ( IBITS(wall_flags_0(k,j,i-1),18,1)  &
-                                         + IBITS(wall_flags_0(k,j,i-1),19,1)  &
-                                         + IBITS(wall_flags_0(k,j,i-1),20,1)  &
-                                         )                                    &
-                  ) * rho_air(k) * ddx                                        &
-               +  ( v_comp                                                    &
-                                       * ( ibit21 + ibit22 + ibit23 )         &
-                - ( v(k,j,i)     + v(k,j-1,i) )                               &
-                                       * ( IBITS(wall_flags_0(k,j-1,i),21,1)  &
-                                         + IBITS(wall_flags_0(k,j-1,i),22,1)  &
-                                         + IBITS(wall_flags_0(k,j-1,i),23,1)  &
-                                         )                                    &
-                  ) * rho_air(k) * ddy                                        &
-               +  ( w_comp * rho_air_zw(k)                                    &
-                                       * ( ibit24 + ibit25 + ibit26 )         &
-                - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
-                                       * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
-                                         + IBITS(wall_flags_0(k-1,j,i),25,1)  &
-                                         + IBITS(wall_flags_0(k-1,j,i),26,1)  &
-                                         )                                    &
-                   ) * ddzw(k)   &
-                ) * 0.5_wp
-
-
-                tend(k,j,i) = - (                                              &
-                               ( flux_r + diss_r - flux_l - diss_l ) * ddx     &
-                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy     &
-                             + ( ( flux_t + diss_t ) -                         &
-                                 ( flux_d + diss_d )                           &
-                               ) * drho_air(k) * ddzw(k)                       &
-                                ) + div * v(k,j,i)
-
-
-!++
-!--             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                                           &
-!                      * ( v_comp - 2.0_wp * hom(k,1,2,0) )                    &
-!                      / ( v_comp - gv + 1.0E-20_wp )                       &
-!                      +   diss_n                                           &
-!                      *   ABS( v_comp - 2.0_wp * hom(k,1,2,0) )               &
-!                      / ( ABS( v_comp - gv ) +1.0E-20_wp ) )               &
-!                      *   weight_substep(intermediate_timestep_count)
-!
-!--              Statistical Evaluation of w'v'.
-!                 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                &
-!                              + ( flux_t + diss_t )                         &
-!                              *   weight_substep(intermediate_timestep_count)
-
-             ENDDO
-          ENDDO
-       ENDDO
-       !$acc end kernels
-
-!++
-!       sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
-
-    END SUBROUTINE advec_v_ws_acc
     
     
@@ -6756,352 +5504,3 @@
 
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Advection of w - Call for all grid points - accelerator version
-!------------------------------------------------------------------------------!
-    SUBROUTINE advec_w_ws_acc
-
-       USE arrays_3d,                                                          &
-           ONLY:  ddzu, drho_air_zw, tend, u, v, w, rho_air, rho_air_zw
-
-       USE constants,                                                          &
-           ONLY:  adv_mom_1, adv_mom_3, adv_mom_5
-
-       USE control_parameters,                                                 &
-           ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
-
-       USE grid_variables,                                                     &
-           ONLY:  ddx, ddy
-
-       USE indices,                                                            &
-           ONLY:  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb,  &
-                  nzb_max, nzt, wall_flags_0, wall_flags_00
-           
-       USE kinds
-       
-!        USE statistics,                                                       &
-!            ONLY:  hom, sums_ws2_ws_l, weight_substep
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i      !<
-       INTEGER(iwp) ::  ibit27 !<
-       INTEGER(iwp) ::  ibit28 !<
-       INTEGER(iwp) ::  ibit29 !<
-       INTEGER(iwp) ::  ibit30 !<
-       INTEGER(iwp) ::  ibit31 !<
-       INTEGER(iwp) ::  ibit32 !<
-       INTEGER(iwp) ::  ibit33 !<
-       INTEGER(iwp) ::  ibit34 !<
-       INTEGER(iwp) ::  ibit35 !<
-       INTEGER(iwp) ::  j      !<
-       INTEGER(iwp) ::  k      !<
-       INTEGER(iwp) ::  k_mmm  !<
-       INTEGER(iwp) ::  k_mm   !<
-       INTEGER(iwp) ::  k_pp   !<
-       INTEGER(iwp) ::  k_ppp  !<
-       INTEGER(iwp) ::  tn = 0 !<
-
-       REAL(wp)    ::  diss_d   !<
-       REAL(wp)    ::  diss_l   !<
-       REAL(wp)    ::  diss_n   !<
-       REAL(wp)    ::  diss_r   !<
-       REAL(wp)    ::  diss_s   !<
-       REAL(wp)    ::  diss_t   !<
-       REAL(wp)    ::  div      !<
-       REAL(wp)    ::  flux_d   !<
-       REAL(wp)    ::  flux_l   !<
-       REAL(wp)    ::  flux_n   !<
-       REAL(wp)    ::  flux_r   !<
-       REAL(wp)    ::  flux_s   !<
-       REAL(wp)    ::  flux_t   !<
-       REAL(wp)    ::  gu       !<
-       REAL(wp)    ::  gv       !<
-       REAL(wp)    ::  u_comp   !<
-       REAL(wp)    ::  u_comp_l !<
-       REAL(wp)    ::  v_comp   !<
-       REAL(wp)    ::  v_comp_s !<
-       REAL(wp)    ::  w_comp   !<
-
-       gu = 2.0_wp * u_gtrans
-       gv = 2.0_wp * v_gtrans
-
-
-!
-!--    Computation of fluxes and tendency terms
-       !$acc kernels present( ddzu, tend, u, v, w, wall_flags_0, wall_flags_00 )
-       DO i = i_left, i_right
-          DO  j = j_south, j_north
-             DO  k = nzb+1, nzt
-
-                ibit27 = IBITS(wall_flags_0(k,j,i-1),27,1)
-                ibit28 = IBITS(wall_flags_0(k,j,i-1),28,1)
-                ibit29 = IBITS(wall_flags_0(k,j,i-1),29,1)
-
-                u_comp_l                 = u(k+1,j,i) + u(k,j,i) - gu
-                flux_l                   = u_comp_l * (                        &
-                                      ( 37.0_wp * ibit29 * adv_mom_5              &
-                                   +     7.0_wp * ibit28 * adv_mom_3              &
-                                   +              ibit27 * adv_mom_1              &
-                                      ) *                                      &
-                                     ( w(k,j,i)   + w(k,j,i-1) )               &
-                               -      (  8.0_wp * ibit29 * adv_mom_5              &
-                                   +              ibit28 * adv_mom_3              &
-                                      ) *                                      &
-                                     ( w(k,j,i+1) + w(k,j,i-2) )               &
-                               +      (           ibit29 * adv_mom_5              &
-                                      ) *                                      &
-                                     ( w(k,j,i+2) + w(k,j,i-3) )               &
-                                                 )
-
-                diss_l                    = - ABS( u_comp_l ) * (              &
-                                        ( 10.0_wp * ibit29 * adv_mom_5            &
-                                     +     3.0_wp * ibit28 * adv_mom_3            &
-                                     +              ibit27 * adv_mom_1            &
-                                        ) *                                    &
-                                     ( w(k,j,i)   - w(k,j,i-1) )               &
-                                 -      (  5.0_wp * ibit29 * adv_mom_5            &
-                                     +              ibit28 * adv_mom_3            &
-                                        ) *                                    &
-                                     ( w(k,j,i+1) - w(k,j,i-2) )               &
-                                 +      (           ibit29 * adv_mom_5            &
-                                        ) *                                    &
-                                     ( w(k,j,i+2) - w(k,j,i-3) )               &
-                                                            )
-
-                ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
-                ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
-                ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
-
-                u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
-                flux_r    = u_comp * (                                       &
-                          ( 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    = - ABS( u_comp ) * (                              &
-                          ( 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(wall_flags_00(k,j-1,i),0,1)
-                ibit31 = IBITS(wall_flags_0(k,j-1,i),31,1)
-                ibit30 = IBITS(wall_flags_0(k,j-1,i),30,1)
-
-                v_comp_s               = v(k+1,j,i) + v(k,j,i) - gv
-                flux_s                 = v_comp_s * (                         &
-                                    ( 37.0_wp * ibit32 * adv_mom_5               &
-                                 +     7.0_wp * ibit31 * adv_mom_3               &
-                                 +              ibit30 * adv_mom_1               &
-                                    ) *                                       &
-                                     ( w(k,j,i)   + w(k,j-1,i) )              &
-                             -      (  8.0_wp * ibit32 * adv_mom_5               &
-                                 +              ibit31 * adv_mom_3               &
-                                    ) *                                       &
-                                     ( w(k,j+1,i) + w(k,j-2,i) )              &
-                             +      (           ibit32 * adv_mom_5               &
-                                    ) *                                       &
-                                     ( w(k,j+2,i) + w(k,j-3,i) )              &
-                                               )
-
-                diss_s                 = - ABS( v_comp_s ) * (                &
-                                    ( 10.0_wp * ibit32 * adv_mom_5               &
-                                 +     3.0_wp * ibit31 * adv_mom_3               &
-                                 +              ibit30 * adv_mom_1               &
-                                    ) *                                       &
-                                     ( w(k,j,i)   - w(k,j-1,i) )              &
-                             -      (  5.0_wp * ibit32 * adv_mom_5               &
-                                 +              ibit31 * adv_mom_3               &
-                                    ) *                                       &
-                                     ( w(k,j+1,i) - w(k,j-2,i) )              &
-                             +      (           ibit32 * adv_mom_5               &
-                                    ) *                                       &
-                                     ( w(k,j+2,i) - w(k,j-3,i) )              &
-                                                        )
-
-                ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
-                ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
-                ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
-
-                v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
-                flux_n    = v_comp * (                                       &
-                          ( 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    = - ABS( v_comp ) * (                              &
-                          ( 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) )                 &
-                                              )
-
-                ibit35 = IBITS(wall_flags_00(k-1,j,i),3,1)
-                ibit34 = IBITS(wall_flags_00(k-1,j,i),2,1)
-                ibit33 = IBITS(wall_flags_00(k-1,j,i),1,1)
-
-                k_pp  = k + 2 * ibit35
-                k_mm  = k - 2 * ( ibit34 + ibit35 )
-                k_mmm = k - 3 * ibit35
-
-                w_comp    = w(k,j,i) + w(k-1,j,i)
-                flux_d    = w_comp  * (                                      &
-                          ( 37.0_wp * ibit35 * adv_mom_5                        &
-                       +     7.0_wp * ibit34 * adv_mom_3                        &
-                       +              ibit33 * adv_mom_1                        &
-                          ) *                                                &
-                             ( w(k,j,i)    + w(k-1,j,i)   )                  &
-                   -      (  8.0_wp * ibit35 * adv_mom_5                        &
-                       +              ibit34 * adv_mom_3                        &
-                          ) *                                                &
-                             ( w(k+1,j,i)  + w(k_mm,j,i)  )                  &
-                   +      (           ibit35 * adv_mom_5                        &
-                          ) *                                                &
-                             ( w(k_pp,j,i) + w(k_mmm,j,i) )                  &
-                                       )
-
-                diss_d    = - ABS( w_comp ) * (                              &
-                          ( 10.0_wp * ibit35 * adv_mom_5                        &
-                       +     3.0_wp * ibit34 * adv_mom_3                        &
-                       +              ibit33 * adv_mom_1                        &
-                          ) *                                                &
-                             ( w(k,j,i)    - w(k-1,j,i)   )                  &
-                   -      (  5.0_wp * ibit35 * adv_mom_5                        &
-                       +              ibit34 * adv_mom_3                        &
-                          ) *                                                &
-                             ( w(k+1,j,i)  - w(k_mm,j,i)  )                  &
-                   +      (           ibit35 * adv_mom_5                        &
-                          ) *                                                &
-                             ( w(k_pp,j,i) - w(k_mmm,j,i) )                  &
-                                               )
-
-!
-!--             k index has to be modified near bottom and top, else array
-!--             subscripts will be exceeded.
-                ibit35 = IBITS(wall_flags_00(k,j,i),3,1)
-                ibit34 = IBITS(wall_flags_00(k,j,i),2,1)
-                ibit33 = IBITS(wall_flags_00(k,j,i),1,1)
-
-                k_ppp = k + 3 * ibit35
-                k_pp  = k + 2 * ( 1 - ibit33  )
-                k_mm  = k - 2 * ibit35
-
-                w_comp    = w(k+1,j,i) + w(k,j,i)
-                flux_t    = w_comp * 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    = - ABS( w_comp ) * 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) )                  &
-                                               )
-!
-!--             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_comp + gu ) * ( ibit27 + ibit28 + ibit29 )      &
-                  - ( u(k+1,j,i) + u(k,j,i)   )                               & 
-                                    * ( IBITS(wall_flags_0(k,j,i-1),27,1)     &
-                                      + IBITS(wall_flags_0(k,j,i-1),28,1)     &
-                                      + IBITS(wall_flags_0(k,j,i-1),29,1)     &
-                                      )                                       &
-                  ) * rho_air_zw(k) * ddx                                     &
-              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
-                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
-                                    * ( IBITS(wall_flags_0(k,j-1,i),30,1)     &
-                                      + IBITS(wall_flags_0(k,j-1,i),31,1)     &
-                                      + IBITS(wall_flags_00(k,j-1,i),0,1)     &
-                                      )                                       &
-                  ) * rho_air_zw(k) * ddy                                     &
-              +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
-                - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
-                                    * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
-                                      + IBITS(wall_flags_00(k-1,j,i),2,1)     &
-                                      + IBITS(wall_flags_00(k-1,j,i),3,1)     &
-                                      )                                       & 
-                  ) * ddzu(k+1)   &
-                ) * 0.5_wp
-
-
-                tend(k,j,i) = - (                                                &
-                               ( flux_r + diss_r - flux_l - diss_l ) * ddx       &
-                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy       &
-                             + ( ( flux_t + diss_t ) -                           &
-                                 ( flux_d + diss_d ) * rho_air(k)                &
-                               ) * drho_air_zw(k) * ddzu(k+1)                    &
-                                 ) + div * w(k,j,i)
-
-
-!++
-!--             Statistical Evaluation of w'w'.
-!                sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                    &
-!                               + ( flux_t + diss_t )                    &
-!                               *   weight_substep(intermediate_timestep_count)
-
-             ENDDO
-          ENDDO
-       ENDDO
-       !$acc end kernels
-
-    END SUBROUTINE advec_w_ws_acc
-
  END MODULE advec_ws
