Index: /palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- /palm/trunk/SOURCE/advec_ws.f90	(revision 1014)
+++ /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
Index: /palm/trunk/SOURCE/buoyancy.f90
===================================================================
--- /palm/trunk/SOURCE/buoyancy.f90	(revision 1014)
+++ /palm/trunk/SOURCE/buoyancy.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Currrent revisions:
 ! -----------------
-!
+! accelerator version (*_acc) added
 !
 ! Former revisions:
@@ -55,5 +55,5 @@
 
     PRIVATE
-    PUBLIC buoyancy, calc_mean_profile
+    PUBLIC buoyancy, buoyancy_acc, calc_mean_profile
 
     INTERFACE buoyancy
@@ -61,4 +61,8 @@
        MODULE PROCEDURE buoyancy_ij
     END INTERFACE buoyancy
+
+    INTERFACE buoyancy_acc
+       MODULE PROCEDURE buoyancy_acc
+    END INTERFACE buoyancy_acc
 
     INTERFACE calc_mean_profile
@@ -164,4 +168,107 @@
 
     END SUBROUTINE buoyancy
+
+
+!------------------------------------------------------------------------------!
+! Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE buoyancy_acc( var, var_reference, wind_component, pr )
+
+       USE arrays_3d
+       USE control_parameters
+       USE indices
+       USE pegrid
+       USE statistics
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k, pr, wind_component
+       REAL    ::  var_reference
+#if defined( __nopointer )
+       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
+#else
+       REAL, DIMENSION(:,:,:), POINTER ::  var
+#endif
+
+
+       IF ( .NOT. sloping_surface )  THEN
+!
+!--       Normal case: horizontal surface
+          IF ( use_reference )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
+                                                            (                  &
+                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
+                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
+                                                            )
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             !$acc kernels present( nzb_s_inner, hom, tend, var )
+             !$acc loop
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   !$acc loop vector(32)
+                   DO  k = 1, nzt-1
+                      IF ( k > nzb_s_inner(j,i) )  THEN
+                         tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
+                                                               (                  &
+                             ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
+                             ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
+                                                               )
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDDO
+             !$acc end kernels
+          ENDIF
+
+       ELSE
+!
+!--       Buoyancy term for a surface with a slope in x-direction. The equations
+!--       for both the u and w velocity-component contain proportionate terms.
+!--       Temperature field at time t=0 serves as environmental temperature.
+!--       Reference temperature (pt_surface) is the one at the lower left corner
+!--       of the total domain.
+          IF ( wind_component == 1 )  THEN
+
+             DO  i = nxlu, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
+                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
+                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
+                                 ) / pt_surface
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          ELSEIF ( wind_component == 3 )  THEN
+
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
+                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
+                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
+                                 ) / pt_surface
+                   ENDDO
+                ENDDO
+            ENDDO
+
+          ELSE
+
+             WRITE( message_string, * ) 'no term for component "',&
+                                       wind_component,'"'
+             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
+
+          ENDIF
+
+       ENDIF
+
+    END SUBROUTINE buoyancy_acc
 
 
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 1014)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 1015)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-!
+! acc allowed for loop optimization,
+! checks for adjustment of mixing length to the Prandtl mixing length removed
 !
 ! Former revisions:
@@ -481,10 +482,16 @@
        ENDIF
     ENDIF
-    IF ( loop_optimization /= 'noopt'  .AND.  loop_optimization /= 'cache' &
-         .AND.  loop_optimization /= 'vector' )  THEN
-       message_string = 'illegal value given for loop_optimization: "' // &
-                        TRIM( loop_optimization ) // '"'
-       CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 )
-    ENDIF
+
+    SELECT CASE ( TRIM( loop_optimization ) )
+
+       CASE ( 'acc', 'cache', 'noopt', 'vector' )
+          CONTINUE
+
+       CASE DEFAULT
+          message_string = 'illegal value given for loop_optimization: "' // &
+                           TRIM( loop_optimization ) // '"'
+          CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 )
+
+    END SELECT
 
 !
@@ -1394,12 +1401,8 @@
     IF ( bc_e_b == 'neumann' )  THEN
        ibc_e_b = 1
-       IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-          message_string = 'adjust_mixing_length = TRUE and bc_e_b = "neumann"'
-          CALL message( 'check_parameters', 'PA0055', 0, 1, 0, 6, 0 )
-       ENDIF
     ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
        ibc_e_b = 2
-       IF ( .NOT. adjust_mixing_length  .AND.  prandtl_layer )  THEN
-          message_string = 'adjust_mixing_length = FALSE and bc_e_b = "' // &
+       IF ( prandtl_layer )  THEN
+          message_string = 'adjust mixing length = FALSE and bc_e_b = "' // &
                            TRIM( bc_e_b ) // '"'
           CALL message( 'check_parameters', 'PA0056', 0, 1, 0, 6, 0 )
Index: /palm/trunk/SOURCE/coriolis.f90
===================================================================
--- /palm/trunk/SOURCE/coriolis.f90	(revision 1014)
+++ /palm/trunk/SOURCE/coriolis.f90	(revision 1015)
@@ -4,6 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
-!
+! accelerator version (*_acc) added
 !
 ! Former revisions:
@@ -36,5 +35,5 @@
 
     PRIVATE
-    PUBLIC coriolis
+    PUBLIC coriolis, coriolis_acc
 
     INTERFACE coriolis
@@ -42,4 +41,8 @@
        MODULE PROCEDURE coriolis_ij
     END INTERFACE coriolis
+
+    INTERFACE coriolis_acc
+       MODULE PROCEDURE coriolis_acc
+    END INTERFACE coriolis_acc
 
  CONTAINS
@@ -116,4 +119,94 @@
 
     END SUBROUTINE coriolis
+
+
+!------------------------------------------------------------------------------!
+! Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE coriolis_acc( component )
+
+       USE arrays_3d
+       USE control_parameters
+       USE indices
+       USE pegrid
+
+       IMPLICIT NONE
+
+       INTEGER ::  component, i, j, k
+
+
+!
+!--    Compute Coriolis terms for the three velocity components
+       SELECT CASE ( component )
+
+!
+!--       u-component
+          CASE ( 1 )
+             !$acc  kernels present( nzb_u_inner, tend, v, vg, w )
+             !$acc  loop
+             DO  i = nxlu, nxr
+                DO  j = nys, nyn
+                   !$acc loop vector( 32 )
+                   DO  k = 1, nzt
+                      IF  ( k > nzb_u_inner(j,i) )  THEN
+                         tend(k,j,i) = tend(k,j,i) + f  *    ( 0.25 *          &
+                                      ( v(k,j,i-1) + v(k,j,i) + v(k,j+1,i-1) + &
+                                        v(k,j+1,i) ) - vg(k) )                 &
+                                                - fs *    ( 0.25 *             &
+                                      ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) &
+                                        + w(k,j,i)   )                         &
+                                                             )
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDDO
+             !$acc end kernels
+
+!
+!--       v-component
+          CASE ( 2 )
+             !$acc  kernels present( nzb_v_inner, tend, u, ug )
+             !$acc  loop
+             DO  i = nxl, nxr
+                DO  j = nysv, nyn
+                   !$acc loop vector( 32 )
+                   DO  k = 1, nzt
+                      IF  ( k > nzb_v_inner(j,i) )  THEN
+                         tend(k,j,i) = tend(k,j,i) - f *     ( 0.25 *          &
+                                      ( u(k,j-1,i) + u(k,j,i) + u(k,j-1,i+1) + &
+                                        u(k,j,i+1) ) - ug(k) )
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDDO
+             !$acc end kernels
+
+!
+!--       w-component
+          CASE ( 3 )
+             !$acc  kernels present( nzb_w_inner, tend, u )
+             !$acc  loop
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   !$acc loop vector( 32 )
+                   DO  k = 1, nzt
+                      IF  ( k > nzb_w_inner(j,i) )  THEN
+                         tend(k,j,i) = tend(k,j,i) + fs * 0.25 *             &
+                                      ( u(k,j,i) + u(k+1,j,i) + u(k,j,i+1) + &
+                                        u(k+1,j,i+1) )
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDDO
+             !$acc end kernels
+
+          CASE DEFAULT
+
+             WRITE( message_string, * ) ' wrong component: ', component
+             CALL message( 'coriolis', 'PA0173', 1, 2, 0, 6, 0 )
+
+       END SELECT
+
+    END SUBROUTINE coriolis_acc
 
 
Index: /palm/trunk/SOURCE/cpu_statistics.f90
===================================================================
--- /palm/trunk/SOURCE/cpu_statistics.f90	(revision 1014)
+++ /palm/trunk/SOURCE/cpu_statistics.f90	(revision 1015)
@@ -4,9 +4,12 @@
 ! Current revisions:
 ! -----------------
-! output of handling of ghostpoint exchange
+! output of accelerator board information
 !
 ! Former revisions:
 ! -----------------
 ! $Id$
+!
+! 683 2011-02-09 14:25:15Z raasch
+! output of handling of ghostpoint exchange
 !
 ! 622 2010-12-10 08:08:13Z raasch
@@ -135,8 +138,12 @@
                           numprocs * threads_per_task, pdims(1), pdims(2), &
                           threads_per_task
+       IF ( num_acc_per_node /= 0 )  WRITE ( 18, 108 )  num_acc_per_node
+       WRITE ( 18, 110 )
 #else
        WRITE ( 18, 100 )  TRIM( run_description_header ),        &
                           numprocs * threads_per_task, 1, 1, &
                           threads_per_task
+       IF ( num_acc_per_node /= 0 )  WRITE ( 18, 109 )  num_acc_per_node
+       WRITE ( 18, 110 )
 #endif
        DO
