Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 1011)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! ------------------
-!
+! accelerator versions (*_acc) added
 !
 ! Former revisions:
@@ -88,5 +88,6 @@
 
     PRIVATE
-    PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, &
+    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_statistics
 
@@ -109,4 +110,8 @@
     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
@@ -114,8 +119,16 @@
     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
@@ -2335,4 +2348,286 @@
 
 !------------------------------------------------------------------------------!
+! Scalar advection - Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE advec_s_ws_acc ( sk, sk_char )
+
+       USE arrays_3d
+       USE constants
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE statistics
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN = *), INTENT(IN)    :: sk_char
+
+       INTEGER ::  i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6,        &
+                   ibit7, ibit8, j, k, k_mm, k_mmm, k_pp, k_ppp, tn = 0
+
+       REAL    :: diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div, flux_d, &
+                  flux_l, flux_n, flux_r, flux_s, flux_t, u_comp, v_comp
+
+       REAL, 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 )
+       !$acc loop
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             !$acc loop vector( 32 )
+             DO  k = nzb+1, nzt
+
+                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) - u_gtrans
+                flux_l              = u_comp * (                           &
+                                               ( 37.0 * ibit2 * adv_sca_5  &
+                                            +     7.0 * ibit1 * adv_sca_3  &
+                                            +           ibit0 * adv_sca_1  &
+                                               ) *                         &
+                                         ( sk(k,j,i)   + sk(k,j,i-1)    )  &
+                                        -      (  8.0 * 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 * ibit2 * adv_sca_5  &
+                                            +     3.0 * ibit1 * adv_sca_3  &
+                                            +           ibit0 * adv_sca_1  &
+                                               ) *                         &
+                                         ( sk(k,j,i)   - sk(k,j,i-1)    )  &
+                                        -      (  5.0 * 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) )     &
+                                                    )
+
+                u_comp    = u(k,j,i+1) - u_gtrans
+                flux_r    = u_comp * (                                      &
+                          ( 37.0 * ibit2 * adv_sca_5                        &
+                      +      7.0 * ibit1 * adv_sca_3                        &
+                      +            ibit0 * adv_sca_1                        &
+                          ) *                                               &
+                             ( sk(k,j,i+1) + sk(k,j,i)   )                  &
+                   -      (  8.0 * 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 * ibit2 * adv_sca_5                        &
+                       +     3.0 * ibit1 * adv_sca_3                        &
+                       +           ibit0 * adv_sca_1                        &
+                          ) *                                               &
+                             ( sk(k,j,i+1) - sk(k,j,i)  )                   &
+                   -      (  5.0 * 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,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,i) - v_gtrans
+                flux_s               = v_comp * (                             &
+                                                  ( 37.0 * ibit5 * adv_sca_5  &
+                                               +     7.0 * ibit4 * adv_sca_3  &
+                                               +           ibit3 * adv_sca_1  &
+                                                  ) *                         &
+                                           ( sk(k,j,i)  + sk(k,j-1,i)     )   &
+                                            -     (  8.0 * 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 * ibit5 * adv_sca_5  &
+                                               +     3.0 * ibit4 * adv_sca_3  &
+                                               +           ibit3 * adv_sca_1  &
+                                                  ) *                         &
+                                            ( sk(k,j,i)   - sk(k,j-1,i)    )  &
+                                           -      (  5.0 * 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) )     &
+                                                     )
+
+
+                v_comp    = v(k,j+1,i) - v_gtrans
+                flux_n    = v_comp * (                                      &
+                          ( 37.0 * ibit5 * adv_sca_5                        &
+                       +     7.0 * ibit4 * adv_sca_3                        &
+                       +           ibit3 * adv_sca_1                        &
+                          ) *                                               &
+                             ( sk(k,j+1,i) + sk(k,j,i)   )                  &
+                   -      (  8.0 * 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 * ibit5 * adv_sca_5                        &
+                       +     3.0 * ibit4 * adv_sca_3                        &
+                       +           ibit3 * adv_sca_1                        &
+                          ) *                                               &
+                             ( sk(k,j+1,i) - sk(k,j,i)    )                 &
+                   -      (  5.0 * 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 * ( 1 - ibit6 )
+                k_mm  = k - 2 * ( 1 - ibit6 )
+                k_mmm = k - 3 * ibit8
+
+                flux_d    = w(k-1,j,i) * (                                    &
+                           ( 37.0 * ibit8 * adv_sca_5                         &
+                        +     7.0 * ibit7 * adv_sca_3                         &
+                        +           ibit6 * adv_sca_1                         &
+                           ) *                                                &
+                                   ( sk(k,j,i)    + sk(k-1,j,i) )             &
+                          -      (  8.0 * 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 * ibit8 * adv_sca_5                         &
+                        +     3.0 * ibit7 * adv_sca_3                         &
+                        +           ibit6 * adv_sca_1                         &
+                           ) *                                                &
+                                   ( sk(k,j,i)    - sk(k-1,j,i)   )           &
+                    -      (  5.0 * 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) * (                                      &
+                           ( 37.0 * ibit8 * adv_sca_5                         &
+                        +     7.0 * ibit7 * adv_sca_3                         &
+                        +           ibit6 * adv_sca_1                         &
+                           ) *                                                &
+                                   ( sk(k+1,j,i)  + sk(k,j,i)   )             &
+                          -      (  8.0 * 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) ) * (                              &
+                           ( 10.0 * ibit8 * adv_sca_5                         &
+                        +     3.0 * ibit7 * adv_sca_3                         &
+                        +           ibit6 * adv_sca_1                         &
+                           ) *                                                &
+                                   ( sk(k+1,j,i)   - sk(k,j,i)    )           &
+                    -      (  5.0 * 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) )           &
+                                         )
+!
+!--             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) - u(k,j,i)   ) * ddx             &
+                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
+                              + ( w(k,j,i)   - w(k-1,j,i) ) * 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 ) * 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)
+!
+!                END SELECT
+
+             ENDDO
+         ENDDO
+      ENDDO
+      !$acc end kernels
+
+    END SUBROUTINE advec_s_ws_acc
+
+
+!------------------------------------------------------------------------------!
 ! Advection of u - Call for all grid points
 !------------------------------------------------------------------------------!