@@ -276,13 +283,5 @@
 
 100 FORMAT (A/11('-')//'CPU measures for ',I5,' PEs (',I5,'(x) * ',I5,'(y', &
-            &') tasks *',I5,' threads):'/ &
-             &'----------------------------------------------------------', &
-             &'------------'//&
-            &'place:                        mean        counts      min  ', &
-             &'     max       rms'/ &
-            &'                           sec.      %                sec. ', &
-             &'     sec.      sec.'/  &
-            &'-----------------------------------------------------------', &
-             &'-------------------')
+            &') tasks *',I5,' threads):')
 
 101 FORMAT (/'special measures:'/ &
@@ -296,4 +295,14 @@
 106 FORMAT (/'Exchange of ghostpoints via MPI_ISEND/MPI_IRECV')
 107 FORMAT (//)
+108 FORMAT ('Accelerator boards per node: ',I2)
+109 FORMAT ('Accelerator boards: ',I2)
+110 FORMAT ('----------------------------------------------------------',   &
+            &'------------'//&
+            &'place:                        mean        counts      min  ', &
+            &'     max       rms'/ &
+            &'                           sec.      %                sec. ', &
+            &'     sec.      sec.'/  &
+            &'-----------------------------------------------------------', &
+            &'-------------------')
 
  END SUBROUTINE cpu_statistics
Index: /palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_e.f90	(revision 1014)
+++ /palm/trunk/SOURCE/diffusion_e.f90	(revision 1015)
@@ -4,5 +4,7 @@
 ! Current revisions:
 ! -----------------
-!
+! accelerator version (*_acc) added,
+! adjustment of mixing length to the Prandtl mixing length at first grid point
+! above ground removed
 !
 ! Former revisions:
@@ -56,5 +58,5 @@
 
     PRIVATE
-    PUBLIC diffusion_e
+    PUBLIC diffusion_e, diffusion_e_acc
     
 
@@ -64,4 +66,8 @@
     END INTERFACE diffusion_e
  
+    INTERFACE diffusion_e_acc
+       MODULE PROCEDURE diffusion_e_acc
+    END INTERFACE diffusion_e_acc
+
  CONTAINS
 
@@ -81,5 +87,5 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  dvar_dz, l_stable, phi_m, var_reference
+       REAL    ::  dvar_dz, l_stable, var_reference
 
 #if defined( __nopointer )
@@ -98,15 +104,4 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-!
-!--             First, calculate phi-function for eventually adjusting the &
-!--             mixing length to the prandtl mixing length
-                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-                   IF ( rif(j,i) >= 0.0 )  THEN
-                      phi_m = 1.0 + 5.0 * rif(j,i)
-                   ELSE
-                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
-                   ENDIF
-                ENDIF
-
                 DO  k = nzb_s_inner(j,i)+1, nzt
 !
@@ -133,12 +128,4 @@
                       ll(k,j) = l_grid(k)
                    ENDIF
-                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-                      l(k,j)  = MIN( l(k,j),  kappa *                          &
-                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
-                                              / phi_m )
-                      ll(k,j) = MIN( ll(k,j), kappa *                          &
-                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
-                                              / phi_m )
-                   ENDIF
 
                 ENDDO
@@ -188,15 +175,4 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-!
-!--             First, calculate phi-function for eventually adjusting the &
-!--             mixing length to the prandtl mixing length
-                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-                   IF ( rif(j,i) >= 0.0 )  THEN
-                      phi_m = 1.0 + 5.0 * rif(j,i)
-                   ELSE
-                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
-                   ENDIF
-                ENDIF
-
                 DO  k = nzb_s_inner(j,i)+1, nzt
 !
@@ -223,12 +199,4 @@
                       ll(k,j) = l_grid(k)
                    ENDIF
-                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-                      l(k,j)  = MIN( l(k,j),  kappa *                          &
-                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
-                                              / phi_m )
-                      ll(k,j) = MIN( ll(k,j), kappa *                          &
-                                              ( zu(k) - zw(nzb_s_inner(j,i)) ) &
-                                              / phi_m )
-                   ENDIF
 
                 ENDDO
@@ -290,7 +258,7 @@
 
 !------------------------------------------------------------------------------!
-! Call for grid point i,j
-!------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
+! Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE diffusion_e_acc( var, var_reference )
 
        USE arrays_3d
@@ -303,5 +271,189 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  dvar_dz, l_stable, phi_m, var_reference
+       REAL    ::  dissipation, dvar_dz, l, ll, l_stable, var_reference
+
+#if defined( __nopointer )
+       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
+#else
+       REAL, DIMENSION(:,:,:), POINTER ::  var
+#endif
+
+
+!
+!--    This if clause must be outside the k-loop because otherwise
+!--    runtime errors occur with -C hopt on NEC
+       IF ( use_reference )  THEN
+          STOP '+++ use_reference in diffusion_e not implemented'
+!          DO  i = nxl, nxr
+!             DO  j = nys, nyn
+!                DO  k = nzb_s_inner(j,i)+1, nzt
+!
+!--                Calculate the mixing length (for dissipation)
+!                   dvar_dz = atmos_ocean_sign * &
+!                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+!                   IF ( dvar_dz > 0.0 ) THEN
+!                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
+!                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
+!                   ELSE
+!                      l_stable = l_grid(k)
+!                   ENDIF
+!
+!--                Adjustment of the mixing length
+!                   IF ( wall_adjustment )  THEN
+!                      l(k,j)  = MIN( wall_adjustment_factor *          &
+!                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+!                                     l_grid(k), l_stable )
+!                      ll(k,j) = MIN( wall_adjustment_factor *          &
+!                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+!                                     l_grid(k) )
+!                   ELSE
+!                      l(k,j)  = MIN( l_grid(k), l_stable )
+!                      ll(k,j) = l_grid(k)
+!                   ENDIF
+!
+!                ENDDO
+!             ENDDO
+!
+!
+!--          Calculate the tendency terms
+!             DO  j = nys, nyn
+!                DO  k = nzb_s_inner(j,i)+1, nzt
+!
+!                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
+!                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
+!
+!                    tend(k,j,i) = tend(k,j,i)                                  &
+!                                        + (                                    &
+!                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
+!                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
+!                                          ) * ddx2                             &
+!                                        + (                                    &
+!                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
+!                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
+!                                          ) * ddy2                             &
+!                                        + (                                    &
+!               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
+!             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
+!                                          ) * ddzw(k)                          &
+!                             - dissipation(k,j)
+!
+!                ENDDO
+!             ENDDO
+!
+!
+!--          Store dissipation if needed for calculating the sgs particle
+!--          velocities
+!             IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
+!                DO  j = nys, nyn
+!                   DO  k = nzb_s_inner(j,i)+1, nzt
+!                      diss(k,j,i) = dissipation(k,j)
+!                   ENDDO
+!                ENDDO
+!             ENDIF
+!
+!          ENDDO
+!
+       ELSE
+
+          !$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
+          !$acc         present( nzb_s_inner, rif, tend, var, zu, zw )
+          !$acc loop
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                !$acc loop vector( 32 )
+                DO  k = 1, nzt
+
+                   IF ( k > nzb_s_inner(j,i) )  THEN
+!
+!--                   Calculate the mixing length (for dissipation)
+                      dvar_dz = atmos_ocean_sign * &
+                                ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                      IF ( dvar_dz > 0.0 ) THEN
+                         l_stable = 0.76 * SQRT( e(k,j,i) ) / &
+                                           SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
+                      ELSE
+                         l_stable = l_grid(k)
+                      ENDIF
+!
+!--                   Adjustment of the mixing length
+                      IF ( wall_adjustment )  THEN
+                         l  = MIN( wall_adjustment_factor *          &
+                                   ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+                                     l_grid(k), l_stable )
+                         ll = MIN( wall_adjustment_factor *          &
+                                   ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+                                   l_grid(k) )
+                      ELSE
+                         l  = MIN( l_grid(k), l_stable )
+                         ll = l_grid(k)
+                      ENDIF
+!
+!--                   Calculate the tendency terms
+                      dissipation = ( 0.19 + 0.74 * l / ll ) * &
+                                    e(k,j,i) * SQRT( e(k,j,i) ) / l
+
+                      tend(k,j,i) = tend(k,j,i)                                &
+                                        + (                                    &
+                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
+                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
+                                          ) * ddx2                             &
+                                        + (                                    &
+                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
+                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
+                                          ) * ddy2                             &
+                                        + (                                    &
+               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
+             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
+                                          ) * ddzw(k)                          &
+                             - dissipation
+
+!
+!--                   Store dissipation if needed for calculating the sgs
+!--                   particle  velocities
+                      IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
+                         diss(k,j,i) = dissipation
+                      ENDIF
+
+                   ENDIF
+
+                ENDDO
+             ENDDO
+          ENDDO
+          !$acc end kernels
+
+       ENDIF
+
+!
+!--    Boundary condition for dissipation
+       IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
+          !$acc kernels present( diss, nzb_s_inner )
+          !$acc loop
+          DO  i = nxl, nxr
+             !$acc loop vector( 32 )
+             DO  j = nys, nyn
+                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
+             ENDDO
+          ENDDO
+          !$acc end kernels
+       ENDIF
+
+    END SUBROUTINE diffusion_e_acc
+
+
+!------------------------------------------------------------------------------!
+! Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE particle_attributes
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  dvar_dz, l_stable, var_reference
 
 #if defined( __nopointer )
@@ -312,15 +464,4 @@
        REAL, DIMENSION(nzb+1:nzt) ::  dissipation, l, ll
 
-
-!
-!--    First, calculate phi-function for eventually adjusting the mixing length
-!--    to the prandtl mixing length
-       IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-          IF ( rif(j,i) >= 0.0 )  THEN
-             phi_m = 1.0 + 5.0 * rif(j,i)
-          ELSE
-             phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
-          ENDIF
-       ENDIF
 
 !
@@ -352,11 +493,4 @@
              ll(k) = l_grid(k)
           ENDIF
-          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-             l(k)  = MIN( l(k),  kappa * &
-                                 ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
-             ll(k) = MIN( ll(k), kappa * &
-                                 ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
-          ENDIF
-
 !
 !--       Calculate the tendency term
Index: /palm/trunk/SOURCE/diffusion_s.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_s.f90	(revision 1014)
+++ /palm/trunk/SOURCE/diffusion_s.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! ------------------
-!
+! accelerator version (*_acc) added
 !
 ! Former revisions:
@@ -48,5 +48,5 @@
 
     PRIVATE
-    PUBLIC diffusion_s
+    PUBLIC diffusion_s, diffusion_s_acc
 
     INTERFACE diffusion_s
@@ -54,4 +54,8 @@
        MODULE PROCEDURE diffusion_s_ij
     END INTERFACE diffusion_s
+
+    INTERFACE diffusion_s_acc
+       MODULE PROCEDURE diffusion_s_acc
+    END INTERFACE diffusion_s_acc
 
  CONTAINS
@@ -173,7 +177,7 @@
 
 !------------------------------------------------------------------------------!
-! Call for grid point i,j
-!------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_s_ij( i, j, s, s_flux_b, s_flux_t, wall_s_flux )
+! Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE diffusion_s_acc( s, s_flux_b, s_flux_t, wall_s_flux )
 
        USE arrays_3d
@@ -194,4 +198,125 @@
 #endif
 
+       !$acc kernels present( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, kh )        &
+       !$acc         present( nzb_diff_s_inner, nzb_s_inner, nzb_s_outer, s ) &
+       !$acc         present( s_flux_b, s_flux_t, tend, wall_s_flux )         &
+       !$acc         present( wall_w_x, wall_w_y )
+       !$acc loop
+       DO  i = nxl, nxr
+          DO  j = nys,nyn
+!
+!--          Compute horizontal diffusion
+             !$acc loop vector( 32 )
+             DO  k = 1, nzt
+                IF ( k > nzb_s_outer(j,i) )  THEN
+
+                   tend(k,j,i) = tend(k,j,i)                                  &
+                                          + 0.5 * (                           &
+                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
+                      - ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
+                                                  ) * ddx2                    &
+                                          + 0.5 * (                           &
+                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
+                      - ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
+                                                  ) * ddy2
+                ENDIF
+             ENDDO
+
+!
+!--          Apply prescribed horizontal wall heatflux where necessary
+             !$acc loop vector(32)
+             DO  k = 1, nzt
+                IF ( k > nzb_s_inner(j,i)  .AND.  k <= nzb_s_outer(j,i)  .AND. &
+                     ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 ) )    &
+                THEN
+                   tend(k,j,i) = tend(k,j,i)                                  &
+                                                + ( fwxp(j,i) * 0.5 *         &
+                        ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
+                        + ( 1.0 - fwxp(j,i) ) * wall_s_flux(1)                &
+                                                   -fwxm(j,i) * 0.5 *         &
+                        ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
+                        + ( 1.0 - fwxm(j,i) ) * wall_s_flux(2)                &
+                                                  ) * ddx2                    &
+                                                + ( fwyp(j,i) * 0.5 *         &
+                        ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
+                        + ( 1.0 - fwyp(j,i) ) * wall_s_flux(3)                &
+                                                   -fwym(j,i) * 0.5 *         &
+                        ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
+                        + ( 1.0 - fwym(j,i) ) * wall_s_flux(4)                &
+                                                  ) * ddy2
+                ENDIF
+             ENDDO
+
+!
+!--          Compute vertical diffusion. In case that surface fluxes have been
+!--          prescribed or computed at bottom and/or top, index k starts/ends at
+!--          nzb+2 or nzt-1, respectively.
+             !$acc loop vector( 32 )
+             DO  k = 1, nzt_diff
+                IF ( k >= nzb_diff_s_inner(j,i) )  THEN
+                   tend(k,j,i) = tend(k,j,i)                                  &
+                                       + 0.5 * (                              &
+            ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
+          - ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)   &
+                                               ) * ddzw(k)
+                ENDIF
+             ENDDO
+
+!
+!--          Vertical diffusion at the first computational gridpoint along
+!--          z-direction
+             !$acc loop vector( 32 )
+             DO  k = 1, nzt
+                IF ( use_surface_fluxes  .AND.  k == nzb_s_inner(j,i)+1 )  THEN
+                   tend(k,j,i) = tend(k,j,i)                                  &
+                                          + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) &
+                                                  * ( s(k+1,j,i)-s(k,j,i) )   &
+                                                  * ddzu(k+1)                 &
+                                              + s_flux_b(j,i)                 &
+                                            ) * ddzw(k)
+                ENDIF
+
+!
+!--             Vertical diffusion at the last computational gridpoint along
+!--             z-direction
+                IF ( use_top_fluxes  .AND.  k == nzt )  THEN
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                                          + ( - s_flux_t(j,i)                  &
+                                              - 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )&
+                                                    * ( s(k,j,i)-s(k-1,j,i) )  &
+                                                    * ddzu(k)                  &
+                                            ) * ddzw(k)
+                ENDIF
+             ENDDO
+
+          ENDDO
+       ENDDO
+       !$acc end kernels
+
+    END SUBROUTINE diffusion_s_acc
+
+
+!------------------------------------------------------------------------------!
+! Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE diffusion_s_ij( i, j, s, s_flux_b, s_flux_t, wall_s_flux )
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  vertical_gridspace
+       REAL    ::  wall_s_flux(0:4)
+       REAL, DIMENSION(nysg:nyng,nxlg:nxrg) ::  s_flux_b, s_flux_t
+#if defined( __nopointer )
+       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s
+#else
+       REAL, DIMENSION(:,:,:), POINTER ::  s
+#endif
+
 !
 !--    Compute horizontal diffusion
Index: /palm/trunk/SOURCE/diffusion_u.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_u.f90	(revision 1014)
+++ /palm/trunk/SOURCE/diffusion_u.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! accelerator version (*_acc) added
 !
 ! Former revisions:
@@ -58,5 +58,5 @@
 
     PRIVATE
-    PUBLIC diffusion_u
+    PUBLIC diffusion_u, diffusion_u_acc
 
     INTERFACE diffusion_u
@@ -64,4 +64,8 @@
        MODULE PROCEDURE diffusion_u_ij
     END INTERFACE diffusion_u
+
+    INTERFACE diffusion_u_acc
+       MODULE PROCEDURE diffusion_u_acc
+    END INTERFACE diffusion_u_acc
 
  CONTAINS
@@ -220,4 +224,188 @@
 
 !------------------------------------------------------------------------------!
+! Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE diffusion_u_acc
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  kmym, kmyp, kmzm, kmzp
+
+       !$acc declare create ( usvs )
+       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs
+
+!
+!--    First calculate horizontal momentum flux u'v' at vertical walls,
+!--    if neccessary
+       IF ( topography /= 'flat' )  THEN
+          CALL wall_fluxes_acc( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
+                                nzb_u_outer, wall_u )
+       ENDIF
+
+       !$acc kernels present ( u, v, w, km, tend, usws, uswst )   &
+       !$acc         present ( ddzu, ddzw, fym, fyp, wall_u )           &
+       !$acc         present ( nzb_u_inner, nzb_u_outer, nzb_diff_u )
+       !$acc loop
+       DO  i = nxlu, nxr
+          DO  j = nys, nyn
+!
+!--          Compute horizontal diffusion
+             !$acc loop vector(32)
+             DO  k = 1, nzt
+                IF ( k > nzb_u_outer(j,i) )  THEN
+!
+!--                Interpolate eddy diffusivities on staggered gridpoints
+                   kmyp = 0.25 * &
+                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
+                   kmym = 0.25 * &
+                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
+
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                         & + 2.0 * (                                           &
+                         &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   ) &
+                         &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) ) &
+                         &         ) * ddx2                                    &
+                         & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy      &
+                         &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx      &
+                         &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
+                         &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
+                         &   ) * ddy
+                ENDIF
+             ENDDO
+
+!
+!--          Wall functions at the north and south walls, respectively
+             !$acc loop vector(32)
+             DO  k = 1, nzt
+                IF( k > nzb_u_inner(j,i)  .AND.  k <= nzb_u_outer(j,i)  .AND. &
+                    wall_u(j,i) /= 0.0 )  THEN
+
+                   kmyp = 0.25 * &
+                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
+                   kmym = 0.25 * &
+                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
+
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                                 + 2.0 * (                                     &
+                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
+                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
+                                         ) * ddx2                              &
+                                 + (   fyp(j,i) * (                            &
+                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
+                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
+                                                  )                            &
+                                     - fym(j,i) * (                            &
+                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
+                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
+                                                  )                            &
+                                     + wall_u(j,i) * usvs(k,j,i)               &
+                                   ) * ddy
+                ENDIF
+             ENDDO
+
+!
+!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
+!--          index k starts at nzb_u_inner+2.
+             !$acc loop vector(32)
+             DO  k = 1, nzt_diff
+                IF ( k >= nzb_diff_u(j,i) )  THEN
+!
+!--                Interpolate eddy diffusivities on staggered gridpoints
+                   kmzp = 0.25 * &
+                          ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
+                   kmzm = 0.25 * &
+                          ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
+
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                         & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)&
+                         &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx      &
+                         &            )                                        &
+                         &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)&
+                         &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx    &
+                         &            )                                        &
+                         &   ) * ddzw(k)
+                ENDIF
+             ENDDO
+
+          ENDDO
+       ENDDO
+
+!
+!--    Vertical diffusion at the first grid point above the surface,
+!--    if the momentum flux at the bottom is given by the Prandtl law or
+!--    if it is prescribed by the user.
+!--    Difference quotient of the momentum flux is not formed over half
+!--    of the grid spacing (2.0*ddzw(k)) any more, since the comparison
+!--    with other (LES) modell showed that the values of the momentum
+!--    flux becomes too large in this case.
+!--    The term containing w(k-1,..) (see above equation) is removed here
+!--    because the vertical velocity is assumed to be zero at the surface.
+       IF ( use_surface_fluxes )  THEN
+
+          !$acc loop
+          DO  i = nxlu, nxr
+             !$acc loop vector(32)
+             DO  j = nys, nyn
+          
+                k = nzb_u_inner(j,i)+1
+!
+!--             Interpolate eddy diffusivities on staggered gridpoints
+                kmzp = 0.25 * &
+                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
+                kmzm = 0.25 * &
+                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
+
+                tend(k,j,i) = tend(k,j,i)                                    &
+                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
+                      &   ) * ddzw(k)                                        &
+                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
+                      &   + usws(j,i)                                        &
+                      &   ) * ddzw(k)
+             ENDDO
+          ENDDO
+
+       ENDIF
+
+!
+!--    Vertical diffusion at the first gridpoint below the top boundary,
+!--    if the momentum flux at the top is prescribed by the user
+       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
+
+          k = nzt
+
+          !$acc loop
+          DO  i = nxlu, nxr
+             !$acc loop vector(32)
+             DO  j = nys, nyn
+
+!
+!--             Interpolate eddy diffusivities on staggered gridpoints
+                kmzp = 0.25 * &
+                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
+                kmzm = 0.25 * &
+                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
+
+                tend(k,j,i) = tend(k,j,i)                                    &
+                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
+                      &   ) * ddzw(k)                                        &
+                      & + ( -uswst(j,i)                                      &
+                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
+                      &   ) * ddzw(k)
+             ENDDO
+          ENDDO
+
+       ENDIF
+       !$acc end kernels
+
+    END SUBROUTINE diffusion_u_acc
+
+
+!------------------------------------------------------------------------------!
 ! Call for grid point i,j
 !------------------------------------------------------------------------------!
Index: /palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_v.f90	(revision 1014)
+++ /palm/trunk/SOURCE/diffusion_v.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! accelerator version (*_acc) added
 !
 ! Former revisions:
@@ -56,5 +56,5 @@
 
     PRIVATE
-    PUBLIC diffusion_v
+    PUBLIC diffusion_v, diffusion_v_acc
 
     INTERFACE diffusion_v
@@ -62,4 +62,8 @@
        MODULE PROCEDURE diffusion_v_ij
     END INTERFACE diffusion_v
+
+    INTERFACE diffusion_v_acc
+       MODULE PROCEDURE diffusion_v_acc
+    END INTERFACE diffusion_v_acc
 
  CONTAINS
@@ -218,4 +222,188 @@
 
 !------------------------------------------------------------------------------!
+! Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE diffusion_v_acc
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  kmxm, kmxp, kmzm, kmzp
+
+       !$acc declare create ( vsus )
+       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus
+
+!
+!--    First calculate horizontal momentum flux v'u' at vertical walls,
+!--    if neccessary
+       IF ( topography /= 'flat' )  THEN
+          CALL wall_fluxes_acc( vsus, 0.0, 1.0, 0.0, 0.0, nzb_v_inner, &
+                                nzb_v_outer, wall_v )
+       ENDIF
+
+       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )   &
+       !$acc         present ( ddzu, ddzw, fxm, fxp, wall_v )           &
+       !$acc         present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
+       !$acc loop
+       DO  i = nxl, nxr
+          DO  j = nysv, nyn
+!
+!--          Compute horizontal diffusion
+             !$acc loop vector(32)
+             DO  k = 1, nzt
+                IF ( k > nzb_v_outer(j,i) )  THEN
+!
+!--                Interpolate eddy diffusivities on staggered gridpoints
+                   kmxp = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
+                   kmxm = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
+
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                         & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx      &
+                         &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy      &
+                         &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
+                         &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
+                         &   ) * ddx                                           &
+                         & + 2.0 * (                                           &
+                         &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )   &
+                         &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )   &
+                         &         ) * ddy2
+                ENDIF
+             ENDDO
+
+!
+!--          Wall functions at the left and right walls, respectively
+             !$acc loop vector(32)
+             DO  k = 1, nzt
+                IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND. &
+                    wall_v(j,i) /= 0.0 )  THEN
+
+                   kmxp = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
+                   kmxm = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
+                   
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                                 + 2.0 * (                                     &
+                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
+                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
+                                         ) * ddy2                              &
+                                 + (   fxp(j,i) * (                            &
+                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
+                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
+                                                  )                            &
+                                     - fxm(j,i) * (                            &
+                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
+                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
+                                                  )                            &
+                                     + wall_v(j,i) * vsus(k,j,i)               &
+                                   ) * ddx
+                ENDIF
+             ENDDO
+
+!
+!--          Compute vertical diffusion. In case of simulating a Prandtl
+!--          layer, index k starts at nzb_v_inner+2.
+             !$acc loop vector(32)
+             DO  k = 1, nzt_diff
+                IF ( k >= nzb_diff_v(j,i) )  THEN
+!
+!--                Interpolate eddy diffusivities on staggered gridpoints
+                   kmzp = 0.25 * &
+                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+                   kmzm = 0.25 * &
+                          ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                         & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)&
+                         &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy      &
+                         &            )                                        &
+                         &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)&
+                         &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy    &
+                         &            )                                        &
+                         &   ) * ddzw(k)
+                ENDIF
+             ENDDO
+
+          ENDDO
+       ENDDO
+
+!
+!--    Vertical diffusion at the first grid point above the surface,
+!--    if the momentum flux at the bottom is given by the Prandtl law
+!--    or if it is prescribed by the user.
+!--    Difference quotient of the momentum flux is not formed over
+!--    half of the grid spacing (2.0*ddzw(k)) any more, since the
+!--    comparison with other (LES) modell showed that the values of
+!--    the momentum flux becomes too large in this case.
+!--    The term containing w(k-1,..) (see above equation) is removed here
+!--    because the vertical velocity is assumed to be zero at the surface.
+       IF ( use_surface_fluxes )  THEN
+
+          !$acc loop
+          DO  i = nxl, nxr
+             !$acc loop vector(32)
+             DO  j = nysv, nyn
+          
+                k = nzb_v_inner(j,i)+1
+!
+!--             Interpolate eddy diffusivities on staggered gridpoints
+                kmzp = 0.25 * &
+                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+                kmzm = 0.25 * &
+                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+
+                tend(k,j,i) = tend(k,j,i)                                    &
+                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
+                      &   ) * ddzw(k)                                        &
+                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
+                      &   + vsws(j,i)                                        &
+                      &   ) * ddzw(k)
+             ENDDO
+          ENDDO
+
+       ENDIF
+
+!
+!--    Vertical diffusion at the first gridpoint below the top boundary,
+!--    if the momentum flux at the top is prescribed by the user
+       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
+
+          k = nzt
+
+          !$acc loop
+          DO  i = nxl, nxr
+             !$acc loop vector(32)
+             DO  j = nysv, nyn
+
+!
+!--             Interpolate eddy diffusivities on staggered gridpoints
+                kmzp = 0.25 * &
+                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+                kmzm = 0.25 * &
+                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+
+                tend(k,j,i) = tend(k,j,i)                                    &
+                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
+                      &   ) * ddzw(k)                                        &
+                      & + ( -vswst(j,i)                                      &
+                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
+                      &   ) * ddzw(k)
+             ENDDO
+          ENDDO
+
+       ENDIF
+       !$acc end kernels
+
+    END SUBROUTINE diffusion_v_acc
+
+
+!------------------------------------------------------------------------------!
 ! Call for grid point i,j
 !------------------------------------------------------------------------------!
Index: /palm/trunk/SOURCE/diffusion_w.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_w.f90	(revision 1014)
+++ /palm/trunk/SOURCE/diffusion_w.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! accelerator version (*_acc) added
 !
 ! Former revisions:
@@ -50,5 +50,5 @@
 
     PRIVATE
-    PUBLIC diffusion_w
+    PUBLIC diffusion_w, diffusion_w_acc
 
     INTERFACE diffusion_w
@@ -56,4 +56,8 @@
        MODULE PROCEDURE diffusion_w_ij
     END INTERFACE diffusion_w
+
+    INTERFACE diffusion_w_acc
+       MODULE PROCEDURE diffusion_w_acc
+    END INTERFACE diffusion_w_acc
 
  CONTAINS
@@ -170,4 +174,124 @@
 
 !------------------------------------------------------------------------------!
+! Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE diffusion_w_acc
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+       REAL    ::  kmxm, kmxp, kmym, kmyp
+
+       !$acc declare create ( wsus, wsvs )
+       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus, wsvs
+
+
+!
+!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
+!--    walls, if neccessary
+       IF ( topography /= 'flat' )  THEN
+          CALL wall_fluxes_acc( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
+                                nzb_w_outer, wall_w_x )
+          CALL wall_fluxes_acc( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
+                                nzb_w_outer, wall_w_y )
+       ENDIF
+
+       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )    &
+       !$acc         present ( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y )           &
+       !$acc         present ( nzb_w_inner, nzb_w_outer )
+       !$acc loop
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             !$acc loop vector( 32 )
+             DO  k = 1, nzt
+                IF ( k > nzb_w_outer(j,i) )  THEN
+!
+!--                Interpolate eddy diffusivities on staggered gridpoints
+                   kmxp = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
+                   kmxm = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
+                   kmyp = 0.25 * &
+                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
+                   kmym = 0.25 * &
+                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+
+                   tend(k,j,i) = tend(k,j,i)                                     &
+                         & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx        &
+                         &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
+                         &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx          &
+                         &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)    &
+                         &   ) * ddx                                             &
+                         & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy        &
+                         &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
+                         &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy          &
+                         &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)    &
+                         &   ) * ddy                                             &
+                         & + 2.0 * (                                             &
+                         &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
+                         & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
+                         &         ) * ddzu(k+1)
+                ENDIF
+             ENDDO
+
+!
+!--          Wall functions at all vertical walls, where necessary
+             !$acc loop vector( 32 )
+             DO  k = 1,nzt
+
+                IF ( k > nzb_w_inner(j,i)  .AND.  k <= nzb_w_outer(j,i)  .AND. &
+                     wall_w_x(j,i) /= 0.0  .AND.  wall_w_y(j,i) /= 0.0 )  THEN
+!
+!--                Interpolate eddy diffusivities on staggered gridpoints
+                   kmxp = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
+                   kmxm = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
+                   kmyp = 0.25 * &
+                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
+                   kmym = 0.25 * &
+                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+
+                   tend(k,j,i) = tend(k,j,i)                                   &
+                                 + (   fwxp(j,i) * (                           &
+                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
+                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
+                                                   )                           &
+                                     - fwxm(j,i) * (                           &
+                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
+                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
+                                                   )                           &
+                                     + wall_w_x(j,i) * wsus(k,j,i)             &
+                                   ) * ddx                                     &
+                                 + (   fwyp(j,i) * (                           &
+                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
+                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
+                                                   )                           &
+                                     - fwym(j,i) * (                           &
+                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
+                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
+                                                   )                           &
+                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
+                                   ) * ddy                                     &
+                                 + 2.0 * (                                     &
+                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
+                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
+                                         ) * ddzu(k+1)
+                ENDIF
+             ENDDO
+
+          ENDDO
+       ENDDO
+       !$acc end kernels
+
+    END SUBROUTINE diffusion_w_acc
+
+
+!------------------------------------------------------------------------------!
 ! Call for grid point i,j
 !------------------------------------------------------------------------------!