@@ -2752,4 +3047,293 @@
     
     
+!------------------------------------------------------------------------------!
+! Advection of u - Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE advec_u_ws_acc
+
+       USE arrays_3d
+       USE constants
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE statistics
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, ibit9, ibit10, ibit11, ibit12, ibit13, ibit14, ibit15,   &
+                   ibit16, ibit17, j, k, k_mmm, k_mm, k_pp, k_ppp, tn = 0
+
+       REAL    ::  diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div,    &
+                   flux_d, flux_l, flux_n, flux_r, flux_s, flux_t, gu, gv, &
+                   u_comp, u_comp_l, v_comp, v_comp_s, w_comp
+
+
+       gu = 2.0 * u_gtrans
+       gv = 2.0 * v_gtrans
+
+!
+!--    Computation of fluxes and tendency terms
+       !$acc  kernels present( ddzw, tend, u, v, w, wall_flags_0 )
+       !$acc  loop
+       DO i = nxlu, nxr
+          DO  j = nys, nyn
+             !$acc  loop vector( 32 )
+             DO  k = nzb+1, nzt
+
+                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_l           = u(k,j,i) + u(k,j,i-1) - gu
+                flux_l             = u_comp_l * (                          &
+                                    ( 37.0 * ibit11 * adv_mom_5             &
+                                 +     7.0 * ibit10 * adv_mom_3             &
+                                 +           ibit9  * adv_mom_1             &
+                                    ) *                                     &
+                                  ( u(k,j,i)   + u(k,j,i-1) )               &
+                             -      (  8.0 * 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 * ibit11 * adv_mom_5             &
+                                +     3.0 * ibit10 * adv_mom_3             &
+                                +           ibit9  * adv_mom_1             &
+                                   ) *                                     &
+                                 ( u(k,j,i)   - u(k,j,i-1) )               &
+                            -      (  5.0 * 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) )               &
+                                                         )
+
+                u_comp    = u(k,j,i+1) + u(k,j,i)
+                flux_r    = ( u_comp   - gu ) * (                           &
+                          ( 37.0 * ibit11 * adv_mom_5                        &
+                       +     7.0 * ibit10 * adv_mom_3                        &
+                       +           ibit9  * adv_mom_1                        &
+                          ) *                                                &
+                                 ( u(k,j,i+1) + u(k,j,i)   )                 &
+                   -      (  8.0 * 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 * ibit11 * adv_mom_5                        &
+                       +     3.0 * ibit10 * adv_mom_3                        &
+                       +           ibit9  * adv_mom_1                        &
+                          ) *                                                &
+                                 ( u(k,j,i+1) - u(k,j,i)  )                  &
+                   -      (  5.0 * 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,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_s                 = v(k,j,i) + v(k,j,i-1) - gv
+                flux_s                   = v_comp_s * (                       &
+                                   ( 37.0 * ibit14 * adv_mom_5                &
+                                +     7.0 * ibit13 * adv_mom_3                &
+                                +           ibit12 * adv_mom_1                &
+                                   ) *                                        &
+                                     ( u(k,j,i)   + u(k,j-1,i) )              &
+                            -      (  8.0 * 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 * ibit14 * adv_mom_5                &
+                                +     3.0 * ibit13 * adv_mom_3                &
+                                +           ibit12 * adv_mom_1                &
+                                   ) *                                        &
+                                     ( u(k,j,i)   - u(k,j-1,i) )              &
+                            -      (  5.0 * 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) )              &
+                                                         )
+
+
+                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
+                flux_n    = v_comp * (                                       &
+                          ( 37.0 * ibit14 * adv_mom_5                        &
+                       +     7.0 * ibit13 * adv_mom_3                        &
+                       +           ibit12 * adv_mom_1                        &
+                          ) *                                                &
+                                 ( u(k,j+1,i) + u(k,j,i)   )                 &
+                   -      (  8.0 * 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 * ibit14 * adv_mom_5                        &
+                       +     3.0 * ibit13 * adv_mom_3                        &
+                       +           ibit12 * adv_mom_1                        &
+                          ) *                                                &
+                                 ( u(k,j+1,i) - u(k,j,i)  )                  &
+                   -      (  5.0 * 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 * ( 1 - ibit15 )
+                k_mm  = k - 2 * ( 1 - ibit15 )
+                k_mmm = k - 3 * ibit17
+
+                w_comp    = w(k-1,j,i) + w(k-1,j,i-1)
+                flux_d    = w_comp  * (                                      &
+                          ( 37.0 * ibit17 * adv_mom_5                        &
+                       +     7.0 * ibit16 * adv_mom_3                        &
+                       +           ibit15 * adv_mom_1                        &
+                          ) *                                                &
+                             ( u(k,j,i)    + u(k-1,j,i)   )                  &
+                   -      (  8.0 * 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 * ibit17 * adv_mom_5                        &
+                       +     3.0 * ibit16 * adv_mom_3                        &
+                       +           ibit15 * adv_mom_1                        &
+                          ) *                                                &
+                             ( u(k,j,i)     - u(k-1,j,i)  )                  &
+                   -      (  5.0 * 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  * (                                      &
+                          ( 37.0 * ibit17 * adv_mom_5                        &
+                       +     7.0 * ibit16 * adv_mom_3                        &
+                       +           ibit15 * adv_mom_1                        &
+                          ) *                                                &
+                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
+                   -      (  8.0 * 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 ) * (                              &
+                          ( 10.0 * ibit17 * adv_mom_5                        &
+                       +     3.0 * ibit16 * adv_mom_3                        &
+                       +           ibit15 * adv_mom_1                        &
+                          ) *                                                &
+                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
+                   -      (  5.0 * 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      - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx  &
+                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy  &
+                     +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
+                                                                    * ddzw(k)  &
+                      ) * 0.5
+
+                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 ) * 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 * hom(k,1,1,0) )            &
+!                              / ( u_comp    - gu + 1.0E-20      )             &
+!                              +   diss_r    *                                 &
+!                                  ABS( u_comp    - 2.0 * hom(k,1,1,0) )       &
+!                              / ( ABS( u_comp    - gu ) + 1.0E-20 ) )         &
+!                              *   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
+
+
 !------------------------------------------------------------------------------!
 ! Advection of v - Call for all grid points
@@ -3181,4 +3765,295 @@
     
 !------------------------------------------------------------------------------!
+! Advection of v - Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE advec_v_ws_acc
+
+       USE arrays_3d
+       USE constants
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE statistics
+
+       IMPLICIT NONE
+
+
+       INTEGER ::  i, ibit18, ibit19, ibit20, ibit21, ibit22, ibit23, ibit24, &
+                    ibit25, ibit26, j, k, k_mm, k_mmm, k_pp, k_ppp, tn = 0
+
+       REAL    ::  diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div,    &
+                   flux_d, flux_l, flux_n, flux_r, flux_s, flux_t, gu, gv, &
+                   u_comp, u_comp_l, v_comp, v_comp_s, w_comp
+
+       gu = 2.0 * u_gtrans
+       gv = 2.0 * v_gtrans
+
+!
+!--    Computation of fluxes and tendency terms
+       !$acc kernels present( ddzw, tend, u, v, w, wall_flags_0 )
+       !$acc loop
+       DO  i = nxl, nxr
+          DO  j = nysv, nyn
+             !$acc loop vector( 32 )
+             DO  k = nzb+1, nzt
+
+                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_l                 = u(k,j-1,i) + u(k,j,i) - gu
+                flux_l                   = u_comp_l * (                          &
+                                      ( 37.0 * ibit20 * adv_mom_5              &
+                                   +     7.0 * ibit19 * adv_mom_3              &
+                                   +           ibit18 * adv_mom_1              &
+                                      ) *                                      &
+                                     ( v(k,j,i)   + v(k,j,i-1) )               &
+                               -      (  8.0 * 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 * ibit20 * adv_mom_5              &
+                                   +     3.0 * ibit19 * adv_mom_3              &
+                                   +           ibit18 * adv_mom_1              &
+                                      ) *                                      &
+                                     ( v(k,j,i)   - v(k,j,i-1) )               &
+                               -      (  5.0 * 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) )               &
+                                                           )
+
+                u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
+                flux_r    = u_comp * (                                       &
+                          ( 37.0 * ibit20 * adv_mom_5                        &
+                       +     7.0 * ibit19 * adv_mom_3                        &
+                       +           ibit18 * adv_mom_1                        &
+                          ) *                                                &
+                                 ( v(k,j,i+1) + v(k,j,i)   )                 &
+                   -      (  8.0 * 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 * ibit20 * adv_mom_5                        &
+                       +     3.0 * ibit19 * adv_mom_3                        &
+                       +           ibit18 * adv_mom_1                        &
+                          ) *                                                &
+                                 ( v(k,j,i+1) - v(k,j,i)  )                  &
+                   -      (  5.0 * 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,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_s              = v(k,j,i) + v(k,j-1,i) - gv
+                flux_s                = v_comp_s    * (                       &
+                                   ( 37.0 * ibit23 * adv_mom_5                &
+                                +     7.0 * ibit22 * adv_mom_3                &
+                                +           ibit21 * adv_mom_1                &
+                                   ) *                                        &
+                                     ( v(k,j,i)   + v(k,j-1,i) )              &
+                            -      (  8.0 * 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 * ibit23 * adv_mom_5                &
+                                +     3.0 * ibit22 * adv_mom_3                &
+                                +           ibit21 * adv_mom_1                &
+                                   ) *                                        &
+                                     ( v(k,j,i)   - v(k,j-1,i) )              &
+                            -      (  5.0 * 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) )              &
+                                                          )
+
+                v_comp = v(k,j+1,i) + v(k,j,i)
+                flux_n = ( v_comp - gv ) * (                                 &
+                          ( 37.0 * ibit23 * adv_mom_5                        &
+                       +     7.0 * ibit22 * adv_mom_3                        &
+                       +           ibit21 * adv_mom_1                        &
+                          ) *                                                &
+                                 ( v(k,j+1,i) + v(k,j,i)   )                 &
+                   -      (  8.0 * 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 * ibit23 * adv_mom_5                        &
+                       +     3.0 * ibit22 * adv_mom_3                        &
+                       +           ibit21 * adv_mom_1                        &
+                          ) *                                                &
+                                 ( v(k,j+1,i) - v(k,j,i)  )                  &
+                   -      (  5.0 * 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 * ( 1 - ibit24 )
+                k_mm  = k - 2 * ( 1 - ibit24 )
+                k_mmm = k - 3 * ibit26
+
+                w_comp    = w(k-1,j-1,i) + w(k-1,j,i)
+                flux_d    = w_comp  * (                                      &
+                          ( 37.0 * ibit26 * adv_mom_5                        &
+                       +     7.0 * ibit25 * adv_mom_3                        &
+                       +           ibit24 * adv_mom_1                        &
+                          ) *                                                &
+                             ( v(k,j,i)     + v(k-1,j,i)  )                  &
+                   -      (  8.0 * 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 * ibit26 * adv_mom_5                        &
+                       +     3.0 * ibit25 * adv_mom_3                        &
+                       +           ibit24 * adv_mom_1                        &
+                          ) *                                                &
+                             ( v(k,j,i)     - v(k-1,j,i)  )                  &
+                   -      (  5.0 * 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  * (                                      &
+                          ( 37.0 * ibit26 * adv_mom_5                        &
+                       +     7.0 * ibit25 * adv_mom_3                        &
+                       +           ibit24 * adv_mom_1                        &
+                          ) *                                                &
+                             ( v(k+1,j,i)   + v(k,j,i)    )                  &
+                   -      (  8.0 * 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 ) * (                              &
+                          ( 10.0 * ibit26 * adv_mom_5                        &
+                       +     3.0 * ibit25 * adv_mom_3                        &
+                       +           ibit24 * adv_mom_1                        &
+                          ) *                                                &
+                             ( v(k+1,j,i)   - v(k,j,i)    )                  &
+                   -      (  5.0 * 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 - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx &
+                     +  ( v_comp      - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
+                     +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) )         &
+                                                                  ) * ddzw(k) &
+                      ) * 0.5
+
+                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 ) * 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 * hom(k,1,2,0) )                    &
+!                      / ( v_comp - gv + 1.0E-20 )                          &
+!                      +   diss_n                                           &
+!                      *   ABS( v_comp - 2.0 * hom(k,1,2,0) )               &
+!                      / ( ABS( v_comp - gv ) +1.0E-20 ) )                  &
+!                      *   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
+    
+    
+!------------------------------------------------------------------------------!
 ! Advection of w - Call for all grid points
 !------------------------------------------------------------------------------!
@@ -3584,3 +4459,281 @@
     END SUBROUTINE advec_w_ws
 
+
+!------------------------------------------------------------------------------!
+! Advection of w - Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE advec_w_ws_acc
+
+       USE arrays_3d
+       USE constants
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE statistics
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, ibit27, ibit28, ibit29, ibit30, ibit31, ibit32, ibit33, &
+                   ibit34, ibit35, j, k, k_mmm, k_mm, k_pp, k_ppp, tn = 0
+
+       REAL    ::  diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div,    &
+                   flux_d, flux_l, flux_n, flux_r, flux_s, flux_t, gu, gv, &
+                   u_comp, u_comp_l, v_comp, v_comp_s, w_comp
+
+       gu = 2.0 * u_gtrans
+       gv = 2.0 * v_gtrans
+
+!
+!--    Computation of fluxes and tendency terms
+       !$acc kernels present( ddzu, tend, u, v, w, wall_flags_0 )
+       !$acc loop
+       DO i = nxl, nxr
+          DO  j = nys, nyn
+             !$acc loop vector( 32 )
+             DO  k = nzb+1, nzt
+
+                ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
+                ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
+                ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
+
+
+                u_comp_l                 = u(k+1,j,i) + u(k,j,i) - gu
+                flux_l                   = u_comp_l * (                        &
+                                      ( 37.0 * ibit29 * adv_mom_5              &
+                                   +     7.0 * ibit28 * adv_mom_3              &
+                                   +           ibit27 * adv_mom_1              &
+                                      ) *                                      &
+                                     ( w(k,j,i)   + w(k,j,i-1) )               &
+                               -      (  8.0 * 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 * ibit29 * adv_mom_5            &
+                                     +     3.0 * ibit28 * adv_mom_3            &
+                                     +           ibit27 * adv_mom_1            &
+                                        ) *                                    &
+                                     ( w(k,j,i)   - w(k,j,i-1) )               &
+                                 -      (  5.0 * 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) )               &
+                                                            )
+
+                u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
+                flux_r    = u_comp * (                                       &
+                          ( 37.0 * ibit29 * adv_mom_5                        &
+                       +     7.0 * ibit28 * adv_mom_3                        &
+                       +           ibit27 * adv_mom_1                        &
+                          ) *                                                &
+                                 ( w(k,j,i+1) + w(k,j,i)   )                 &
+                   -      (  8.0 * 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 * ibit29 * adv_mom_5                        &
+                       +     3.0 * ibit28 * adv_mom_3                        &
+                       +           ibit27 * adv_mom_1                        &
+                          ) *                                                &
+                                 ( w(k,j,i+1) - w(k,j,i)  )                  &
+                   -      (  5.0 * 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_0(k,j,i),32,1)
+                ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
+                ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
+
+                v_comp_s               = v(k+1,j,i) + v(k,j,i) - gv
+                flux_s                 = v_comp_s * (                         &
+                                    ( 37.0 * ibit32 * adv_mom_5               &
+                                 +     7.0 * ibit31 * adv_mom_3               &
+                                 +           ibit30 * adv_mom_1               &
+                                    ) *                                       &
+                                     ( w(k,j,i)   + w(k,j-1,i) )              &
+                             -      (  8.0 * 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 * ibit32 * adv_mom_5               &
+                                 +     3.0 * ibit31 * adv_mom_3               &
+                                 +           ibit30 * adv_mom_1               &
+                                    ) *                                       &
+                                     ( w(k,j,i)   - w(k,j-1,i) )              &
+                             -      (  5.0 * 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) )              &
+                                                        )
+
+                v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
+                flux_n    = v_comp * (                                       &
+                          ( 37.0 * ibit32 * adv_mom_5                        &
+                       +     7.0 * ibit31 * adv_mom_3                        &
+                       +           ibit30 * adv_mom_1                        &
+                          ) *                                                &
+                                 ( w(k,j+1,i) + w(k,j,i)   )                 &
+                   -      (  8.0 * 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 * ibit32 * adv_mom_5                        &
+                       +     3.0 * ibit31 * adv_mom_3                        &
+                       +           ibit30 * adv_mom_1                        &
+                          ) *                                                &
+                                 ( w(k,j+1,i) - w(k,j,i)  )                  &
+                   -      (  5.0 * 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_0(k-1,j,i),35,1)
+                ibit34 = IBITS(wall_flags_0(k-1,j,i),34,1)
+                ibit33 = IBITS(wall_flags_0(k-1,j,i),33,1)
+
+                k_pp  = k + 2 * ( 1 - ibit33 )
+                k_mm  = k - 2 * ( 1 - ibit33 )
+                k_mmm = k - 3 * ibit35
+
+                w_comp    = w(k,j,i) + w(k-1,j,i)
+                flux_d    = w_comp  * (                                      &
+                          ( 37.0 * ibit35 * adv_mom_5                        &
+                       +     7.0 * ibit34 * adv_mom_3                        &
+                       +           ibit33 * adv_mom_1                        &
+                          ) *                                                &
+                             ( w(k,j,i)    + w(k-1,j,i)   )                  &
+                   -      (  8.0 * 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 * ibit35 * adv_mom_5                        &
+                       +     3.0 * ibit34 * adv_mom_3                        &
+                       +           ibit33 * adv_mom_1                        &
+                          ) *                                                &
+                             ( w(k,j,i)    - w(k-1,j,i)   )                  &
+                   -      (  5.0 * 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_0(k,j,i),35,1)
+                ibit34 = IBITS(wall_flags_0(k,j,i),34,1)
+                ibit33 = IBITS(wall_flags_0(k,j,i),33,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  * (                                      &
+                          ( 37.0 * ibit35 * adv_mom_5                        &
+                       +     7.0 * ibit34 * adv_mom_3                        &
+                       +           ibit33 * adv_mom_1                        &
+                          ) *                                                &
+                             ( w(k+1,j,i)  + w(k,j,i)     )                  &
+                   -      (  8.0 * 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 ) * (                              &
+                          ( 10.0 * ibit35 * adv_mom_5                        &
+                       +     3.0 * ibit34 * adv_mom_3                        &
+                       +           ibit33 * adv_mom_1                        &
+                          ) *                                                &
+                             ( w(k+1,j,i)   - w(k,j,i)    )                  &
+                   -      (  5.0 * 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 - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx  &
+                    +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy  &
+                    +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) )        &
+                                                                 * ddzu(k+1) &
+                      ) * 0.5
+
+                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 ) * 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