Index: /palm/trunk/SOURCE/diffusivities.f90
===================================================================
--- /palm/trunk/SOURCE/diffusivities.f90	(revision 1014)
+++ /palm/trunk/SOURCE/diffusivities.f90	(revision 1015)
@@ -4,4 +4,7 @@
 ! Current revisions:
 ! -----------------
+! OpenACC statements added + code changes required for GPU optimization,
+! adjustment of mixing length to the Prandtl mixing length at first grid point
+! above ground removed
 ! 
 ! Former revisions:
@@ -52,11 +55,7 @@
     INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
 
-    REAL    ::  dvar_dz, l_stable, var_reference
-
-    REAL, SAVE ::  phi_m = 1.0
+    REAL    ::  dvar_dz, l, ll, l_stable, sqrt_e, var_reference
 
     REAL    ::  var(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
-
-    REAL, DIMENSION(1:nzt) ::  l, ll, sqrt_e
 
 
@@ -71,94 +70,90 @@
 !
 !-- Compute the turbulent diffusion coefficient for momentum
-    !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,phi_m,sqrt_e,sr,tn)
+    !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,sqrt_e,sr,tn)
 !$  tn = omp_get_thread_num()
 
+!
+!-- Data declerations for accelerators
+    !$acc data present( dd2zu, e, km, kh, l_grid, l_wall, nzb_s_inner, rif, var )
+    !$acc kernels
+
+!
+!-- Introduce an optional minimum tke
+    IF ( e_min > 0.0 )  THEN
+       !$OMP DO
+       !$acc loop
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             !$acc loop vector( 32 )
+             DO  k = 1, nzt
+                IF ( k > nzb_s_inner(j,i) )  THEN
+                   e(k,j,i) = MAX( e(k,j,i), e_min )
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+
     !$OMP DO
+    !$acc loop
     DO  i = nxlg, nxrg
        DO  j = nysg, nyng
+          !$acc loop vector( 32 )
+          DO  k = 1, nzt
 
+             IF ( k > nzb_s_inner(j,i) )  THEN
+
+                sqrt_e = SQRT( e(k,j,i) )
 !
-!--       Compute the Phi-function for a possible adaption of the mixing length
-!--       to the Prandtl mixing length
-          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-             IF ( rif(j,i) >= 0.0 )  THEN
-                phi_m = 1.0 + 5.0 * rif(j,i)
-             ELSE
-                phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
-             ENDIF
-          ENDIF
-          
+!--             Determine the mixing length
+                dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
+                          ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                IF ( dvar_dz > 0.0 ) THEN
+                   IF ( use_reference )  THEN
+                      l_stable = 0.76 * sqrt_e / &
+                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
+                   ELSE
+                      l_stable = 0.76 * sqrt_e / &
+                                 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
+                   ENDIF
+                ELSE
+                   l_stable = l_grid(k)
+                ENDIF
 !
-!--       Introduce an optional minimum tke
-          IF ( e_min > 0.0 )  THEN
-             DO  k = nzb_s_inner(j,i)+1, nzt
-                e(k,j,i) = MAX( e(k,j,i), e_min )
-             ENDDO
-          ENDIF
+!--             Adjustment of the mixing length
+                IF ( wall_adjustment )  THEN
+                   l  = MIN( l_wall(k,j,i), l_grid(k), l_stable )
+                   ll = MIN( l_wall(k,j,i), l_grid(k) )
+                ELSE
+                   l  = MIN( l_grid(k), l_stable )
+                   ll = l_grid(k)
+                ENDIF
 
-!
-!--       Calculate square root of e in a seperate loop, because it is used
-!--       twice in the next loop (better vectorization)
-          DO  k = nzb_s_inner(j,i)+1, nzt
-             sqrt_e(k) = SQRT( e(k,j,i) )
-          ENDDO
+      !
+      !--       Compute diffusion coefficients for momentum and heat
+                km(k,j,i) = 0.1 * l * sqrt_e
+                kh(k,j,i) = ( 1.0 + 2.0 * l / ll ) * km(k,j,i)
 
-!
-!--       Determine the mixing length
-          DO  k = nzb_s_inner(j,i)+1, nzt
-             dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
-                       ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-             IF ( dvar_dz > 0.0 ) THEN
-                IF ( use_reference )  THEN
-                   l_stable = 0.76 * sqrt_e(k) / &
-                                     SQRT( g / var_reference * dvar_dz ) + 1E-5
-                ELSE
-                   l_stable = 0.76 * sqrt_e(k) / &
-                                     SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
-                ENDIF
-             ELSE
-                l_stable = l_grid(k)
-             ENDIF
-!
-!--          Adjustment of the mixing length
-             IF ( wall_adjustment )  THEN
-                l(k)  = MIN( l_wall(k,j,i), l_grid(k), l_stable )
-                ll(k) = MIN( l_wall(k,j,i), l_grid(k) )
-             ELSE
-                l(k)  = MIN( l_grid(k), l_stable )
-                ll(k) = l_grid(k)
-             ENDIF
-             IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-                l(k)  = MIN( l(k),  kappa * &
-                                    ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
-                ll(k) = MIN( ll(k), kappa * &
-                                    ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
              ENDIF
 
+#if ! defined( __openacc )
 !
-!--          Compute diffusion coefficients for momentum and heat
-             km(k,j,i) = 0.1 * l(k) * sqrt_e(k)
-             kh(k,j,i) = ( 1.0 + 2.0 * l(k) / ll(k) ) * km(k,j,i)
+!++          Statistics still have to be realized for accelerators
+!--          Summation for averaged profile (cf. flow_statistics)
+             DO  sr = 0, statistic_regions
+                sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)
+             ENDDO
+#endif
 
           ENDDO
-
-!
-!--       Summation for averaged profile (cf. flow_statistics)
-!--       (the IF statement still requires a performance check on NEC machines)
-          DO  sr = 0, statistic_regions
-             IF ( rmask(j,i,sr) /= 0.0 .AND.  &
-                  i >= nxl .AND. i <= nxr .AND. j >= nys .AND. j <= nyn )  THEN
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l(k)
-                ENDDO
-             ENDIF
-          ENDDO
-
        ENDDO
     ENDDO
 
+#if ! defined( __openacc )
+!
+!++ Statistics still have to be realized for accelerators
     sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
                                                   ! data output
-
+#endif
     !$OMP END PARALLEL
 
@@ -169,4 +164,5 @@
 !-- values of the diffusivities are not needed
     !$OMP PARALLEL DO
+    !$acc loop
     DO  i = nxlg, nxrg
        DO  j = nysg, nyng
@@ -198,4 +194,6 @@
     ENDIF
 
+    !$acc end kernels
+    !$acc end data
 
  END SUBROUTINE diffusivities
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 1014)
+++ /palm/trunk/SOURCE/header.f90	(revision 1015)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-!
+! output of Aajustment of mixing length to the Prandtl mixing length at first
+! grid point above ground removed
 !
 ! Former revisions:
@@ -1482,5 +1483,4 @@
        IF ( e_min > 0.0 )  WRITE ( io, 454 )  e_min
        IF ( wall_adjustment )  WRITE ( io, 453 )  wall_adjustment_factor
-       IF ( adjust_mixing_length  .AND.  prandtl_layer )  WRITE ( io, 452 )
     ENDIF
 
@@ -1936,5 +1936,4 @@
 451 FORMAT ('    Diffusion coefficients are constant:'/ &
             '    Km = ',F6.2,' m**2/s   Kh = ',F6.2,' m**2/s   Pr = ',F5.2)
-452 FORMAT ('    Mixing length is limited to the Prandtl mixing lenth.')
 453 FORMAT ('    Mixing length is limited to ',F4.2,' * z')
 454 FORMAT ('    TKE is not allowed to fall below ',E9.2,' (m/s)**2')
Index: /palm/trunk/SOURCE/init_1d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_1d_model.f90	(revision 1014)
+++ /palm/trunk/SOURCE/init_1d_model.f90	(revision 1015)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-!
+! adjustment of mixing length to the Prandtl mixing length at first grid point
+! above ground removed
 !
 ! Former revisions:
@@ -113,12 +114,4 @@
           l_black(nzt+1) = l_black(nzt)
 
-       ENDIF
-
-!
-!--    Adjust mixing length to the prandtl mixing length (within the prandtl
-!--    layer)
-       IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-          k = nzb+1
-          l_black(k) = MIN( l_black(k), kappa * zu(k) )
        ENDIF
     ENDIF
@@ -603,17 +596,4 @@
 
 !
-!--          Adjust mixing length to the prandtl mixing length
-             IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-                k = nzb+1
-                IF ( rif1d(k) >= 0.0 )  THEN
-                   l1d(k) = MIN( l1d(k), kappa * zu(k) / ( 1.0 + 5.0 * &
-                                                           rif1d(k) ) )
-                ELSE
-                   l1d(k) = MIN( l1d(k), kappa * zu(k) *          &
-                                  SQRT( SQRT( 1.0 - 16.0 * rif1d(k) ) ) )
-                ENDIF
-             ENDIF
-
-!
 !--          Compute the diffusion coefficients for momentum via the
 !--          corresponding Prandtl-layer relationship and according to
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 1014)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 1015)
@@ -7,5 +7,5 @@
 ! Current revisions:
 ! ------------------
-!
+! mask is set to zero for ghost boundaries
 !
 ! Former revisions:
@@ -1530,5 +1530,9 @@
 !
 !-- Pre-set masks for regional statistics. Default is the total model domain.
+!-- Ghost points are excluded because counting values at the ghost boundaries
+!-- would bias the statistics
     rmask = 1.0
+    rmask(nxlg:nxl-1,:,:) = 0.0;  rmask(nxr+1:nxrg,:,:) = 0.0
+    rmask(:,nysg:nys-1,:) = 0.0;  rmask(:,nyn+1:nyng,:) = 0.0
 
 !
Index: /palm/trunk/SOURCE/init_grid.f90
===================================================================
--- /palm/trunk/SOURCE/init_grid.f90	(revision 1014)
+++ /palm/trunk/SOURCE/init_grid.f90	(revision 1015)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-!
+! lower index for calculating wall_flags_0 set to nzb_w_inner instead of
+! nzb_w_inner+1
 !
 ! Former revisions:
@@ -1253,5 +1254,5 @@
        DO  i = nxl, nxr
           DO  j = nys, nyn
-             DO  k = nzb_w_inner(j,i)+1, nzt
+             DO  k = nzb_w_inner(j,i), nzt
 !
 !--             w component - x-direction
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 1014)
+++ /palm/trunk/SOURCE/modules.f90	(revision 1015)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-!
+! +acc_rank, num_acc_per_node,
+! -adjust_mixing_length
 !
 ! Former revisions:
@@ -592,5 +593,5 @@
                 mask_i_global, mask_j_global, mask_k_global
 
-    LOGICAL ::  adjust_mixing_length = .FALSE., avs_output = .FALSE., &
+    LOGICAL ::  avs_output = .FALSE., &
                 bc_lr_cyc =.TRUE., bc_lr_dirneu = .FALSE., &
                 bc_lr_dirrad = .FALSE., bc_lr_neudir = .FALSE., &
@@ -1335,5 +1336,6 @@
 #endif
     CHARACTER(LEN=5)       ::  myid_char = ''
-    INTEGER                ::  id_inflow = 0, id_recycling = 0, myid = 0,      &
+    INTEGER                ::  acc_rank, id_inflow = 0, id_recycling = 0,      &
+                               myid = 0, num_acc_per_node = 0,                 &
                                target_id, npex = -1, npey = -1, numprocs = 1,  &
                                numprocs_previous_run = -1,                     &
Index: /palm/trunk/SOURCE/palm.f90
===================================================================
--- /palm/trunk/SOURCE/palm.f90	(revision 1014)
+++ /palm/trunk/SOURCE/palm.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! OpenACC statements added + code changes required for GPU optimization
 !
 ! Former revisions:
@@ -77,4 +77,8 @@
     USE statistics
 
+#if defined( __openacc )
+    USE OPENACC
+#endif
+
     IMPLICIT NONE
 
@@ -84,4 +88,7 @@
     CHARACTER (LEN=1) ::  cdum
     INTEGER           ::  i, run_description_header_i(80)
+#if defined( __openacc )
+    REAL, DIMENSION(100) ::  acc_dum
+#endif
 
     version = 'PALM 3.8a'
@@ -102,4 +109,23 @@
 #endif
 
+#if defined( __openacc )
+!
+!-- Get the number of accelerator boards per node and assign the MPI processes
+!-- to these boards
+    num_acc_per_node  = ACC_GET_NUM_DEVICES( ACC_DEVICE_NVIDIA )
+    acc_rank = MOD( myid, num_acc_per_node )
+    CALL ACC_SET_DEVICE_NUM ( acc_rank, ACC_DEVICE_NVIDIA )
+!
+!-- Test output (to be removed later)
+    WRITE (*,'(A,I4,A,I3,A,I3,A,I3)') '*** Connect MPI-Task ', myid,' to CPU ',&
+                                      acc_rank, ' Devices: ', num_acc_per_node,&
+                                      ' connected to:',                        &
+                                      ACC_GET_DEVICE_NUM( ACC_DEVICE_NVIDIA )
+#endif
+!
+!-- Ensure that OpenACC first attaches the GPU devices by copying a dummy data
+!-- region
+    !$acc data copyin( acc_dum )
+
 !
 !-- Initialize measuring of the CPU-time remaining to the run
@@ -177,4 +203,14 @@
     ENDIF
 
+!
+!-- Declare and initialize variables in the accelerator memory with their
+!-- host values
+    !$acc  data copyin( diss, e, e_p, kh, km, pt, pt_p, q, ql, tend, te_m, tpt_m, tu_m, tv_m, tw_m, u, u_p, v, vpt, v_p, w, w_p )          &
+    !$acc       copyin( ddzu, ddzw, dd2zu, l_grid, l_wall, ptdf_x, ptdf_y, pt_init, rdf, rdf_sc, ug, vg, zu, zw )   &
+    !$acc       copyin( hom, qs, qsws, qswst, rif, rif_wall, shf, ts, tswst, us, usws, uswst, vsws, vswst, z0, z0h )      &
+    !$acc       copyin( fxm, fxp, fym, fyp, fwxm, fwxp, fwym, fwyp, nzb_diff_s_inner, nzb_diff_s_outer, nzb_diff_u )       &
+    !$acc       copyin( nzb_diff_v, nzb_s_inner, nzb_s_outer, nzb_u_inner )    &
+    !$acc       copyin( nzb_u_outer, nzb_v_inner, nzb_v_outer, nzb_w_inner )   &
+    !$acc       copyin( nzb_w_outer, wall_heatflux, wall_e_x, wall_e_y, wall_u, wall_v, wall_w_x, wall_w_y, wall_flags_0 )
 !
 !-- Integration of the model equations using timestep-scheme
@@ -242,4 +278,9 @@
 
 !
+!-- Close the OpenACC dummy data region
+    !$acc end data
+    !$acc end data
+
+!
 !-- Take final CPU-time for CPU-time analysis
     CALL cpu_log( log_point(1), 'total', 'stop' )
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 1014)
+++ /palm/trunk/SOURCE/parin.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! -adjust_mixing_length
 !
 ! Former revisions:
@@ -164,5 +164,5 @@
 
 
-    NAMELIST /inipar/  adjust_mixing_length, alpha_surface, bc_e_b, bc_lr, &
+    NAMELIST /inipar/  alpha_surface, bc_e_b, bc_lr, &
                        bc_ns, bc_p_b, bc_p_t, bc_pt_b, bc_pt_t, bc_q_b, &
              bc_q_t,bc_s_b, bc_s_t, bc_sa_t, bc_uv_b, bc_uv_t, &
Index: /palm/trunk/SOURCE/prandtl_fluxes.f90
===================================================================
--- /palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 1014)
+++ /palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! OpenACC statements added
 !
 ! Former revisions:
@@ -65,6 +65,11 @@
 
     INTEGER ::  i, j, k
+    LOGICAL ::  coupled_run
     REAL    ::  a, b, e_q, rifm, uv_total, z_p
 
+!
+!-- Data information for accelerators
+    !$acc data present( e, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt, q, qs ) &
+    !$acc      present( qsws, rif, shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h )
 ! 
 !-- Compute theta*
@@ -74,4 +79,5 @@
 !--    for u* use the value from the previous time step
        !$OMP PARALLEL DO
+       !$acc kernels do
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -90,4 +96,5 @@
 !--    (the Richardson number is still the one from the previous time step)
        !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
+       !$acc kernels do
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -122,4 +129,5 @@
     IF ( .NOT. humidity )  THEN
        !$OMP PARALLEL DO PRIVATE( k, z_p )
+       !$acc kernels do
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -140,4 +148,5 @@
     ELSE
        !$OMP PARALLEL DO PRIVATE( k, z_p )
+       !$acc kernels do
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -162,4 +171,5 @@
 !-- Compute u* at the scalars' grid points
     !$OMP PARALLEL DO PRIVATE( a, b, k, uv_total, z_p )
+    !$acc kernels do
     DO  i = nxl, nxr
        DO  j = nys, nyn
@@ -203,9 +213,13 @@
 !-- Values of us at ghost point locations are needed for the evaluation of usws 
 !-- and vsws.
+    !$acc update host( us )
     CALL exchange_horiz_2d( us )
+    !$acc update device( us )
+
 !
 !-- Compute u'w' for the total model domain.
 !-- First compute the corresponding component of u* and square it.
     !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
+    !$acc kernels do
     DO  i = nxl, nxr
        DO  j = nys, nyn
@@ -245,4 +259,5 @@
 !-- First compute the corresponding component of u* and square it.
     !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
+    !$acc kernels do
     DO  i = nxl, nxr
        DO  j = nys, nyn
@@ -285,4 +300,5 @@
 !--       For a given water flux in the Prandtl layer:
           !$OMP PARALLEL DO
+          !$acc kernels do
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
@@ -291,6 +307,8 @@
           ENDDO
           
-       ELSE          
+       ELSE
+          coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. run_coupled )
           !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
+          !$acc kernels do
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
@@ -302,6 +320,5 @@
 !--             Assume saturation for atmosphere coupled to ocean (but not
 !--             in case of precursor runs)
-                IF ( coupling_mode == 'atmosphere_to_ocean' .AND. run_coupled )&
-                THEN
+                IF ( coupled_run )  THEN
                    e_q = 6.1 * &
                         EXP( 0.07 * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15 ) )
@@ -334,7 +351,13 @@
 !-- Exchange the boundaries for the momentum fluxes (only for sake of 
 !-- completeness)
+    !$acc update host( usws, vsws )
     CALL exchange_horiz_2d( usws )
     CALL exchange_horiz_2d( vsws )
-    IF ( humidity  .OR.  passive_scalar )  CALL exchange_horiz_2d( qsws )
+    !$acc update device( usws, vsws )
+    IF ( humidity  .OR.  passive_scalar )  THEN
+       !$acc update host( qsws )
+       CALL exchange_horiz_2d( qsws )
+       !$acc update device( qsws )
+    ENDIF
 
 !
@@ -342,4 +365,5 @@
     IF ( .NOT. constant_heatflux )  THEN
        !$OMP PARALLEL DO
+       !$acc kernels do
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -353,4 +377,5 @@
     IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar ) ) THEN
        !$OMP PARALLEL DO
+       !$acc kernels do
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -364,4 +389,5 @@
     IF ( ibc_e_b == 2 )  THEN
        !$OMP PARALLEL DO
+       !$acc kernels do
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -375,4 +401,5 @@
     ENDIF
 
+    !$acc end data
 
  END SUBROUTINE prandtl_fluxes
Index: /palm/trunk/SOURCE/production_e.f90
===================================================================
--- /palm/trunk/SOURCE/production_e.f90	(revision 1014)
+++ /palm/trunk/SOURCE/production_e.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! accelerator version (*_acc) added
 !
 ! Former revisions:
@@ -83,5 +83,5 @@
 
     PRIVATE
-    PUBLIC production_e, production_e_init
+    PUBLIC production_e, production_e_acc, production_e_init
 
     LOGICAL, SAVE ::  first_call = .TRUE.
@@ -94,4 +94,8 @@
     END INTERFACE production_e
     
+    INTERFACE production_e_acc
+       MODULE PROCEDURE production_e_acc
+    END INTERFACE production_e_acc
+
     INTERFACE production_e_init
        MODULE PROCEDURE production_e_init
@@ -675,4 +679,655 @@
 
 !------------------------------------------------------------------------------!
+! Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE production_e_acc
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE statistics
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+
+       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
+                   k1, k2, km_neutral, theta, temp
+
+       !$acc declare create ( usvs, vsus, wsus, wsvs )
+       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs, vsus, wsus, wsvs
+!       REAL, DIMENSION(nzb:nzt+1) ::   usvs, vsus, wsus, wsvs
+
+!
+!--    First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
+!--    vertical walls, if neccessary
+!--    CAUTION: results are slightly different from the ij-version!!
+!--    ij-version should be called further below within the ij-loops!!
+       IF ( topography /= 'flat' )  THEN
+          CALL wall_fluxes_e_acc( usvs, 1.0, 0.0, 0.0, 0.0, wall_e_y )
+          CALL wall_fluxes_e_acc( wsvs, 0.0, 0.0, 1.0, 0.0, wall_e_y )
+          CALL wall_fluxes_e_acc( vsus, 0.0, 1.0, 0.0, 0.0, wall_e_x )
+          CALL wall_fluxes_e_acc( wsus, 0.0, 0.0, 0.0, 1.0, wall_e_x )
+       ENDIF
+
+
+!
+!--    Calculate TKE production by shear
+       !$acc kernels present( ddzw, dd2zu, kh, km, nzb_diff_s_inner, nzb_diff_s_outer ) &
+       !$acc         present( nzb_s_inner, nzb_s_outer, pt, q, ql, qsws, qswst, rho )   &
+       !$acc         present( shf, tend, tswst, u, v, vpt, w, wall_e_x, wall_e_y )      &
+       !$acc         copyin( u_0, v_0 )
+       !$acc loop
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             !$acc loop vector( 32 )
+             DO  k = 1, nzt
+
+                IF ( k >= nzb_diff_s_outer(j,i) )  THEN
+
+                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
+                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                   IF ( def < 0.0 )  def = 0.0
+
+                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+       IF ( prandtl_layer )  THEN
+
+!
+!--       Position beneath wall
+!--       (2) - Will allways be executed.
+!--       'bottom and wall: use u_0,v_0 and wall functions'
+          !$acc loop
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                !$acc loop vector( 32 )
+                DO  k = 1, nzt
+
+                   IF ( ( wall_e_x(j,i) /= 0.0 ).OR.( wall_e_y(j,i) /= 0.0 ) ) &
+                   THEN
+
+                      IF ( k == nzb_diff_s_inner(j,i) - 1 )  THEN
+                         dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
+                         dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                        u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+                         dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
+                         dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                        v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+                         dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
+
+                         IF ( wall_e_y(j,i) /= 0.0 )  THEN
+!
+!--                         Inconsistency removed: as the thermal stratification is
+!--                         not taken into account for the evaluation of the wall
+!--                         fluxes at vertical walls, the eddy viscosity km must not
+!--                         be used for the evaluation of the velocity gradients dudy
+!--                         and dwdy
+!--                         Note: The validity of the new method has not yet been
+!--                               shown, as so far no suitable data for a validation
+!--                               has been available
+!                            CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
+!                                                usvs, 1.0, 0.0, 0.0, 0.0 )
+!                            CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
+!                                                wsvs, 0.0, 0.0, 1.0, 0.0 )
+                            km_neutral = kappa *                                    &
+                                        ( usvs(k,j,i)**2 + wsvs(k,j,i)**2 )**0.25 * &
+                                         0.5 * dy
+                            IF ( km_neutral > 0.0 )  THEN
+                               dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral
+                               dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral
+                            ELSE
+                               dudy = 0.0
+                               dwdy = 0.0
+                            ENDIF
+                         ELSE
+                            dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                            dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                         ENDIF
+
+                         IF ( wall_e_x(j,i) /= 0.0 )  THEN
+!
+!--                         Inconsistency removed: as the thermal stratification is
+!--                         not taken into account for the evaluation of the wall
+!--                         fluxes at vertical walls, the eddy viscosity km must not
+!--                         be used for the evaluation of the velocity gradients dvdx
+!--                         and dwdx
+!--                         Note: The validity of the new method has not yet been
+!--                               shown, as so far no suitable data for a validation
+!--                               has been available
+!                            CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
+!                                                vsus, 0.0, 1.0, 0.0, 0.0 )
+!                            CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
+!                                                wsus, 0.0, 0.0, 0.0, 1.0 )
+                            km_neutral = kappa *                                     &
+                                         ( vsus(k,j,i)**2 + wsus(k,j,i)**2 )**0.25 * &
+                                         0.5 * dx
+                            IF ( km_neutral > 0.0 )  THEN
+                               dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral
+                               dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral
+                            ELSE
+                               dvdx = 0.0
+                               dwdx = 0.0
+                            ENDIF
+                         ELSE
+                            dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                            dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                         ENDIF
+
+                         def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                               dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
+                               dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                         IF ( def < 0.0 )  def = 0.0
+
+                         tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
+
+                      ENDIF
+!
+!--                   (3) - will be executed only, if there is at least one level
+!--                   between (2) and (4), i.e. the topography must have a
+!--                   minimum height of 2 dz. Wall fluxes for this case have
+!--                   already been calculated for (2).
+!--                   'wall only: use wall functions'
+
+                      IF ( k >= nzb_diff_s_inner(j,i)  .AND.  &
+                           k <= nzb_diff_s_outer(j,i)-2 )  THEN
+
+                         dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
+                         dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                        u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+                         dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                         dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                        v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+                         dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
+
+                         IF ( wall_e_y(j,i) /= 0.0 )  THEN
+!
+!--                         Inconsistency removed: as the thermal stratification
+!--                         is not taken into account for the evaluation of the
+!--                         wall fluxes at vertical walls, the eddy viscosity km
+!--                         must not be used for the evaluation of the velocity
+!--                         gradients dudy and dwdy
+!--                         Note: The validity of the new method has not yet
+!--                               been shown, as so far no suitable data for a
+!--                               validation has been available
+                            km_neutral = kappa * ( usvs(k,j,i)**2 + &
+                                                   wsvs(k,j,i)**2 )**0.25 * 0.5 * dy
+                            IF ( km_neutral > 0.0 )  THEN
+                               dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral
+                               dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral
+                            ELSE
+                               dudy = 0.0
+                               dwdy = 0.0
+                            ENDIF
+                         ELSE
+                            dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                            dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                         ENDIF
+
+                         IF ( wall_e_x(j,i) /= 0.0 )  THEN
+!
+!--                         Inconsistency removed: as the thermal stratification
+!--                         is not taken into account for the evaluation of the
+!--                         wall fluxes at vertical walls, the eddy viscosity km
+!--                         must not be used for the evaluation of the velocity
+!--                         gradients dvdx and dwdx
+!--                         Note: The validity of the new method has not yet
+!--                               been shown, as so far no suitable data for a
+!--                               validation has been available
+                            km_neutral = kappa * ( vsus(k,j,i)**2 + &
+                                                   wsus(k,j,i)**2 )**0.25 * 0.5 * dx
+                            IF ( km_neutral > 0.0 )  THEN
+                               dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral
+                               dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral
+                            ELSE
+                               dvdx = 0.0
+                               dwdx = 0.0
+                            ENDIF
+                         ELSE
+                            dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                            dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                         ENDIF
+
+                         def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                              dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
+                              dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                         IF ( def < 0.0 )  def = 0.0
+
+                         tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
+
+                      ENDIF
+
+!
+!--                   (4) - will allways be executed.
+!--                   'special case: free atmosphere' (as for case (0))
+                      IF ( k == nzb_diff_s_outer(j,i)-1 )  THEN
+
+                         dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                         dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                          u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                         dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                          u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                         dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                          v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                         dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                         dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                          v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                         dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                          w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                         dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                          w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                         dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                         def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                               dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
+                               dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                         IF ( def < 0.0 )  def = 0.0
+
+                         tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
+
+                      ENDIF
+
+                   ENDIF
+
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       Position without adjacent wall
+!--       (1) - will allways be executed.
+!--       'bottom only: use u_0,v_0'
+          !$acc loop
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                !$acc loop vector( 32 )
+                DO  k = 1, nzt
+
+                   IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
+                   THEN
+
+                      IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
+
+                         dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                         dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                          u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                         dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                          u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
+
+                         dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                          v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                         dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                         dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                          v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
+
+                         dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                          w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                         dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                          w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                         dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                         def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                               dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
+                               dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                         IF ( def < 0.0 )  def = 0.0
+
+                         tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
+
+                      ENDIF
+
+                   ENDIF
+
+                ENDDO
+             ENDDO
+          ENDDO
+
+       ELSEIF ( use_surface_fluxes )  THEN
+
+          !$acc loop
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                !$acc loop vector(32)
+                DO  k = 1, nzt
+
+                   IF ( k == nzb_diff_s_outer(j,i)-1 )  THEN
+
+                      dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
+                      dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
+                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
+                      dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
+                                       u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
+
+                      dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
+                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
+                      dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
+                      dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
+                                       v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
+
+                      dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
+                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
+                      dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
+                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
+                      dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
+
+                      def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
+                            dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
+                            dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
+
+                      IF ( def < 0.0 )  def = 0.0
+
+                      tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
+
+                   ENDIF
+
+                ENDDO
+             ENDDO
+          ENDDO
+
+       ENDIF
+
+!
+!--    If required, calculate TKE production by buoyancy
+       IF ( .NOT. neutral )  THEN
+
+          IF ( .NOT. humidity )  THEN
+
+             IF ( use_reference )  THEN
+
+                IF ( ocean )  THEN
+!
+!--                So far in the ocean no special treatment of density flux
+!--                in the bottom and top surface layer
+                   !$acc loop
+                   DO  i = nxl, nxr
+                      DO  j = nys, nyn
+                         !$acc loop vector( 32 )
+                         DO  k = 1, nzt
+                            IF ( k > nzb_s_inner(j,i) )  THEN
+                               tend(k,j,i) = tend(k,j,i) +                     &
+                                             kh(k,j,i) * g / rho_reference *   &
+                                             ( rho(k+1,j,i) - rho(k-1,j,i) ) * &
+                                             dd2zu(k)
+                            ENDIF
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ELSE
+
+                   !$acc loop
+                   DO  i = nxl, nxr
+                      DO  j = nys, nyn
+                         !$acc loop vector( 32 )
+                         DO  k = 1, nzt_diff
+                            IF ( k >= nzb_diff_s_inner(j,i) )  THEN
+                               tend(k,j,i) = tend(k,j,i) -                   &
+                                             kh(k,j,i) * g / pt_reference *  &
+                                             ( pt(k+1,j,i) - pt(k-1,j,i) ) * &
+                                             dd2zu(k)
+                            ENDIF
+
+                            IF ( k == nzb_diff_s_inner(j,i)-1  .AND.  &
+                                 use_surface_fluxes )  THEN
+                               tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
+                                                           shf(j,i)
+                            ENDIF
+
+                            IF ( k == nzt  .AND.  use_top_fluxes )  THEN
+                               tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
+                                                           tswst(j,i)
+                            ENDIF
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ENDIF
+
+             ELSE
+
+                IF ( ocean )  THEN
+!
+!--                So far in the ocean no special treatment of density flux
+!--                in the bottom and top surface layer
+                   !$acc loop
+                   DO  i = nxl, nxr
+                      DO  j = nys, nyn
+                         !$acc loop vector( 32 )
+                         DO  k = 1, nzt
+                            IF ( k > nzb_s_inner(j,i) )  THEN
+                               tend(k,j,i) = tend(k,j,i) +                     &
+                                             kh(k,j,i) * g / rho(k,j,i) *      &
+                                             ( rho(k+1,j,i) - rho(k-1,j,i) ) * &
+                                             dd2zu(k)
+                            ENDIF
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ELSE
+
+                   !$acc loop
+                   DO  i = nxl, nxr
+                      DO  j = nys, nyn
+                         !$acc loop vector( 32 )
+                         DO  k = 1, nzt_diff
+                            IF( k >= nzb_diff_s_inner(j,i) )  THEN
+                               tend(k,j,i) = tend(k,j,i) -                   &
+                                             kh(k,j,i) * g / pt(k,j,i) *     &
+                                             ( pt(k+1,j,i) - pt(k-1,j,i) ) * &
+                                             dd2zu(k)
+                            ENDIF
+
+                            IF (  k == nzb_diff_s_inner(j,i)-1  .AND.  &
+                                  use_surface_fluxes )  THEN
+                               tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
+                                                           shf(j,i)
+                            ENDIF
+
+                            IF ( k == nzt  .AND.  use_top_fluxes )  THEN
+                               tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
+                                                           tswst(j,i)
+                            ENDIF
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ENDIF
+
+             ENDIF
+
+          ELSE
+!
+!++          This part gives the PGI compiler problems in the previous loop
+!++          even without any acc statements????
+!             STOP '+++ production_e problems with acc-directives'
+!             !acc loop
+!             DO  i = nxl, nxr
+!                DO  j = nys, nyn
+!                   !acc loop vector( 32 )
+!                   DO  k = 1, nzt_diff
+!
+!                      IF ( k >= nzb_diff_s_inner(j,i) )  THEN
+!
+!                         IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+!                            k1 = 1.0 + 0.61 * q(k,j,i)
+!                            k2 = 0.61 * pt(k,j,i)
+!                            tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
+!                                            g / vpt(k,j,i) *                      &
+!                                            ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
+!                                              k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
+!                                            ) * dd2zu(k)
+!                         ELSE IF ( cloud_physics )  THEN
+!                            IF ( ql(k,j,i) == 0.0 )  THEN
+!                               k1 = 1.0 + 0.61 * q(k,j,i)
+!                               k2 = 0.61 * pt(k,j,i)
+!                            ELSE
+!                               theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+!                               temp  = theta * t_d_pt(k)
+!                               k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+!                                          ( q(k,j,i) - ql(k,j,i) ) *          &
+!                                    ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+!                                    ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+!                                    ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+!                               k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+!                            ENDIF
+!                            tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
+!                                            g / vpt(k,j,i) *                      &
+!                                            ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
+!                                              k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
+!                                            ) * dd2zu(k)
+!                         ELSE IF ( cloud_droplets )  THEN
+!                            k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
+!                            k2 = 0.61 * pt(k,j,i)
+!                            tend(k,j,i) = tend(k,j,i) -                          &
+!                                          kh(k,j,i) * g / vpt(k,j,i) *           &
+!                                          ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) +  &
+!                                            k2 * ( q(k+1,j,i) -  q(k-1,j,i) ) -  &
+!                                            pt(k,j,i) * ( ql(k+1,j,i) -          &
+!                                            ql(k-1,j,i) ) ) * dd2zu(k)
+!                         ENDIF
+!
+!                      ENDIF
+!
+!                   ENDDO
+!                ENDDO
+!             ENDDO
+!
+
+!!++          Next two loops are probably very inefficiently parallellized
+!!++          and will require better optimization
+!             IF ( use_surface_fluxes )  THEN
+!
+!                !acc loop
+!                DO  i = nxl, nxr
+!                   DO  j = nys, nyn
+!                      !acc loop vector( 32 )
+!                      DO  k = 1, nzt_diff
+!
+!                         IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
+!
+!                            IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+!                               k1 = 1.0 + 0.61 * q(k,j,i)
+!                               k2 = 0.61 * pt(k,j,i)
+!                            ELSE IF ( cloud_physics )  THEN
+!                               IF ( ql(k,j,i) == 0.0 )  THEN
+!                                  k1 = 1.0 + 0.61 * q(k,j,i)
+!                                  k2 = 0.61 * pt(k,j,i)
+!                               ELSE
+!                                  theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+!                                  temp  = theta * t_d_pt(k)
+!                                  k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+!                                             ( q(k,j,i) - ql(k,j,i) ) *          &
+!                                       ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+!                                       ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+!                                       ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+!                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+!                               ENDIF
+!                            ELSE IF ( cloud_droplets )  THEN
+!                               k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
+!                               k2 = 0.61 * pt(k,j,i)
+!                            ENDIF
+!
+!                            tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
+!                                                  ( k1* shf(j,i) + k2 * qsws(j,i) )
+!                         ENDIF
+!
+!                      ENDDO
+!                   ENDDO
+!                ENDDO
+!
+!             ENDIF
+!
+!             IF ( use_top_fluxes )  THEN
+!
+!                !acc loop
+!                DO  i = nxl, nxr
+!                   DO  j = nys, nyn
+!                      !acc loop vector( 32 )
+!                      DO  k = 1, nzt
+!                         IF ( k == nzt )  THEN
+!
+!                            IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
+!                               k1 = 1.0 + 0.61 * q(k,j,i)
+!                               k2 = 0.61 * pt(k,j,i)
+!                            ELSE IF ( cloud_physics )  THEN
+!                               IF ( ql(k,j,i) == 0.0 )  THEN
+!                                  k1 = 1.0 + 0.61 * q(k,j,i)
+!                                  k2 = 0.61 * pt(k,j,i)
+!                               ELSE
+!                                  theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+!                                  temp  = theta * t_d_pt(k)
+!                                  k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+!                                             ( q(k,j,i) - ql(k,j,i) ) *          &
+!                                       ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+!                                       ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+!                                       ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+!                                  k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+!                               ENDIF
+!                            ELSE IF ( cloud_droplets )  THEN
+!                               k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
+!                               k2 = 0.61 * pt(k,j,i)
+!                            ENDIF
+!
+!                            tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
+!                                                  ( k1* tswst(j,i) + k2 * qswst(j,i) )
+!
+!                         ENDIF
+!
+!                      ENDDO
+!                   ENDDO
+!                ENDDO
+!
+!             ENDIF
+
+          ENDIF
+
+       ENDIF
+       !$acc end kernels
+
+    END SUBROUTINE production_e_acc
+
+
+!------------------------------------------------------------------------------!
 ! Call for grid point i,j
 !------------------------------------------------------------------------------!
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 1014)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 1015)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-!
+! new branch prognostic_equations_acc
+! OpenACC statements added + code changes required for GPU optimization
 !
 ! Former revisions:
@@ -136,5 +137,5 @@
     PRIVATE
     PUBLIC prognostic_equations_noopt, prognostic_equations_cache, &
-           prognostic_equations_vector
+           prognostic_equations_vector, prognostic_equations_acc
 
     INTERFACE prognostic_equations_noopt
@@ -149,4 +150,8 @@
        MODULE PROCEDURE prognostic_equations_vector
     END INTERFACE prognostic_equations_vector
+
+    INTERFACE prognostic_equations_acc
+       MODULE PROCEDURE prognostic_equations_acc
+    END INTERFACE prognostic_equations_acc
 
 
@@ -1900,3 +1905,582 @@
 
 
+ SUBROUTINE prognostic_equations_acc
+
+!------------------------------------------------------------------------------!
+! Version for accelerator boards
+!------------------------------------------------------------------------------!
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=9) ::  time_to_string
+    INTEGER ::  i, j, k, runge_step
+    REAL    ::  sbt
+
+!
+!-- Set switch for intermediate Runge-Kutta step
+    runge_step = 0
+    IF ( timestep_scheme(1:5) == 'runge' )  THEN
+       IF ( intermediate_timestep_count == 1 )  THEN
+          runge_step = 1
+       ELSEIF ( intermediate_timestep_count < &
+                intermediate_timestep_count_max )  THEN
+          runge_step = 2
+       ENDIF
+    ENDIF
+
+!
+!-- Calculate those variables needed in the tendency terms which need
+!-- global communication
+    IF ( .NOT. neutral )  CALL calc_mean_profile( pt, 4 )
+    IF ( ocean         )  CALL calc_mean_profile( rho, 64 )
+    IF ( humidity      )  CALL calc_mean_profile( vpt, 44 )
+    IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
+         intermediate_timestep_count == 1 )  CALL ws_statistics
+
+!
+!-- u-velocity component
+!++ Statistics still not ported to accelerators
+    !$acc update device( hom )
+    CALL cpu_log( log_point(5), 'u-equation', 'start' )
+
+    IF ( timestep_scheme(1:5) == 'runge' )  THEN
+       IF ( ws_scheme_mom )  THEN
+          CALL advec_u_ws_acc
+       ELSE
+          tend = 0.0    ! to be removed later??
+          CALL advec_u_pw
+       ENDIF
+    ELSE
+       CALL advec_u_up
+    ENDIF
+    CALL diffusion_u_acc
+    CALL coriolis_acc( 1 )
+    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
+       CALL buoyancy( pt, pt_reference, 1, 4 )
+    ENDIF
+
+!
+!-- Drag by plant canopy
+    IF ( plant_canopy )  CALL plant_canopy_model( 1 )
+
+!
+!-- External pressure gradient
+    IF ( dp_external )  THEN
+       DO  i = nxlu, nxr
+          DO  j = nys, nyn
+             DO  k = dp_level_ind_b+1, nzt
+                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+
+    CALL user_actions( 'u-tendency' )
+
+!
+!-- Prognostic equation for u-velocity component
+    !$acc kernels present( nzb_u_inner, rdf, tend, tu_m, u, ug, u_p )
+    !$acc loop
+    DO  i = nxlu, nxr
+       DO  j = nys, nyn
+          !$acc loop vector( 32 )
+          DO  k = 1, nzt
+             IF ( k > nzb_u_inner(j,i) )  THEN
+                u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
+                                                  tsc(3) * tu_m(k,j,i) )       &
+                                      - tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
+!
+!--             Tendencies for the next Runge-Kutta step
+                IF ( runge_step == 1 )  THEN
+                   tu_m(k,j,i) = tend(k,j,i)
+                ELSEIF ( runge_step == 2 )  THEN
+                   tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
+                ENDIF
+             ENDIF
+          ENDDO
+       ENDDO
+    ENDDO
+    !$acc end kernels
+
+    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
+    !$acc update host( u_p )
+
+!
+!-- v-velocity component
+    CALL cpu_log( log_point(6), 'v-equation', 'start' )
+
+    IF ( timestep_scheme(1:5) == 'runge' )  THEN
+       IF ( ws_scheme_mom )  THEN
+          CALL advec_v_ws_acc
+       ELSE
+          tend = 0.0    ! to be removed later??
+          CALL advec_v_pw
+       END IF
+    ELSE
+       CALL advec_v_up
+    ENDIF
+    CALL diffusion_v_acc
+    CALL coriolis_acc( 2 )
+
+!
+!-- Drag by plant canopy
+    IF ( plant_canopy )  CALL plant_canopy_model( 2 )
+
+!
+!-- External pressure gradient
+    IF ( dp_external )  THEN
+       DO  i = nxl, nxr
+          DO  j = nysv, nyn
+             DO  k = dp_level_ind_b+1, nzt
+                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+
+    CALL user_actions( 'v-tendency' )
+
+!
+!-- Prognostic equation for v-velocity component
+    !$acc kernels present( nzb_v_inner, rdf, tend, tv_m, v, vg, v_p )
+    !$acc loop
+    DO  i = nxl, nxr
+       DO  j = nysv, nyn
+          !$acc loop vector( 32 )
+          DO  k = 1, nzt
+             IF ( k > nzb_v_inner(j,i) )  THEN
+                v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
+                                                  tsc(3) * tv_m(k,j,i) )       &
+                                      - tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
+!
+!--             Tendencies for the next Runge-Kutta step
+                IF ( runge_step == 1 )  THEN
+                   tv_m(k,j,i) = tend(k,j,i)
+                ELSEIF ( runge_step == 2 )  THEN
+                   tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
+                ENDIF
+             ENDIF
+          ENDDO
+       ENDDO
+    ENDDO
+    !$acc end kernels
+
+    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
+    !$acc update host( v_p )
+
+!
+!-- w-velocity component
+    CALL cpu_log( log_point(7), 'w-equation', 'start' )
+
+    IF ( timestep_scheme(1:5) == 'runge' )  THEN
+       IF ( ws_scheme_mom )  THEN
+          CALL advec_w_ws_acc
+       ELSE
+          tend = 0.0    ! to be removed later??
+          CALL advec_w_pw
+       ENDIF
+    ELSE
+       CALL advec_w_up
+    ENDIF
+    CALL diffusion_w_acc
+    CALL coriolis_acc( 3 )
+
+    IF ( .NOT. neutral )  THEN
+       IF ( ocean )  THEN
+          CALL buoyancy( rho, rho_reference, 3, 64 )
+       ELSE
+          IF ( .NOT. humidity )  THEN
+             CALL buoyancy_acc( pt, pt_reference, 3, 4 )
+          ELSE
+             CALL buoyancy( vpt, pt_reference, 3, 44 )
+          ENDIF
+       ENDIF
+    ENDIF
+
+!
+!-- Drag by plant canopy
+    IF ( plant_canopy )  CALL plant_canopy_model( 3 )
+
+    CALL user_actions( 'w-tendency' )
+
+!
+!-- Prognostic equation for w-velocity component
+    !$acc kernels present( nzb_w_inner, rdf, tend, tw_m, w, w_p )
+    !$acc loop
+    DO  i = nxl, nxr
+       DO  j = nys, nyn
+          !$acc loop vector( 32 )
+          DO  k = 1, nzt-1
+             IF ( k > nzb_w_inner(j,i) )  THEN
+                w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
+                                                  tsc(3) * tw_m(k,j,i) )       &
+                                      - tsc(5) * rdf(k) * w(k,j,i)
+   !
+   !--          Tendencies for the next Runge-Kutta step
+                IF ( runge_step == 1 )  THEN
+                   tw_m(k,j,i) = tend(k,j,i)
+                ELSEIF ( runge_step == 2 )  THEN
+                   tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
+                ENDIF
+             ENDIF
+          ENDDO
+       ENDDO
+    ENDDO
+    !$acc end kernels
+
+    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
+    !$acc update host( w_p )
+
+
+!
+!-- If required, compute prognostic equation for potential temperature
+    IF ( .NOT. neutral )  THEN
+
+       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
+
+!
+!--    pt-tendency terms with communication
+       sbt = tsc(2)
+       IF ( scalar_advec == 'bc-scheme' )  THEN
+
+          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--          Bott-Chlond scheme always uses Euler time step. Thus:
+             sbt = 1.0
+          ENDIF
+          tend = 0.0
+          CALL advec_s_bc( pt, 'pt' )
+
+       ENDIF
+
+!
+!--    pt-tendency terms with no communication
+       IF ( scalar_advec /= 'bc-scheme' )  THEN
+          tend = 0.0
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( ws_scheme_sca )  THEN
+                CALL advec_s_ws_acc( pt, 'pt' )
+             ELSE
+                tend = 0.0    ! to be removed later??
+                CALL advec_s_pw( pt )
+             ENDIF
+          ELSE
+             CALL advec_s_up( pt )
+          ENDIF
+       ENDIF
+
+       CALL diffusion_s_acc( pt, shf, tswst, wall_heatflux )
+
+!
+!--    If required compute heating/cooling due to long wave radiation processes
+       IF ( radiation )  THEN
+          CALL calc_radiation
+       ENDIF
+
+!
+!--    If required compute impact of latent heat due to precipitation
+       IF ( precipitation )  THEN
+          CALL impact_of_latent_heat
+       ENDIF
+
+!
+!--    Consideration of heat sources within the plant canopy
+       IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
+          CALL plant_canopy_model( 4 )
+       ENDIF
+
+!
+!--    If required compute influence of large-scale subsidence/ascent
+       IF ( large_scale_subsidence )  THEN
+          CALL subsidence( tend, pt, pt_init )
+       ENDIF
+
+       CALL user_actions( 'pt-tendency' )
+
+!
+!--    Prognostic equation for potential temperature
+       !$acc kernels present( nzb_s_inner, rdf_sc, ptdf_x, ptdf_y, pt_init ) &
+       !$acc         present( tend, tpt_m, pt, pt_p )
+       !$acc loop
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             !$acc loop vector( 32 )
+             DO  k = 1, nzt
+                IF ( k > nzb_s_inner(j,i) )  THEN
+                   pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +        &
+                                                       tsc(3) * tpt_m(k,j,i) )    &
+                                           - tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
+                                             ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
+!
+!--                Tendencies for the next Runge-Kutta step
+                   IF ( runge_step == 1 )  THEN
+                      tpt_m(k,j,i) = tend(k,j,i)
+                   ELSEIF ( runge_step == 2 )  THEN
+                      tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
+                   ENDIF
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+       !$acc end kernels
+
+       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
+       !$acc update host( pt_p )
+
+    ENDIF
+
+!
+!-- If required, compute prognostic equation for salinity
+    IF ( ocean )  THEN
+
+       CALL cpu_log( log_point(37), 'sa-equation', 'start' )
+
+!
+!--    sa-tendency terms with communication
+       sbt = tsc(2)
+       IF ( scalar_advec == 'bc-scheme' )  THEN
+
+          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--          Bott-Chlond scheme always uses Euler time step. Thus:
+             sbt = 1.0
+          ENDIF
+          tend = 0.0
+          CALL advec_s_bc( sa, 'sa' )
+
+       ENDIF
+
+!
+!--    sa-tendency terms with no communication
+       IF ( scalar_advec /= 'bc-scheme' )  THEN
+          tend = 0.0
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( ws_scheme_sca )  THEN
+                 CALL advec_s_ws( sa, 'sa' )
+             ELSE
+                 CALL advec_s_pw( sa )
+             ENDIF
+          ELSE
+             CALL advec_s_up( sa )
+          ENDIF
+       ENDIF
+
+       CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
+
+       CALL user_actions( 'sa-tendency' )
+
+!
+!--    Prognostic equation for salinity
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +        &
+                                                    tsc(3) * tsa_m(k,j,i) )    &
+                                        - tsc(5) * rdf_sc(k) *                 &
+                                          ( sa(k,j,i) - sa_init(k) )
+                IF ( sa_p(k,j,i) < 0.0 )  sa_p(k,j,i) = 0.1 * sa(k,j,i)
+!
+!--             Tendencies for the next Runge-Kutta step
+                IF ( runge_step == 1 )  THEN
+                   tsa_m(k,j,i) = tend(k,j,i)
+                ELSEIF ( runge_step == 2 )  THEN
+                   tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tsa_m(k,j,i)
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+
+       CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
+
+!
+!--    Calculate density by the equation of state for seawater
+       CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
+       CALL eqn_state_seawater
+       CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
+
+    ENDIF
+
+!
+!-- If required, compute prognostic equation for total water content / scalar
+    IF ( humidity  .OR.  passive_scalar )  THEN
+
+       CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
+
+!
+!--    Scalar/q-tendency terms with communication
+       sbt = tsc(2)
+       IF ( scalar_advec == 'bc-scheme' )  THEN
+
+          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--          Bott-Chlond scheme always uses Euler time step. Thus:
+             sbt = 1.0
+          ENDIF
+          tend = 0.0
+          CALL advec_s_bc( q, 'q' )
+
+       ENDIF
+
+!
+!--    Scalar/q-tendency terms with no communication
+       IF ( scalar_advec /= 'bc-scheme' )  THEN
+          tend = 0.0
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( ws_scheme_sca )  THEN
+                CALL advec_s_ws( q, 'q' )
+             ELSE
+                CALL advec_s_pw( q )
+             ENDIF
+          ELSE
+             CALL advec_s_up( q )
+          ENDIF
+       ENDIF
+
+       CALL diffusion_s( q, qsws, qswst, wall_qflux )
+
+!
+!--    If required compute decrease of total water content due to
+!--    precipitation
+       IF ( precipitation )  THEN
+          CALL calc_precipitation
+       ENDIF
+
+!
+!--    Sink or source of scalar concentration due to canopy elements
+       IF ( plant_canopy ) CALL plant_canopy_model( 5 )
+
+!
+!--    If required compute influence of large-scale subsidence/ascent
+       IF ( large_scale_subsidence )  THEN
+         CALL subsidence( tend, q, q_init )
+       ENDIF
+
+       CALL user_actions( 'q-tendency' )
+
+!
+!--    Prognostic equation for total water content / scalar
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
+                                                  tsc(3) * tq_m(k,j,i) )       &
+                                      - tsc(5) * rdf_sc(k) *                   &
+                                        ( q(k,j,i) - q_init(k) )
+                IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
+!
+!--             Tendencies for the next Runge-Kutta step
+                IF ( runge_step == 1 )  THEN
+                   tq_m(k,j,i) = tend(k,j,i)
+                ELSEIF ( runge_step == 2 )  THEN
+                   tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+
+       CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
+
+    ENDIF
+
+!
+!-- If required, compute prognostic equation for turbulent kinetic
+!-- energy (TKE)
+    IF ( .NOT. constant_diffusion )  THEN
+
+       CALL cpu_log( log_point(16), 'tke-equation', 'start' )
+
+!
+!--    TKE-tendency terms with communication
+       CALL production_e_init
+
+       sbt = tsc(2)
+       IF ( .NOT. use_upstream_for_tke )  THEN
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0
+             ENDIF
+             tend = 0.0
+             CALL advec_s_bc( e, 'e' )
+
+          ENDIF
+       ENDIF
+
+!
+!--    TKE-tendency terms with no communication
+       IF ( scalar_advec /= 'bc-scheme'  .OR.  use_upstream_for_tke )  THEN
+          IF ( use_upstream_for_tke )  THEN
+             tend = 0.0
+             CALL advec_s_up( e )
+          ELSE
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws_acc( e, 'e' )
+                ELSE
+                   tend = 0.0    ! to be removed later??
+                   CALL advec_s_pw( e )
+                ENDIF
+             ELSE
+                tend = 0.0    ! to be removed later??
+                CALL advec_s_up( e )
+             ENDIF
+          ENDIF
+       ENDIF
+
+       IF ( .NOT. humidity )  THEN
+          IF ( ocean )  THEN
+             CALL diffusion_e( prho, prho_reference )
+          ELSE
+             CALL diffusion_e_acc( pt, pt_reference )
+          ENDIF
+       ELSE
+          CALL diffusion_e( vpt, pt_reference )
+       ENDIF
+
+       CALL production_e_acc
+
+!
+!--    Additional sink term for flows through plant canopies
+       IF ( plant_canopy )  CALL plant_canopy_model( 6 )
+       CALL user_actions( 'e-tendency' )
+
+!
+!--    Prognostic equation for TKE.
+!--    Eliminate negative TKE values, which can occur due to numerical
+!--    reasons in the course of the integration. In such cases the old TKE
+!--    value is reduced by 90%.
+       !$acc kernels present( e, e_p, nzb_s_inner, tend, te_m )
+       !$acc loop
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             !$acc loop vector( 32 )
+             DO  k = 1, nzt
+                IF ( k > nzb_s_inner(j,i) )  THEN
+                   e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
+                                                     tsc(3) * te_m(k,j,i) )
+                   IF ( e_p(k,j,i) < 0.0 )  e_p(k,j,i) = 0.1 * e(k,j,i)
+!
+!--                Tendencies for the next Runge-Kutta step
+                   IF ( runge_step == 1 )  THEN
+                      te_m(k,j,i) = tend(k,j,i)
+                   ELSEIF ( runge_step == 2 )  THEN
+                      te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
+                   ENDIF
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+       !$acc end kernels
+
+       CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
+       !$acc update host( e_p )
+
+    ENDIF
+
+
+ END SUBROUTINE prognostic_equations_acc
+
+
  END MODULE prognostic_equations_mod
Index: /palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/read_var_list.f90	(revision 1014)
+++ /palm/trunk/SOURCE/read_var_list.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! ------------------
-!
+! -adjust_mixing_length
 !
 ! Former revisions:
@@ -228,6 +228,4 @@
        SELECT CASE ( TRIM( variable_chr ) )
 
-          CASE ( 'adjust_mixing_length' )
-             READ ( 13 )  adjust_mixing_length
           CASE ( 'advected_distance_x' )
              READ ( 13 )  advected_distance_x
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 1014)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! +call of prognostic_equations_acc
 !
 ! Former revisions:
@@ -176,4 +176,6 @@
           IF ( loop_optimization == 'vector' )  THEN
              CALL prognostic_equations_vector
+          ELSEIF ( loop_optimization == 'acc' )  THEN
+             CALL prognostic_equations_acc
           ELSE
              IF ( scalar_advec == 'bc-scheme' )  THEN
@@ -284,12 +286,23 @@
              CALL pres
           ENDIF
+!
+!--       Update device memory for calculating diffusion quantities and for next
+!--       timestep
+          !$acc update device( e, pt, u, v, w )
+          !$acc update device( q )  if ( allocated( q ) )
 
 !
 !--       If required, compute virtuell potential temperature
-          IF ( humidity ) CALL compute_vpt
+          IF ( humidity )  THEN
+             CALL compute_vpt
+             !$acc update device( vpt )
+          ENDIF
 
 !
 !--       If required, compute liquid water content
-          IF ( cloud_physics ) CALL calc_liquid_water_content
+          IF ( cloud_physics )  THEN
+             CALL calc_liquid_water_content
+             !$acc update device( ql )
+          ENDIF
 
 !
@@ -303,4 +316,7 @@
                 CALL prandtl_fluxes
                 CALL cpu_log( log_point(19), 'prandtl_fluxes', 'stop' )
+!
+!++             Statistics still require updates on host
+                !$acc update host( qs, qsws, rif, shf, ts )
              ENDIF
 
@@ -318,4 +334,7 @@
              ENDIF
              CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
+!
+!++          Statistics still require update of diffusivities on host
+             !$acc update host( kh, km )
 
           ENDIF
Index: /palm/trunk/SOURCE/wall_fluxes.f90
===================================================================
--- /palm/trunk/SOURCE/wall_fluxes.f90	(revision 1014)
+++ /palm/trunk/SOURCE/wall_fluxes.f90	(revision 1015)
@@ -3,5 +3,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! accelerator version (*_acc) added
 !
 ! Former revisions:
@@ -32,5 +32,5 @@
 !------------------------------------------------------------------------------!
     PRIVATE
-    PUBLIC wall_fluxes, wall_fluxes_e
+    PUBLIC wall_fluxes, wall_fluxes_acc, wall_fluxes_e, wall_fluxes_e_acc
     
     INTERFACE wall_fluxes
@@ -39,4 +39,8 @@
     END INTERFACE wall_fluxes
     
+    INTERFACE wall_fluxes_acc
+       MODULE PROCEDURE wall_fluxes_acc
+    END INTERFACE wall_fluxes_acc
+
     INTERFACE wall_fluxes_e
        MODULE PROCEDURE wall_fluxes_e
@@ -44,4 +48,8 @@
     END INTERFACE wall_fluxes_e
  
+    INTERFACE wall_fluxes_e_acc
+       MODULE PROCEDURE wall_fluxes_e_acc
+    END INTERFACE wall_fluxes_e_acc
+
  CONTAINS
 
@@ -189,4 +197,153 @@
     END SUBROUTINE wall_fluxes
 
+
+!------------------------------------------------------------------------------!
+! Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE wall_fluxes_acc( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
+                                nzb_uvw_outer, wall )
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE statistics
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k, max_outer, min_inner, wall_index
+
+       INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) ::  nzb_uvw_inner, &
+                                                   nzb_uvw_outer
+       REAL ::  a, b, c1, c2, h1, h2, zp
+       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
+
+       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
+       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
+
+
+       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
+       wall_flux  = 0.0
+       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
+
+       min_inner = MINVAL( nzb_uvw_inner(nys:nyn,nxl:nxr) ) + 1
+       max_outer = MINVAL( nzb_uvw_outer(nys:nyn,nxl:nxr) )
+
+       !$acc kernels present( hom, nzb_uvw_inner, nzb_uvw_outer, pt, rif_wall ) &
+       !$acc         present( u, v, w, wall, wall_flux, z0 )
+       !$acc loop
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             !$acc loop vector( 32 )
+             DO  k = min_inner, max_outer
+!
+!--             All subsequent variables are computed for the respective
+!--             location where the respective flux is defined.
+                IF ( k >= nzb_uvw_inner(j,i)+1  .AND. &
+                     k <= nzb_uvw_outer(j,i)    .AND.  wall(j,i) /= 0.0 )  THEN
+!
+!--                (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
+                   rifs  = rif_wall(k,j,i,wall_index)
+
+                   u_i   = a * u(k,j,i) + c1 * 0.25 * &
+                           ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
+
+                   v_i   = b * v(k,j,i) + c2 * 0.25 * &
+                           ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
+
+                   ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                   &
+                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
+                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
+                                                           )
+                   pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + &
+                                   b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
+
+                   pts   = pt_i - hom(k,1,4,0)
+                   wspts = ws * pts
+
+!
+!--                (2) Compute wall-parallel absolute velocity vel_total
+                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
+
+!
+!--                (3) Compute wall friction velocity us_wall
+                   IF ( rifs >= 0.0 )  THEN
+
+!
+!--                   Stable stratification (and neutral)
+                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
+                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
+                                                    )
+                   ELSE
+
+!
+!--                   Unstable stratification
+                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+                      us_wall = kappa * vel_total / (                          &
+                           LOG( zp / z0(j,i) ) -                               &
+                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                    )
+                   ENDIF
+
+!
+!--                (4) Compute zp/L (corresponds to neutral Richardson flux
+!--                    number rifs)
+                   rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
+                                                        ( us_wall**3 + 1E-30 ) )
+
+!
+!--                Limit the value range of the Richardson numbers.
+!--                This is necessary for very small velocities (u,w --> 0),
+!--                because the absolute value of rif can then become very
+!--                large, which in consequence would result in very large
+!--                shear stresses and very small momentum fluxes (both are
+!--                generally unrealistic).
+                   IF ( rifs < rif_min )  rifs = rif_min
+                   IF ( rifs > rif_max )  rifs = rif_max
+
+!
+!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
+                   IF ( rifs >= 0.0 )  THEN
+
+!
+!--                   Stable stratification (and neutral)
+                      wall_flux(k,j,i) = kappa *                               &
+                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
+                              (  LOG( zp / z0(j,i) ) +                         &
+                                 5.0 * rifs * ( zp - z0(j,i) ) / zp            &
+                              )
+                   ELSE
+
+!
+!--                   Unstable stratification
+                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+                      wall_flux(k,j,i) = kappa *                               &
+                           ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
+                           LOG( zp / z0(j,i) ) -                               &
+                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                                            )
+                   ENDIF
+                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
+
+!
+!--                store rifs for next time step
+                   rif_wall(k,j,i,wall_index) = rifs
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+       ENDDO
+       !$acc end kernels
+
+    END SUBROUTINE wall_fluxes_acc
 
 
@@ -444,4 +601,134 @@
     END SUBROUTINE wall_fluxes_e
 
+
+!------------------------------------------------------------------------------!
+! Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE wall_fluxes_e_acc( wall_flux, a, b, c1, c2, wall )
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Calculates momentum fluxes at vertical walls for routine production_e
+! assuming Monin-Obukhov similarity.
+! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
+!------------------------------------------------------------------------------!
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE statistics
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k, kk, max_outer, min_inner, wall_index
+       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
+                   ws, zp
+
+       REAL ::  rifs
+
+       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
+       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
+
+
+       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
+       wall_flux  = 0.0
+       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
+
+       min_inner = MINVAL( nzb_diff_s_inner(nys:nyn,nxl:nxr) ) - 1
+       max_outer = MAXVAL( nzb_diff_s_outer(nys:nyn,nxl:nxr) ) - 2
+
+       !$acc kernels present( nzb_diff_s_inner, nzb_diff_s_outer, pt, rif_wall ) &
+       !$acc         present( u, v, w, wall, wall_flux, z0 )
+       !$acc loop
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             !$acc loop vector(32)
+             DO  k = min_inner, max_outer
+!
+!--             All subsequent variables are computed for scalar locations
+                IF ( k >= nzb_diff_s_inner(j,i)-1  .AND. &
+                     k <= nzb_diff_s_outer(j,i)-2  .AND.  wall(j,i) /= 0.0 )  THEN
+!
+!--                (1) Compute rifs, u_i, v_i, and ws
+                   IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
+                      kk = nzb_diff_s_inner(j,i)-1
+                   ELSE
+                      kk = k-1
+                   ENDIF
+                   rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                &
+                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
+                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
+                                 )
+
+                   u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
+                   v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
+                   ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
+!
+!--                (2) Compute wall-parallel absolute velocity vel_total and
+!--                interpolate appropriate velocity component vel_zp.
+                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
+                   vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
+!
+!--                (3) Compute wall friction velocity us_wall
+                   IF ( rifs >= 0.0 )  THEN
+
+!
+!--                   Stable stratification (and neutral)
+                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
+                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
+                                                    )
+                   ELSE
+
+!
+!--                   Unstable stratification
+                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+                      us_wall = kappa * vel_total / (                          &
+                           LOG( zp / z0(j,i) ) -                               &
+                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                    )
+                   ENDIF
+
+!
+!--                Skip step (4) of wall_fluxes, because here rifs is already
+!--                available from (1)
+!
+!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
+
+                   IF ( rifs >= 0.0 )  THEN
+
+!
+!--                   Stable stratification (and neutral)
+                      wall_flux(k,j,i) = kappa *  vel_zp / &
+                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
+                   ELSE
+
+!
+!--                   Unstable stratification
+                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
+                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
+
+                      wall_flux(k,j,i) = kappa * vel_zp / (                    &
+                           LOG( zp / z0(j,i) ) -                               &
+                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
+                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
+                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
+                                                          )
+                   ENDIF
+                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+       ENDDO
+       !$acc end kernels
+
+    END SUBROUTINE wall_fluxes_e_acc
 
 
Index: /palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/write_var_list.f90	(revision 1014)
+++ /palm/trunk/SOURCE/write_var_list.f90	(revision 1015)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! -adjust_mixing_length
 !
 ! Former revisions:
@@ -150,6 +150,4 @@
 !--          list in read_var_list.
 
-    WRITE ( 14 )  'adjust_mixing_length          '
-    WRITE ( 14 )  adjust_mixing_length
     WRITE ( 14 )  'advected_distance_x           '
     WRITE ( 14 )  advected_distance_x
