Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 3633)
+++ palm/trunk/SOURCE/Makefile	(revision 3634)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# OpenACC port for SPEC
+# 
+# 3579 2018-11-29 15:32:39Z suehring
 # Dependency for check_parameters on nesting_offl_mod added
 # 
@@ -554,4 +557,5 @@
 	coriolis.f90 \
 	cpulog_mod.f90 \
+	cuda_fft_interfaces.f90 \
 	data_log.f90 \
 	data_output_2d.f90 \
@@ -919,4 +923,8 @@
 	mod_kinds.o \
 	modules.o
+cuda_fft_interfaces.o: \
+	cuda_fft_interfaces.f90 \
+	modules.o \
+	mod_kinds.o
 data_log.o: \
 	mod_kinds.o \
@@ -1042,4 +1050,5 @@
 	pmc_interface_mod.o
 fft_xy_mod.o: \
+	cuda_fft_interfaces.o \
 	mod_kinds.o \
 	modules.o \
Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 3633)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3589 2018-11-30 15:09:51Z suehring
 ! Move the control parameter "salsa" from salsa_mod to control_parameters 
 ! (M. Kurppa)
@@ -1143,4 +1146,6 @@
 !--    beginning of prognostic_equations.
        IF ( ws_scheme_mom )  THEN
+          !$ACC KERNELS PRESENT(sums_wsus_ws_l, sums_wsvs_ws_l) &
+          !$ACC PRESENT(sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l)
           sums_wsus_ws_l = 0.0_wp
           sums_wsvs_ws_l = 0.0_wp
@@ -1148,8 +1153,11 @@
           sums_vs2_ws_l  = 0.0_wp
           sums_ws2_ws_l  = 0.0_wp
+          !$ACC END KERNELS
        ENDIF
 
        IF ( ws_scheme_sca )  THEN
+          !$ACC KERNELS PRESENT(sums_wspts_ws_l)
           sums_wspts_ws_l = 0.0_wp
+          !$ACC END KERNELS
           IF ( humidity       )  sums_wsqs_ws_l = 0.0_wp
           IF ( passive_scalar )  sums_wsss_ws_l = 0.0_wp
@@ -3300,4 +3308,5 @@
 
        CHARACTER (LEN = *), INTENT(IN)    ::  sk_char !< string identifier, used for assign fluxes to the correct dimension in the analysis array 
+       INTEGER(iwp) ::  sk_num !< integer identifier, used for assign fluxes to the correct dimension in the analysis array
        
        INTEGER(iwp) ::  i      !< grid index along x-direction
@@ -3320,7 +3329,17 @@
        REAL(wp) ::  ibit1  !< flag indicating 3rd-order scheme along x-direction
        REAL(wp) ::  ibit2  !< flag indicating 5th-order scheme along x-direction
+#ifdef _OPENACC
+       REAL(wp) ::  ibit0_l  !< flag indicating 1st-order scheme along x-direction
+       REAL(wp) ::  ibit1_l  !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp) ::  ibit2_l  !< flag indicating 5th-order scheme along x-direction
+#endif
        REAL(wp) ::  ibit3  !< flag indicating 1st-order scheme along y-direction
        REAL(wp) ::  ibit4  !< flag indicating 3rd-order scheme along y-direction
        REAL(wp) ::  ibit5  !< flag indicating 5th-order scheme along y-direction
+#ifdef _OPENACC
+       REAL(wp) ::  ibit3_s  !< flag indicating 1st-order scheme along y-direction
+       REAL(wp) ::  ibit4_s  !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp) ::  ibit5_s  !< flag indicating 5th-order scheme along y-direction
+#endif
        REAL(wp) ::  ibit6  !< flag indicating 1st-order scheme along z-direction
        REAL(wp) ::  ibit7  !< flag indicating 3rd-order scheme along z-direction
@@ -3330,20 +3349,57 @@
        REAL(wp) ::  flux_d !< 6th-order flux at grid box bottom
        REAL(wp) ::  u_comp !< advection velocity along x-direction
+#ifdef _OPENACC
+       REAL(wp) ::  u_comp_l !< advection velocity along x-direction
+#endif
        REAL(wp) ::  v_comp !< advection velocity along y-direction
+#ifdef _OPENACC
+       REAL(wp) ::  v_comp_s !< advection velocity along y-direction
+#endif
        
-       REAL(wp), DIMENSION(nzb:nzt)   ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt)   ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt)   ::  diss_t !< discretized artificial dissipation at rightward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt)   ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt)   ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt)   ::  flux_t !< discretized 6th-order flux at rightward-side of the grid box
+       REAL(wp) ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
+       REAL(wp) ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
+       REAL(wp) ::  diss_t !< discretized artificial dissipation at rightward-side of the grid box
+       REAL(wp) ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
+       REAL(wp) ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
+       REAL(wp) ::  flux_t !< discretized 6th-order flux at rightward-side of the grid box
        
+       REAL(wp) ::  diss_s !< discretized artificial dissipation term at southward-side of the grid box
+       REAL(wp) ::  flux_s !< discretized 6th-order flux at northward-side of the grid box
+#ifndef _OPENACC
        REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_diss_y_local !< discretized artificial dissipation term at southward-side of the grid box
        REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_flux_y_local !< discretized 6th-order flux at northward-side of the grid box
+#endif
        
+       REAL(wp) ::  diss_l !< discretized artificial dissipation term at leftward-side of the grid box
+       REAL(wp) ::  flux_l !< discretized 6th-order flux at leftward-side of the grid box
+#ifndef _OPENACC
        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_diss_x_local !< discretized artificial dissipation term at leftward-side of the grid box
        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_flux_x_local !< discretized 6th-order flux at leftward-side of the grid box
+#endif
         
-
+       SELECT CASE ( sk_char )
+
+           CASE ( 'pt' )
+              sk_num = 1
+           CASE ( 'sa' )
+              sk_num = 2
+           CASE ( 'q' )
+              sk_num = 3
+           CASE ( 'qc' )
+              sk_num = 4
+           CASE ( 'qr' )
+              sk_num = 5
+           CASE ( 'nc' )
+              sk_num = 6
+           CASE ( 'nr' )
+              sk_num = 7
+           CASE ( 's' )
+              sk_num = 8
+           CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' )
+              sk_num = 9
+
+       END SELECT
+
+#ifndef _OPENACC
 !
 !--    Compute the fluxes for the whole left boundary of the processor domain.
@@ -3408,7 +3464,30 @@
 
        ENDDO
-
+#endif
+
+       !$ACC PARALLEL LOOP COLLAPSE(2) FIRSTPRIVATE(tn, sk_num) &
+       !$ACC PRIVATE(i, j, k, k_mm, k_pp, k_ppp) &
+       !$ACC PRIVATE(ibit0, ibit1, ibit2, ibit3, ibit4, ibit5) &
+       !$ACC PRIVATE(ibit0_l, ibit1_l, ibit2_l) &
+       !$ACC PRIVATE(ibit3_s, ibit4_s, ibit5_s) &
+       !$ACC PRIVATE(ibit6, ibit7, ibit8) &
+       !$ACC PRIVATE(flux_r, diss_r, flux_l, diss_l) &
+       !$ACC PRIVATE(flux_n, diss_n, flux_s, diss_s) &
+       !$ACC PRIVATE(flux_t, diss_t, flux_d, diss_d) &
+       !$ACC PRIVATE(div, u_comp, u_comp_l, v_comp, v_comp_s) &
+       !$ACC PRESENT(advc_flags_1) &
+       !$ACC PRESENT(sk, u, v, w, u_stokes_zu, v_stokes_zu) &
+       !$ACC PRESENT(drho_air, rho_air_zw, ddzw) &
+       !$ACC PRESENT(tend) &
+       !$ACC PRESENT(hom(nzb+1:nzb_max,1,1:3,0)) &
+       !$ACC PRESENT(weight_substep(intermediate_timestep_count)) &
+       !$ACC PRESENT(sums_wspts_ws_l, sums_wssas_ws_l) &
+       !$ACC PRESENT(sums_wsqs_ws_l, sums_wsqcs_ws_l) &
+       !$ACC PRESENT(sums_wsqrs_ws_l, sums_wsncs_ws_l) &
+       !$ACC PRESENT(sums_wsnrs_ws_l, sums_wsss_ws_l) &
+       !$ACC PRESENT(sums_salsa_ws_l)
        DO  i = nxl, nxr
 
+#ifndef _OPENACC
           j = nys
           DO  k = nzb+1, nzb_max
@@ -3467,9 +3546,8 @@
 
           ENDDO
+#endif
 
           DO  j = nys, nyn
 
-             flux_t(0) = 0.0_wp
-             diss_t(0) = 0.0_wp
              flux_d    = 0.0_wp
              diss_d    = 0.0_wp
@@ -3481,6 +3559,6 @@
                 ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp )
 
-                u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
-                flux_r(k) = u_comp * (                                        &
+                u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
+                flux_r = u_comp * (                                           &
                           ( 37.0_wp * ibit2 * adv_sca_5                       &
                       +      7.0_wp * ibit1 * adv_sca_3                       &
@@ -3497,5 +3575,5 @@
                                      )
 
-                diss_r(k) = -ABS( u_comp ) * (                                &
+                diss_r = -ABS( u_comp ) * (                                   &
                           ( 10.0_wp * ibit2 * adv_sca_5                       &
                        +     3.0_wp * ibit1 * adv_sca_3                       &
@@ -3512,10 +3590,52 @@
                                              )
 
+#ifdef _OPENACC
+!
+!--             Recompute the left fluxes.
+                ibit2_l = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp )
+                ibit1_l = REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp )
+                ibit0_l = REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp )
+
+                u_comp_l = u(k,j,i) - u_gtrans + u_stokes_zu(k)
+                flux_l = u_comp_l * (                                         &
+                                             ( 37.0_wp * ibit2_l * adv_sca_5  &
+                                          +     7.0_wp * ibit1_l * adv_sca_3  &
+                                          +              ibit0_l * adv_sca_1  &
+                                             ) *                              &
+                                          ( sk(k,j,i)   + sk(k,j,i-1)    )    &
+                                      -      (  8.0_wp * ibit2_l * adv_sca_5  &
+                                          +              ibit1_l * adv_sca_3  &
+                                             ) *                              &
+                                          ( sk(k,j,i+1) + sk(k,j,i-2)    )    &
+                                      +      (           ibit2_l * adv_sca_5  & 
+                                             ) *                              &
+                                          ( sk(k,j,i+2) + sk(k,j,i-3)    )    &
+                                               )
+
+                 diss_l = -ABS( u_comp_l ) * (                                &
+                                             ( 10.0_wp * ibit2_l * adv_sca_5  &
+                                          +     3.0_wp * ibit1_l * adv_sca_3  &
+                                          +              ibit0_l * adv_sca_1  &
+                                             ) *                              &
+                                          ( sk(k,j,i)   - sk(k,j,i-1) )       &
+                                      -      (  5.0_wp * ibit2_l * adv_sca_5  &
+                                          +              ibit1_l * adv_sca_3  &
+                                             ) *                              &
+                                         ( sk(k,j,i+1) - sk(k,j,i-2)  )       &
+                                      +      (           ibit2_l * adv_sca_5  &
+                                             ) *                              &
+                                          ( sk(k,j,i+2) - sk(k,j,i-3) )       &
+                                                        )
+#else
+                flux_l = swap_flux_x_local(k,j)
+                diss_l = swap_diss_x_local(k,j)
+#endif
+
                 ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp )
                 ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp )
                 ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp )
 
-                v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
-                flux_n(k) = v_comp * (                                        &
+                v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
+                flux_n = v_comp * (                                           &
                           ( 37.0_wp * ibit5 * adv_sca_5                       &
                        +     7.0_wp * ibit4 * adv_sca_3                       &
@@ -3532,5 +3652,5 @@
                                      )
 
-                diss_n(k) = -ABS( v_comp ) * (                                &
+                diss_n = -ABS( v_comp ) * (                                   &
                           ( 10.0_wp * ibit5 * adv_sca_5                       &
                        +     3.0_wp * ibit4 * adv_sca_3                       &
@@ -3546,4 +3666,47 @@
                              ( sk(k,j+3,i) - sk(k,j-2,i) )                    &
                                              )
+
+#ifdef _OPENACC
+!
+!--             Recompute the south fluxes.
+                ibit5_s = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp )
+                ibit4_s = REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp )
+                ibit3_s = REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp )
+
+                v_comp_s = v(k,j,i) - v_gtrans + v_stokes_zu(k)
+                flux_s = v_comp_s * (                                         &
+                                             ( 37.0_wp * ibit5_s * adv_sca_5  &
+                                          +     7.0_wp * ibit4_s * adv_sca_3  &
+                                          +              ibit3_s * adv_sca_1  &
+                                             ) *                              &
+                                         ( sk(k,j,i)  + sk(k,j-1,i)     )     &
+                                       -     (  8.0_wp * ibit5_s * adv_sca_5  &
+                                          +              ibit4_s * adv_sca_3  &
+                                              ) *                             &
+                                         ( sk(k,j+1,i) + sk(k,j-2,i)    )     &
+                                      +      (           ibit5_s * adv_sca_5  &
+                                             ) *                              &
+                                        ( sk(k,j+2,i) + sk(k,j-3,i)     )     &
+                                             )
+
+                diss_s = -ABS( v_comp_s ) * (                                 &
+                                             ( 10.0_wp * ibit5_s * adv_sca_5  &
+                                          +     3.0_wp * ibit4_s * adv_sca_3  &
+                                          +              ibit3_s * adv_sca_1  &
+                                             ) *                              &
+                                          ( sk(k,j,i)   - sk(k,j-1,i)    )    &
+                                      -      (  5.0_wp * ibit5_s * adv_sca_5  &
+                                          +              ibit4_s * adv_sca_3  &
+                                             ) *                              &
+                                          ( sk(k,j+1,i) - sk(k,j-2,i)    )    &
+                                      +      (           ibit5_s * adv_sca_5  &
+                                             ) *                              &
+                                          ( sk(k,j+2,i) - sk(k,j-3,i)    )    &
+                                                     )
+#else
+                flux_s = swap_flux_y_local(k)
+                diss_s = swap_diss_y_local(k)
+#endif
+
 !
 !--             k index has to be modified near bottom and top, else array
@@ -3558,5 +3721,5 @@
 
 
-                flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                      &
+                flux_t = w(k,j,i) * rho_air_zw(k) * (                         &
                            ( 37.0_wp * ibit8 * adv_sca_5                      &
                         +     7.0_wp * ibit7 * adv_sca_3                      &
@@ -3572,5 +3735,5 @@
                                        )
 
-                diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (              &
+                diss_t = -ABS( w(k,j,i) ) * rho_air_zw(k) * (                 &
                            ( 10.0_wp * ibit8 * adv_sca_5                      &
                         +     3.0_wp * ibit7 * adv_sca_3                      &
@@ -3616,19 +3779,117 @@
 
                 tend(k,j,i) = tend(k,j,i) - (                                 &
-                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
-                          swap_diss_x_local(k,j)            ) * ddx           &
-                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
-                          swap_diss_y_local(k)              ) * ddy           &
-                      + ( ( flux_t(k) + diss_t(k) ) -                         &
-                          ( flux_d    + diss_d    )                           &
-                                                    ) * drho_air(k) * ddzw(k) &
+                        ( ( flux_r + diss_r )                                 &
+                      -   ( flux_l + diss_l ) ) * ddx                         &
+                      + ( ( flux_n + diss_n )                                 &
+                      -   ( flux_s + diss_s ) ) * ddy                         &
+                      + ( ( flux_t + diss_t )                                 &
+                      -   ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k)       &
                                             ) + sk(k,j,i) * div
 
-                swap_flux_y_local(k)   = flux_n(k)
-                swap_diss_y_local(k)   = diss_n(k)
-                swap_flux_x_local(k,j) = flux_r(k)
-                swap_diss_x_local(k,j) = diss_r(k)
-                flux_d                 = flux_t(k)
-                diss_d                 = diss_t(k)
+#ifndef _OPENACC
+                swap_flux_y_local(k)   = flux_n
+                swap_diss_y_local(k)   = diss_n
+                swap_flux_x_local(k,j) = flux_r
+                swap_diss_x_local(k,j) = diss_r
+#endif
+                flux_d                 = flux_t
+                diss_d                 = diss_t
+
+!
+!--             Evaluation of statistics. 
+                SELECT CASE ( sk_num )
+
+                    CASE ( 1 )
+                       !$ACC ATOMIC
+                       sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)           &
+                          + ( flux_t                                           &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t                                           &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
+                    CASE ( 2 )
+                       !$ACC ATOMIC
+                       sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)           &
+                          + ( flux_t                                           &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t                                           &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
+                    CASE ( 3 )
+                       !$ACC ATOMIC
+                       sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn)            &
+                          + ( flux_t                                           &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t                                           &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
+                    CASE ( 4 )
+                       !$ACC ATOMIC
+                       sums_wsqcs_ws_l(k,tn)  = sums_wsqcs_ws_l(k,tn)          &
+                          + ( flux_t                                           &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t                                           &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
+                    CASE ( 5 )
+                       !$ACC ATOMIC
+                       sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn)          &
+                          + ( flux_t                                           &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t                                           &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
+                    CASE ( 6 )
+                       !$ACC ATOMIC
+                       sums_wsncs_ws_l(k,tn)  = sums_wsncs_ws_l(k,tn)          &
+                          + ( flux_t                                           &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t                                           &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
+                    CASE ( 7 )
+                       !$ACC ATOMIC
+                       sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn)          &
+                          + ( flux_t                                           &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t                                           &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
+                    CASE ( 8 )
+                       !$ACC ATOMIC
+                       sums_wsss_ws_l(k,tn)  = sums_wsss_ws_l(k,tn)            &
+                          + ( flux_t                                           &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t                                           &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
+                    CASE ( 9 )
+                        !$ACC ATOMIC
+                        sums_salsa_ws_l(k,tn)  = sums_salsa_ws_l(k,tn)         &
+                          + ( flux_t                                           &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t                                           &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
+
+                END SELECT
 
              ENDDO
@@ -3636,23 +3897,63 @@
              DO  k = nzb_max+1, nzt
 
-                u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
-                flux_r(k) = u_comp * (                                        &
+                u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
+                flux_r = u_comp * (                                           &
                       37.0_wp * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
                     -  8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
                     +           ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
-                diss_r(k) = -ABS( u_comp ) * (                                &
+                diss_r = -ABS( u_comp ) * (                                   &
                       10.0_wp * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
                     -  5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
                     +           ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
 
-                v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
-                flux_n(k) = v_comp * (                                        &
+#ifdef _OPENACC
+!
+!--             Recompute the left fluxes.
+                u_comp_l = u(k,j,i) - u_gtrans + u_stokes_zu(k)
+                flux_l = u_comp_l * (                                         &
+                                      37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) ) &
+                                    -  8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
+                                    +           ( sk(k,j,i+2) + sk(k,j,i-3) ) &
+                                               ) * adv_sca_5
+
+                diss_l = -ABS( u_comp_l ) * (                                 &
+                                      10.0_wp * ( sk(k,j,i)   - sk(k,j,i-1) ) &
+                                    -  5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
+                                    +           ( sk(k,j,i+2) - sk(k,j,i-3) ) &
+                                                       ) * adv_sca_5
+#else
+                flux_l = swap_flux_x_local(k,j)
+                diss_l = swap_diss_x_local(k,j)
+#endif
+
+                v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
+                flux_n = v_comp * (                                           &
                       37.0_wp * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
                     -  8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
                     +           ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
-                diss_n(k) = -ABS( v_comp ) * (                                &
+                diss_n = -ABS( v_comp ) * (                                   &
                       10.0_wp * ( sk(k,j+1,i) - sk(k,j,i)   )                 &
                     -  5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
                     +           ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
+
+#ifdef _OPENACC
+!
+!--             Recompute the south fluxes.
+                v_comp_s = v(k,j,i) - v_gtrans + v_stokes_zu(k)
+                flux_s = v_comp_s * (                                        &
+                                    37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )  &
+                                  -  8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) )  &
+                                  +           ( sk(k,j+2,i) + sk(k,j-3,i) )  &
+                                             ) * adv_sca_5
+                diss_s = -ABS( v_comp_s ) * (                                &
+                                    10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )  &
+                                  -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )  &
+                                  +             sk(k,j+2,i) - sk(k,j-3,i)    &
+                                                      ) * adv_sca_5
+#else
+                flux_s = swap_flux_y_local(k)
+                diss_s = swap_diss_y_local(k)
+#endif
+
 !
 !--             k index has to be modified near bottom and top, else array
@@ -3667,5 +3968,5 @@
 
 
-                flux_t(k) = w(k,j,i) * rho_air_zw(k) * (                      &
+                flux_t = w(k,j,i) * rho_air_zw(k) * (                      &
                            ( 37.0_wp * ibit8 * adv_sca_5                      &
                         +     7.0_wp * ibit7 * adv_sca_3                      &
@@ -3681,5 +3982,5 @@
                                        )
 
-                diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (              &
+                diss_t = -ABS( w(k,j,i) ) * rho_air_zw(k) * (              &
                            ( 10.0_wp * ibit8 * adv_sca_5                      &
                         +     3.0_wp * ibit7 * adv_sca_3                      &
@@ -3706,128 +4007,119 @@
 
                 tend(k,j,i) = tend(k,j,i) - (                                 &
-                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
-                          swap_diss_x_local(k,j)            ) * ddx           &
-                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
-                          swap_diss_y_local(k)              ) * ddy           &
-                      + ( ( flux_t(k) + diss_t(k) ) -                         &
-                          ( flux_d    + diss_d    )                           &
-                                                    ) * drho_air(k) * ddzw(k) &
+                        ( ( flux_r + diss_r )                                 &
+                      -   ( flux_l + diss_l ) ) * ddx                         &
+                      + ( ( flux_n + diss_n )                                 &
+                      -   ( flux_s + diss_s ) ) * ddy                         &
+                      + ( ( flux_t + diss_t )                                 &
+                      -   ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k)       &
                                             ) + sk(k,j,i) * div
 
-                swap_flux_y_local(k)   = flux_n(k)
-                swap_diss_y_local(k)   = diss_n(k)
-                swap_flux_x_local(k,j) = flux_r(k)
-                swap_diss_x_local(k,j) = diss_r(k)
-                flux_d                 = flux_t(k)
-                diss_d                 = diss_t(k)
-
-             ENDDO
-!
-!--          Evaluation of statistics. 
-             SELECT CASE ( sk_char )
-
-                 CASE ( 'pt' )
-                    DO  k = nzb, nzt
+#ifndef _OPENACC
+                swap_flux_y_local(k)   = flux_n
+                swap_diss_y_local(k)   = diss_n
+                swap_flux_x_local(k,j) = flux_r
+                swap_diss_x_local(k,j) = diss_r
+#endif
+                flux_d                 = flux_t
+                diss_d                 = diss_t
+
+!
+!--             Evaluation of statistics. 
+                SELECT CASE ( sk_num )
+
+                    CASE ( 1 )
+                       !$ACC ATOMIC
                        sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)           &
-                          + ( flux_t(k)                                        &
+                          + ( flux_t                                           &
                                 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
-                            + diss_t(k)                                        &
+                            + diss_t                                           &
                                 / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
                                 *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
                             ) * weight_substep(intermediate_timestep_count)
-                    ENDDO
-                 CASE ( 'sa' )
-                    DO  k = nzb, nzt
+                    CASE ( 2 )
+                       !$ACC ATOMIC
                        sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)           &
-                          + ( flux_t(k)                                        &
+                          + ( flux_t                                           &
                                 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
-                            + diss_t(k)                                        &
+                            + diss_t                                           &
                                 / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
                                 *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
                             ) * weight_substep(intermediate_timestep_count)
-                    ENDDO
-                 CASE ( 'q' )
-                    DO  k = nzb, nzt
+                    CASE ( 3 )
+                       !$ACC ATOMIC
                        sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn)            &
-                          + ( flux_t(k)                                        &
+                          + ( flux_t                                           &
                                 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
-                            + diss_t(k)                                        &
+                            + diss_t                                           &
                                 / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
                                 *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
                             ) * weight_substep(intermediate_timestep_count)
-                    ENDDO
-                 CASE ( 'qc' )
-                    DO  k = nzb, nzt
+                    CASE ( 4 )
+                       !$ACC ATOMIC
                        sums_wsqcs_ws_l(k,tn)  = sums_wsqcs_ws_l(k,tn)          &
-                          + ( flux_t(k)                                        &
+                          + ( flux_t                                           &
                                 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
-                            + diss_t(k)                                        &
+                            + diss_t                                           &
                                 / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
                                 *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
                             ) * weight_substep(intermediate_timestep_count)
-                    ENDDO
-                 CASE ( 'qr' )
-                    DO  k = nzb, nzt
+                    CASE ( 5 )
+                       !$ACC ATOMIC
                        sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn)          &
-                          + ( flux_t(k)                                        &
+                          + ( flux_t                                           &
                                 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
-                            + diss_t(k)                                        &
+                            + diss_t                                           &
                                 / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
                                 *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
                             ) * weight_substep(intermediate_timestep_count)
-                    ENDDO
-                 CASE ( 'nc' )
-                    DO  k = nzb, nzt
+                    CASE ( 6 )
+                       !$ACC ATOMIC
                        sums_wsncs_ws_l(k,tn)  = sums_wsncs_ws_l(k,tn)          &
-                          + ( flux_t(k)                                        &
+                          + ( flux_t                                           &
                                 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
-                            + diss_t(k)                                        &
+                            + diss_t                                           &
                                 / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
                                 *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
                             ) * weight_substep(intermediate_timestep_count)
-                    ENDDO
-                 CASE ( 'nr' )
-                    DO  k = nzb, nzt
+                    CASE ( 7 )
+                       !$ACC ATOMIC
                        sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn)          &
-                          + ( flux_t(k)                                        &
+                          + ( flux_t                                           &
                                 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
-                            + diss_t(k)                                        &
+                            + diss_t                                           &
                                 / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
                                 *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
                             ) * weight_substep(intermediate_timestep_count)
-                    ENDDO
-                 CASE ( 's' )
-                    DO  k = nzb, nzt
+                    CASE ( 8 )
+                       !$ACC ATOMIC
                        sums_wsss_ws_l(k,tn)  = sums_wsss_ws_l(k,tn)            &
-                          + ( flux_t(k)                                        &
+                          + ( flux_t                                           &
                                 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
-                            + diss_t(k)                                        &
+                            + diss_t                                           &
                                 / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
                                 *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
                             ) * weight_substep(intermediate_timestep_count)
-                    ENDDO  
-     
-                 CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' )
-                     DO  k = nzb, nzt
+                    CASE ( 9 )
+                        !$ACC ATOMIC
                         sums_salsa_ws_l(k,tn)  = sums_salsa_ws_l(k,tn)         &
-                          + ( flux_t(k)                                        &
+                          + ( flux_t                                           &
                                 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
-                            + diss_t(k)                                        &
+                            + diss_t                                           &
                                 / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
                                 *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
                             ) * weight_substep(intermediate_timestep_count)
-                     ENDDO                    
-                                   
-
-              END SELECT
+
+                END SELECT
+
+             ENDDO
 
          ENDDO
@@ -3874,7 +4166,17 @@
        REAL(wp)    ::  ibit10 !< flag indicating 3rd-order scheme along x-direction
        REAL(wp)    ::  ibit11 !< flag indicating 5th-order scheme along x-direction
+#ifdef _OPENACC
+       REAL(wp)    ::  ibit9_l  !< flag indicating 1st-order scheme along x-direction
+       REAL(wp)    ::  ibit10_l !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp)    ::  ibit11_l !< flag indicating 5th-order scheme along x-direction
+#endif
        REAL(wp)    ::  ibit12 !< flag indicating 1st-order scheme along y-direction
        REAL(wp)    ::  ibit13 !< flag indicating 3rd-order scheme along y-direction
        REAL(wp)    ::  ibit14 !< flag indicating 5th-order scheme along y-direction
+#ifdef _OPENACC
+       REAL(wp)    ::  ibit12_s !< flag indicating 1st-order scheme along y-direction
+       REAL(wp)    ::  ibit13_s !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp)    ::  ibit14_s !< flag indicating 5th-order scheme along y-direction
+#endif
        REAL(wp)    ::  ibit15 !< flag indicating 1st-order scheme along z-direction
        REAL(wp)    ::  ibit16 !< flag indicating 3rd-order scheme along z-direction
@@ -3886,23 +4188,38 @@
        REAL(wp)    ::  gv     !< Galilei-transformation velocity along y 
        REAL(wp)    ::  v_comp !< advection velocity along y
+#ifdef _OPENACC
+       REAL(wp)    ::  v_comp_s !< advection velocity along y
+#endif
        REAL(wp)    ::  w_comp !< advection velocity along z
        
+       REAL(wp)    :: diss_s  !< discretized artificial dissipation at southward-side of the grid box
+       REAL(wp)    :: flux_s  !< discretized 6th-order flux at southward-side of the grid box
+#ifndef _OPENACC
        REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_diss_y_local_u !< discretized artificial dissipation at southward-side of the grid box
        REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_flux_y_local_u !< discretized 6th-order flux at southward-side of the grid box
+#endif
        
+       REAL(wp)    :: diss_l  !< discretized artificial dissipation at leftward-side of the grid box
+       REAL(wp)    :: flux_l  !< discretized 6th-order flux at leftward-side of the grid box
+#ifndef _OPENACC
        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_diss_x_local_u !< discretized artificial dissipation at leftward-side of the grid box
        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_flux_x_local_u !< discretized 6th-order flux at leftward-side of the grid box
+#endif
        
-       REAL(wp), DIMENSION(nzb:nzt) ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt) ::  diss_r !< discretized artificial dissipation at leftward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt) ::  diss_t !< discretized artificial dissipation at top of the grid box
-       REAL(wp), DIMENSION(nzb:nzt) ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt) ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt) ::  flux_t !< discretized 6th-order flux at top of the grid box
-       REAL(wp), DIMENSION(nzb:nzt) ::  u_comp !< advection velocity along x
+       REAL(wp) ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
+       REAL(wp) ::  diss_r !< discretized artificial dissipation at leftward-side of the grid box
+       REAL(wp) ::  diss_t !< discretized artificial dissipation at top of the grid box
+       REAL(wp) ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
+       REAL(wp) ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
+       REAL(wp) ::  flux_t !< discretized 6th-order flux at top of the grid box
+       REAL(wp) ::  u_comp !< advection velocity along x
+#ifdef _OPENACC
+       REAL(wp)    ::  u_comp_l !<
+#endif
  
        gu = 2.0_wp * u_gtrans
        gv = 2.0_wp * v_gtrans
 
+#ifndef _OPENACC
 !
 !--    Compute the fluxes for the whole left boundary of the processor domain.
@@ -3915,6 +4232,6 @@
              ibit9  = REAL( IBITS(advc_flags_1(k,j,i-1),9,1),  KIND = wp )
 
-             u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
-             swap_flux_x_local_u(k,j) = u_comp(k) * (                          &
+             u_comp                   = u(k,j,i) + u(k,j,i-1) - gu
+             swap_flux_x_local_u(k,j) = u_comp * (                             &
                                        ( 37.0_wp * ibit11 * adv_mom_5          &
                                     +     7.0_wp * ibit10 * adv_mom_3          &
@@ -3931,5 +4248,5 @@
                                                    )
 
-              swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * (                &
+              swap_diss_x_local_u(k,j) = - ABS( u_comp ) * (                   &
                                        ( 10.0_wp * ibit11 * adv_mom_5          &
                                     +     3.0_wp * ibit10 * adv_mom_3          &
@@ -3950,10 +4267,10 @@
           DO  k = nzb_max+1, nzt
 
-             u_comp(k)         = u(k,j,i) + u(k,j,i-1) - gu
-             swap_flux_x_local_u(k,j) = u_comp(k) * (                          &
+             u_comp            = u(k,j,i) + u(k,j,i-1) - gu
+             swap_flux_x_local_u(k,j) = u_comp * (                             &
                              37.0_wp * ( u(k,j,i) + u(k,j,i-1)   )             &
                            -  8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) )             &
                            +           ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
-             swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * (                   &
+             swap_diss_x_local_u(k,j) = - ABS(u_comp) * (                      &
                              10.0_wp * ( u(k,j,i) - u(k,j,i-1)   )             &
                            -  5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) )             &
@@ -3962,6 +4279,25 @@
           ENDDO
        ENDDO
-
+#endif
+
+       !$ACC PARALLEL LOOP COLLAPSE(2) FIRSTPRIVATE(tn, gu, gv) &
+       !$ACC PRIVATE(i, j, k, k_mm, k_pp, k_ppp) &
+       !$ACC PRIVATE(ibit9, ibit10, ibit11, ibit12, ibit13, ibit14) &
+       !$ACC PRIVATE(ibit9_l, ibit10_l, ibit11_l) &
+       !$ACC PRIVATE(ibit12_s, ibit13_s, ibit14_s) &
+       !$ACC PRIVATE(ibit15, ibit16, ibit17) &
+       !$ACC PRIVATE(flux_r, diss_r, flux_l, diss_l) &
+       !$ACC PRIVATE(flux_n, diss_n, flux_s, diss_s) &
+       !$ACC PRIVATE(flux_t, diss_t, flux_d, diss_d) &
+       !$ACC PRIVATE(div, u_comp, u_comp_l, v_comp, v_comp_s, w_comp) &
+       !$ACC PRESENT(advc_flags_1) &
+       !$ACC PRESENT(u, v, w) &
+       !$ACC PRESENT(drho_air, rho_air_zw, ddzw) &
+       !$ACC PRESENT(tend) &
+       !$ACC PRESENT(hom(nzb+1:nzb_max,1,1:3,0)) &
+       !$ACC PRESENT(weight_substep(intermediate_timestep_count)) &
+       !$ACC PRESENT(sums_us2_ws_l, sums_wsus_ws_l)
        DO i = nxlu, nxr
+#ifndef _OPENACC
 !       
 !--       The following loop computes the fluxes for the south boundary points
@@ -4019,10 +4355,10 @@
 
           ENDDO
+#endif
+
 !
 !--       Computation of interior fluxes and tendency terms
           DO  j = nys, nyn
 
-             flux_t(0) = 0.0_wp
-             diss_t(0) = 0.0_wp
              flux_d    = 0.0_wp
              diss_d    = 0.0_wp
@@ -4034,6 +4370,6 @@
                 ibit9  = REAL( IBITS(advc_flags_1(k,j,i),9,1),  KIND = wp )
 
-                u_comp(k) = u(k,j,i+1) + u(k,j,i)
-                flux_r(k) = ( u_comp(k) - gu ) * (                           &
+                u_comp = u(k,j,i+1) + u(k,j,i)
+                flux_r = ( u_comp - gu ) * (                                 &
                           ( 37.0_wp * ibit11 * adv_mom_5                     &
                        +     7.0_wp * ibit10 * adv_mom_3                     &
@@ -4050,5 +4386,5 @@
                                                  )
 
-                diss_r(k) = - ABS( u_comp(k) - gu ) * (                      &
+                diss_r = - ABS( u_comp - gu ) * (                            &
                           ( 10.0_wp * ibit11 * adv_mom_5                     &
                        +     3.0_wp * ibit10 * adv_mom_3                     & 
@@ -4065,10 +4401,52 @@
                                                      )
 
+#ifdef _OPENACC
+!
+!--             Recompute the left fluxes.
+                ibit11_l = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )
+                ibit10_l = REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )
+                ibit9_l  = REAL( IBITS(advc_flags_1(k,j,i-1),9,1),  KIND = wp )
+
+                u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
+                flux_l   = u_comp_l * (                                        &
+                                       ( 37.0_wp * ibit11_l * adv_mom_5           &
+                                    +     7.0_wp * ibit10_l * adv_mom_3           &
+                                    +              ibit9_l  * adv_mom_1           &
+                                       ) *                                     &
+                                     ( u(k,j,i)   + u(k,j,i-1) )               &
+                                -      (  8.0_wp * ibit11_l * adv_mom_5           &
+                                    +              ibit10_l * adv_mom_3           &
+                                       ) *                                     &
+                                     ( u(k,j,i+1) + u(k,j,i-2) )               &
+                                +      (           ibit11_l * adv_mom_5           &
+                                       ) *                                     &
+                                     ( u(k,j,i+2) + u(k,j,i-3) )               &
+                                                   )
+
+                diss_l   = - ABS( u_comp_l ) * (                               &
+                                       ( 10.0_wp * ibit11_l * adv_mom_5           &
+                                    +     3.0_wp * ibit10_l * adv_mom_3           &
+                                    +              ibit9_l  * adv_mom_1           &
+                                       ) *                                     &
+                                     ( u(k,j,i)   - u(k,j,i-1) )               &
+                                -      (  5.0_wp * ibit11_l * adv_mom_5           &
+                                    +              ibit10_l * adv_mom_3           &
+                                       ) *                                     &
+                                     ( u(k,j,i+1) - u(k,j,i-2) )               &
+                                +      (           ibit11_l * adv_mom_5           &
+                                       ) *                                     &
+                                     ( u(k,j,i+2) - u(k,j,i-3) )               &
+                                                             )
+#else
+                flux_l = swap_flux_x_local_u(k,j)
+                diss_l = swap_diss_x_local_u(k,j)
+#endif
+
                 ibit14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp )
                 ibit13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp )
                 ibit12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp )
 
-                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
-                flux_n(k) = v_comp * (                                       &
+                v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
+                flux_n = v_comp * (                                          &
                           ( 37.0_wp * ibit14 * adv_mom_5                     &
                        +     7.0_wp * ibit13 * adv_mom_3                     &
@@ -4085,5 +4463,5 @@
                                                  )
 
-                diss_n(k) = - ABS ( v_comp ) * (                             &
+                diss_n = - ABS ( v_comp ) * (                                &
                           ( 10.0_wp * ibit14 * adv_mom_5                     &
                        +     3.0_wp * ibit13 * adv_mom_3                     &
@@ -4099,4 +4477,47 @@
                                  ( u(k,j+3,i) - u(k,j-2,i) )                 &
                                                       )
+
+#ifdef _OPENACC
+!
+!--             Recompute the south fluxes.
+                ibit14_s = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )
+                ibit13_s = REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp )
+                ibit12_s = REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp )
+
+                v_comp_s = v(k,j,i) + v(k,j,i-1) - gv
+                flux_s   = v_comp_s * (                                      &
+                                   ( 37.0_wp * ibit14_s * adv_mom_5             &
+                                +     7.0_wp * ibit13_s * adv_mom_3             &
+                                +              ibit12_s * adv_mom_1             &
+                                   ) *                                       &
+                                     ( u(k,j,i)   + u(k,j-1,i) )             &
+                            -      (  8.0_wp * ibit14_s * adv_mom_5             &
+                            +                  ibit13_s * adv_mom_3             &
+                                   ) *                                       &
+                                     ( u(k,j+1,i) + u(k,j-2,i) )             &
+                        +      (               ibit14_s * adv_mom_5             &
+                               ) *                                           &
+                                     ( u(k,j+2,i) + u(k,j-3,i) )             &
+                                               )
+
+                diss_s   = - ABS ( v_comp_s ) * (                            &
+                                   ( 10.0_wp * ibit14_s * adv_mom_5              &
+                                +     3.0_wp * ibit13_s * adv_mom_3              &
+                                +              ibit12_s * adv_mom_1              &
+                                   ) *                                        &
+                                     ( u(k,j,i)   - u(k,j-1,i) )              &
+                            -      (  5.0_wp * ibit14_s * adv_mom_5              &
+                                +              ibit13_s * adv_mom_3              &
+                                   ) *                                        &
+                                     ( u(k,j+1,i) - u(k,j-2,i) )              &
+                            +      (           ibit14_s * adv_mom_5              &
+                                   ) *                                        &
+                                     ( u(k,j+2,i) - u(k,j-3,i) )              &
+                                                         )
+#else
+                flux_s = swap_flux_y_local_u(k)
+                diss_s = swap_diss_y_local_u(k)
+#endif
+
 !
 !--             k index has to be modified near bottom and top, else array
@@ -4110,6 +4531,6 @@
                 k_mm  = k - 2 * ibit17
 
-                w_comp    = w(k,j,i) + w(k,j,i-1)
-                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
+                w_comp = w(k,j,i) + w(k,j,i-1)
+                flux_t = w_comp * rho_air_zw(k) * (                          &
                           ( 37.0_wp * ibit17 * adv_mom_5                     &
                        +     7.0_wp * ibit16 * adv_mom_3                     &
@@ -4126,5 +4547,5 @@
                                       )
 
-                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
+                diss_t = - ABS( w_comp ) * rho_air_zw(k) * (                 &
                           ( 10.0_wp * ibit17 * adv_mom_5                     &
                        +     3.0_wp * ibit16 * adv_mom_3                     &
@@ -4144,5 +4565,5 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 )             &
+                div = ( ( u_comp * ( ibit9 + ibit10 + ibit11 )                &
                 - ( u(k,j,i)   + u(k,j,i-1)   )                               &
                                     * (                                       &
@@ -4173,37 +4594,40 @@
 
                 tend(k,j,i) = tend(k,j,i) - (                                  &
-                 ( flux_r(k) + diss_r(k)                                       &
-               -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
-               + ( flux_n(k) + diss_n(k)                                       &
-               -   swap_flux_y_local_u(k)   - swap_diss_y_local_u(k)   ) * ddy &
-               + ( ( flux_t(k) + diss_t(k) )                                   &
-               -   ( flux_d    + diss_d    )                                   &
-                                                    ) * drho_air(k) * ddzw(k)  &
+                 ( ( flux_r + diss_r )                                         &
+               -   ( flux_l + diss_l ) ) * ddx                                 &
+               + ( ( flux_n + diss_n )                                         &
+               -   ( flux_s + diss_s ) ) * ddy                                 &
+               + ( ( flux_t + diss_t )                                         &
+               -   ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k)               &
                                            ) + div * u(k,j,i)
 
-                swap_flux_x_local_u(k,j) = flux_r(k)
-                swap_diss_x_local_u(k,j) = diss_r(k)
-                swap_flux_y_local_u(k)   = flux_n(k)
-                swap_diss_y_local_u(k)   = diss_n(k)
-                flux_d                   = flux_t(k)
-                diss_d                   = diss_t(k)
+#ifndef _OPENACC
+                swap_flux_x_local_u(k,j) = flux_r
+                swap_diss_x_local_u(k,j) = diss_r
+                swap_flux_y_local_u(k)   = flux_n
+                swap_diss_y_local_u(k)   = diss_n
+#endif
+                flux_d                   = flux_t
+                diss_d                   = diss_t
 !
 !--             Statistical Evaluation of u'u'. The factor has to be applied
 !--             for right evaluation when gallilei_trans = .T. .
+                !$ACC ATOMIC
                 sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                      &
-                + ( flux_r(k)                                                  &
-                    * ( u_comp(k) - 2.0_wp * hom(k,1,1,0)                   )  &
-                    / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )  &
-                  + diss_r(k)                                                  &
-                    *   ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0)              )  &
-                    / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp                  )  &
+                + ( flux_r                                                     &
+                    * ( u_comp - 2.0_wp * hom(k,1,1,0)                   )     &
+                    / ( u_comp - gu + SIGN( 1.0E-20_wp, u_comp - gu )    )     &
+                  + diss_r                                                     &
+                    *   ABS( u_comp - 2.0_wp * hom(k,1,1,0)              )     &
+                    / ( ABS( u_comp - gu ) + 1.0E-20_wp                  )     &
                   ) *   weight_substep(intermediate_timestep_count)
 !
 !--             Statistical Evaluation of w'u'.
+                !$ACC ATOMIC
                 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                    &
-                + ( flux_t(k)                                                  &
+                + ( flux_t                                                     &
                     * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
                     / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
-                  + diss_t(k)                                                  &
+                  + diss_t                                                     &
                     *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
                     / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
@@ -4214,23 +4638,58 @@
              DO  k = nzb_max+1, nzt
 
-                u_comp(k) = u(k,j,i+1) + u(k,j,i)
-                flux_r(k) = ( u_comp(k) - gu ) * (                            &
+                u_comp = u(k,j,i+1) + u(k,j,i)
+                flux_r = ( u_comp - gu ) * (                                  &
                          37.0_wp * ( u(k,j,i+1) + u(k,j,i)   )                   &
                        -  8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) )                   &
                        +           ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
-                diss_r(k) = - ABS( u_comp(k) - gu ) * (                       &
+                diss_r = - ABS( u_comp - gu ) * (                             &
                          10.0_wp * ( u(k,j,i+1) - u(k,j,i)   )                   &
                        -  5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) )                   &
                        +           ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
 
-                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
-                flux_n(k) = v_comp * (                                        &
+#ifdef _OPENACC
+!
+!--             Recompute the left fluxes.
+                u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
+                flux_l   = u_comp_l * (                                       &
+                             37.0_wp * ( u(k,j,i) + u(k,j,i-1)   )                &
+                           -  8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) )                &
+                           +           ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
+                diss_l   = - ABS(u_comp_l) * (                                &
+                             10.0_wp * ( u(k,j,i) - u(k,j,i-1)   )                &
+                           -  5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) )                &
+                           +           ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
+#else
+                flux_l = swap_flux_x_local_u(k,j)
+                diss_l = swap_diss_x_local_u(k,j)
+#endif
+
+                v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
+                flux_n = v_comp * (                                           &
                          37.0_wp * ( u(k,j+1,i) + u(k,j,i)   )                   &
                        -  8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) )                   &
                        +           ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
-                diss_n(k) = - ABS( v_comp ) * (                               &
+                diss_n = - ABS( v_comp ) * (                                  &
                          10.0_wp * ( u(k,j+1,i) - u(k,j,i)   )                   &
                        -  5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) )                   &
                        +           ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
+
+#ifdef _OPENACC
+!
+!--             Recompute the south fluxes.
+                v_comp_s = v(k,j,i) + v(k,j,i-1) - gv
+                flux_s   = v_comp_s * (                                       &
+                           37.0_wp * ( u(k,j,i) + u(k,j-1,i)   )                 &
+                         -  8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) )                 &
+                         +           ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
+                diss_s   = - ABS( v_comp_s ) * (                              &
+                           10.0_wp * ( u(k,j,i) - u(k,j-1,i)   )                 &
+                         -  5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) )                 &
+                         +           ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
+#else
+                flux_s = swap_flux_y_local_u(k)
+                diss_s = swap_diss_y_local_u(k)
+#endif
+
 !
 !--             k index has to be modified near bottom and top, else array
@@ -4244,6 +4703,6 @@
                 k_mm  = k - 2 * ibit17
 
-                w_comp    = w(k,j,i) + w(k,j,i-1)
-                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
+                w_comp = w(k,j,i) + w(k,j,i-1)
+                flux_t = w_comp * rho_air_zw(k) * (                          &
                           ( 37.0_wp * ibit17 * adv_mom_5                        &
                        +     7.0_wp * ibit16 * adv_mom_3                        &
@@ -4260,5 +4719,5 @@
                                       )
 
-                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
+                diss_t = - ABS( w_comp ) * rho_air_zw(k) * (                 &
                           ( 10.0_wp * ibit17 * adv_mom_5                        &
                        +     3.0_wp * ibit16 * adv_mom_3                        &
@@ -4278,5 +4737,5 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx &
+                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                      * rho_air_zw(k) -     &
@@ -4286,37 +4745,40 @@
 
                 tend(k,j,i) = tend(k,j,i) - (                                  &
-                 ( flux_r(k) + diss_r(k)                                       &
-               -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
-               + ( flux_n(k) + diss_n(k)                                       &
-               -   swap_flux_y_local_u(k)   - swap_diss_y_local_u(k)   ) * ddy &
-               + ( ( flux_t(k) + diss_t(k) )                                   &
-               -   ( flux_d    + diss_d    )                                   &
-                                                    ) * drho_air(k) * ddzw(k)  &
+                 ( ( flux_r + diss_r )                                         &
+               -   ( flux_l + diss_l ) ) * ddx                                 &
+               + ( ( flux_n + diss_n )                                         &
+               -   ( flux_s + diss_s ) ) * ddy                                 &
+               + ( ( flux_t + diss_t )                                         &
+               -   ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k)               &
                                            ) + div * u(k,j,i)
 
-                swap_flux_x_local_u(k,j) = flux_r(k)
-                swap_diss_x_local_u(k,j) = diss_r(k)
-                swap_flux_y_local_u(k)   = flux_n(k)
-                swap_diss_y_local_u(k)   = diss_n(k)
-                flux_d                   = flux_t(k)
-                diss_d                   = diss_t(k)
+#ifndef _OPENACC
+                swap_flux_x_local_u(k,j) = flux_r
+                swap_diss_x_local_u(k,j) = diss_r
+                swap_flux_y_local_u(k)   = flux_n
+                swap_diss_y_local_u(k)   = diss_n
+#endif
+                flux_d                   = flux_t
+                diss_d                   = diss_t
 !
 !--             Statistical Evaluation of u'u'. The factor has to be applied
 !--             for right evaluation when gallilei_trans = .T. .
+                !$ACC ATOMIC
                 sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                      &
-                + ( flux_r(k)                                                  &
-                    * ( u_comp(k) - 2.0_wp * hom(k,1,1,0)                   )  &
-                    / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )  &
-                  + diss_r(k)                                                  &
-                    *   ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0)              )  &
-                    / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp                  )  &
+                + ( flux_r                                                     &
+                    * ( u_comp - 2.0_wp * hom(k,1,1,0)                   )     &
+                    / ( u_comp - gu + SIGN( 1.0E-20_wp, u_comp - gu )    )     &
+                  + diss_r                                                     &
+                    *   ABS( u_comp - 2.0_wp * hom(k,1,1,0)              )     &
+                    / ( ABS( u_comp - gu ) + 1.0E-20_wp                  )     &
                   ) *   weight_substep(intermediate_timestep_count)
 !
 !--             Statistical Evaluation of w'u'.
+                !$ACC ATOMIC
                 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                    &
-                + ( flux_t(k)                                                  &
+                + ( flux_t                                                     &
                     * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
                     / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
-                  + diss_t(k)                                                  &
+                  + diss_t                                                     &
                     *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
                     / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
@@ -4369,7 +4831,17 @@
        REAL(wp)    ::  ibit19 !< flag indicating 3rd-order scheme along x-direction
        REAL(wp)    ::  ibit20 !< flag indicating 5th-order scheme along x-direction
+#ifdef _OPENACC
+       REAL(wp)    ::  ibit18_l !< flag indicating 1st-order scheme along x-direction
+       REAL(wp)    ::  ibit19_l !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp)    ::  ibit20_l !< flag indicating 5th-order scheme along x-direction
+#endif
        REAL(wp)    ::  ibit21 !< flag indicating 1st-order scheme along y-direction
        REAL(wp)    ::  ibit22 !< flag indicating 3rd-order scheme along y-direction
        REAL(wp)    ::  ibit23 !< flag indicating 5th-order scheme along y-direction
+#ifdef _OPENACC
+       REAL(wp)    ::  ibit21_s !< flag indicating 1st-order scheme along y-direction
+       REAL(wp)    ::  ibit22_s !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp)    ::  ibit23_s !< flag indicating 5th-order scheme along y-direction
+#endif
        REAL(wp)    ::  ibit24 !< flag indicating 1st-order scheme along z-direction
        REAL(wp)    ::  ibit25 !< flag indicating 3rd-order scheme along z-direction
@@ -4381,22 +4853,38 @@
        REAL(wp)    ::  gv     !< Galilei-transformation velocity along y
        REAL(wp)    ::  u_comp !< advection velocity along x
+#ifdef _OPENACC
+       REAL(wp)    ::  u_comp_l !< advection velocity along x
+#endif
        REAL(wp)    ::  w_comp !< advection velocity along z
        
+       REAL(wp)    ::  diss_s !< discretized artificial dissipation at southward-side of the grid box
+       REAL(wp)    ::  flux_s !< discretized 6th-order flux at southward-side of the grid box
+#ifndef _OPENACC
        REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_diss_y_local_v !< discretized artificial dissipation at southward-side of the grid box
        REAL(wp), DIMENSION(nzb+1:nzt) ::  swap_flux_y_local_v !< discretized 6th-order flux at southward-side of the grid box
+#endif
        
+       REAL(wp)    ::  diss_l !< discretized artificial dissipation at leftward-side of the grid box
+       REAL(wp)    ::  flux_l !< discretized 6th-order flux at leftward-side of the grid box
+#ifndef _OPENACC
        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_diss_x_local_v !< discretized artificial dissipation at leftward-side of the grid box
        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_flux_x_local_v !< discretized 6th-order flux at leftward-side of the grid box
+#endif
        
-       REAL(wp), DIMENSION(nzb:nzt) ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt) ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt) ::  diss_t !< discretized artificial dissipation at top of the grid box
-       REAL(wp), DIMENSION(nzb:nzt) ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt) ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
-       REAL(wp), DIMENSION(nzb:nzt) ::  flux_t !< discretized 6th-order flux at top of the grid box
-       REAL(wp), DIMENSION(nzb:nzt) ::  v_comp !< advection velocity along y
+       REAL(wp)    ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
+       REAL(wp)    ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
+       REAL(wp)    ::  diss_t !< discretized artificial dissipation at top of the grid box
+       REAL(wp)    ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
+       REAL(wp)    ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
+       REAL(wp)    ::  flux_t !< discretized 6th-order flux at top of the grid box
+       REAL(wp)    ::  v_comp !< advection velocity along y
+#ifdef _OPENACC
+       REAL(wp)    ::  v_comp_s !<
+#endif
 
        gu = 2.0_wp * u_gtrans
        gv = 2.0_wp * v_gtrans
+
+#ifndef _OPENACC
 !
 !--    First compute the whole left boundary of the processor domain
@@ -4457,7 +4945,26 @@
 
        ENDDO
-
+#endif
+
+       !$ACC PARALLEL LOOP COLLAPSE(2) FIRSTPRIVATE(tn, gu, gv) &
+       !$ACC PRIVATE(i, j, k, k_mm, k_pp, k_ppp) &
+       !$ACC PRIVATE(ibit18, ibit19, ibit20, ibit21, ibit22, ibit23) &
+       !$ACC PRIVATE(ibit18_l, ibit19_l, ibit20_l) &
+       !$ACC PRIVATE(ibit21_s, ibit22_s, ibit23_s) &
+       !$ACC PRIVATE(ibit24, ibit25, ibit26) &
+       !$ACC PRIVATE(flux_r, diss_r, flux_l, diss_l) &
+       !$ACC PRIVATE(flux_n, diss_n, flux_s, diss_s) &
+       !$ACC PRIVATE(flux_t, diss_t, flux_d, diss_d) &
+       !$ACC PRIVATE(div, u_comp, u_comp_l, v_comp, v_comp_s, w_comp) &
+       !$ACC PRESENT(advc_flags_1) &
+       !$ACC PRESENT(u, v, w) &
+       !$ACC PRESENT(drho_air, rho_air_zw, ddzw) &
+       !$ACC PRESENT(tend) &
+       !$ACC PRESENT(hom(nzb+1:nzb_max,1,2:3,0)) &
+       !$ACC PRESENT(weight_substep(intermediate_timestep_count)) &
+       !$ACC PRESENT(sums_vs2_ws_l, sums_wsvs_ws_l)
        DO i = nxl, nxr
 
+#ifndef _OPENACC
           j = nysv
           DO  k = nzb+1, nzb_max
@@ -4467,6 +4974,6 @@
              ibit21 = REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )
 
-             v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
-             swap_flux_y_local_v(k) = v_comp(k) * (                           &
+             v_comp                 = v(k,j,i) + v(k,j-1,i) - gv
+             swap_flux_y_local_v(k) = v_comp * (                              &
                                    ( 37.0_wp * ibit23 * adv_mom_5                &
                                 +     7.0_wp * ibit22 * adv_mom_3                &
@@ -4483,5 +4990,5 @@
                                                  )
 
-             swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * (                  &
+             swap_diss_y_local_v(k) = - ABS( v_comp ) * (                     &
                                    ( 10.0_wp * ibit23 * adv_mom_5                &
                                 +     3.0_wp * ibit22 * adv_mom_3                &
@@ -4502,10 +5009,10 @@
           DO  k = nzb_max+1, nzt
 
-             v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
-             swap_flux_y_local_v(k) = v_comp(k) * (                           &
+             v_comp                 = v(k,j,i) + v(k,j-1,i) - gv
+             swap_flux_y_local_v(k) = v_comp * (                              &
                            37.0_wp * ( v(k,j,i) + v(k,j-1,i)   )                 &
                          -  8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) )                 &
                          +           ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
-             swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * (                  &
+             swap_diss_y_local_v(k) = - ABS( v_comp ) * (                     &
                            10.0_wp * ( v(k,j,i) - v(k,j-1,i)   )                 &
                          -  5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) )                 &
@@ -4513,9 +5020,8 @@
 
           ENDDO
+#endif
 
           DO  j = nysv, nyn
 
-             flux_t(0) = 0.0_wp
-             diss_t(0) = 0.0_wp
              flux_d    = 0.0_wp
              diss_d    = 0.0_wp
@@ -4527,6 +5033,6 @@
                 ibit18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp )
 
-                u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
-                flux_r(k) = u_comp * (                                       &
+                u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
+                flux_r = u_comp * (                                          &
                           ( 37.0_wp * ibit20 * adv_mom_5                        &
                        +     7.0_wp * ibit19 * adv_mom_3                        &
@@ -4543,5 +5049,5 @@
                                      )
 
-                diss_r(k) = - ABS( u_comp ) * (                              &
+                diss_r = - ABS( u_comp ) * (                                 &
                           ( 10.0_wp * ibit20 * adv_mom_5                        &
                        +     3.0_wp * ibit19 * adv_mom_3                        &
@@ -4558,10 +5064,52 @@
                                               )
 
+#ifdef _OPENACC
+!
+!--             Recompute the left fluxes.
+                ibit20_l = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )
+                ibit19_l = REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp )
+                ibit18_l = REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp )
+
+                u_comp_l  = u(k,j-1,i) + u(k,j,i) - gu
+                flux_l    = u_comp_l * (                                     &
+                                      ( 37.0_wp * ibit20_l * adv_mom_5          &
+                                   +     7.0_wp * ibit19_l * adv_mom_3          &
+                                   +              ibit18_l * adv_mom_1          &
+                                      ) *                                    &
+                                     ( v(k,j,i)   + v(k,j,i-1) )             &
+                               -      (  8.0_wp * ibit20_l * adv_mom_5          &
+                                   +              ibit19_l * adv_mom_3          &
+                                      ) *                                    &
+                                     ( v(k,j,i+1) + v(k,j,i-2) )             &
+                               +      (           ibit20_l * adv_mom_5          &
+                                      ) *                                    &
+                                     ( v(k,j,i+2) + v(k,j,i-3) )             &
+                                                 )
+
+                 diss_l   = - ABS( u_comp_l ) * (                            &
+                                      ( 10.0_wp * ibit20_l * adv_mom_5          &
+                                   +     3.0_wp * ibit19_l * adv_mom_3          &
+                                   +              ibit18_l * adv_mom_1          &
+                                      ) *                                    &
+                                     ( v(k,j,i)   - v(k,j,i-1) )             &
+                               -      (  5.0_wp * ibit20_l * adv_mom_5          &
+                                   +              ibit19_l * adv_mom_3          &
+                                      ) *                                    &
+                                     ( v(k,j,i+1) - v(k,j,i-2) )             &
+                               +      (           ibit20_l * adv_mom_5          &
+                                      ) *                                    &
+                                     ( v(k,j,i+2) - v(k,j,i-3) )             &
+                                                           )
+#else
+                flux_l = swap_flux_x_local_v(k,j)
+                diss_l = swap_diss_x_local_v(k,j)
+#endif
+
                 ibit23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp )
                 ibit22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp )
                 ibit21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp )
 
-                v_comp(k) = v(k,j+1,i) + v(k,j,i)
-                flux_n(k) = ( v_comp(k) - gv ) * (                           &
+                v_comp = v(k,j+1,i) + v(k,j,i)
+                flux_n = ( v_comp - gv ) * (                                 &
                           ( 37.0_wp * ibit23 * adv_mom_5                        &
                        +     7.0_wp * ibit22 * adv_mom_3                        &
@@ -4578,5 +5126,5 @@
                                      )
 
-                diss_n(k) = - ABS( v_comp(k) - gv ) * (                      &
+                diss_n = - ABS( v_comp - gv ) * (                            &
                           ( 10.0_wp * ibit23 * adv_mom_5                        &
                        +     3.0_wp * ibit22 * adv_mom_3                        &
@@ -4592,4 +5140,47 @@
                                  ( v(k,j+3,i) - v(k,j-2,i) )                 &
                                                       )
+
+#ifdef _OPENACC
+!
+!--             Recompute the south fluxes.
+                ibit23_s = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )
+                ibit22_s = REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp )
+                ibit21_s = REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )
+
+                v_comp_s = v(k,j,i) + v(k,j-1,i) - gv
+                flux_s   = v_comp_s * (                                      &
+                                   ( 37.0_wp * ibit23_s * adv_mom_5             &
+                                +     7.0_wp * ibit22_s * adv_mom_3             &
+                                +              ibit21_s * adv_mom_1             &
+                                   ) *                                       &
+                                     ( v(k,j,i)   + v(k,j-1,i) )             &
+                            -      (  8.0_wp * ibit23_s * adv_mom_5             &
+                                +              ibit22_s * adv_mom_3             &
+                                   ) *                                       &
+                                     ( v(k,j+1,i) + v(k,j-2,i) )             &
+                            +      (           ibit23_s * adv_mom_5             &
+                                   ) *                                       &
+                                     ( v(k,j+2,i) + v(k,j-3,i) )             &
+                                                 )
+
+                diss_s   = - ABS( v_comp_s ) * (                             &
+                                   ( 10.0_wp * ibit23_s * adv_mom_5             &
+                                +     3.0_wp * ibit22_s * adv_mom_3             &
+                                +              ibit21_s * adv_mom_1             &
+                                   ) *                                       &
+                                     ( v(k,j,i)   - v(k,j-1,i) )             &
+                            -      (  5.0_wp * ibit23_s * adv_mom_5             &
+                                +              ibit22_s * adv_mom_3             &
+                                   ) *                                       &
+                                     ( v(k,j+1,i) - v(k,j-2,i) )             &
+                            +      (           ibit23_s * adv_mom_5             &
+                                   ) *                                       &
+                                     ( v(k,j+2,i) - v(k,j-3,i) )             &
+                                                          )
+#else
+               flux_s = swap_flux_y_local_v(k)
+               diss_s = swap_diss_y_local_v(k)
+#endif
+
 !
 !--             k index has to be modified near bottom and top, else array
@@ -4603,6 +5194,6 @@
                 k_mm  = k - 2 * ibit26
 
-                w_comp    = w(k,j-1,i) + w(k,j,i)
-                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
+                w_comp = w(k,j-1,i) + w(k,j,i)
+                flux_t = w_comp * rho_air_zw(k) * (                          &
                           ( 37.0_wp * ibit26 * adv_mom_5                        &
                        +     7.0_wp * ibit25 * adv_mom_3                        &
@@ -4619,5 +5210,5 @@
                                       )
 
-                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
+                diss_t = - ABS( w_comp ) * rho_air_zw(k) * (                 &
                           ( 10.0_wp * ibit26 * adv_mom_5                        &
                        +     3.0_wp * ibit25 * adv_mom_3                        &
@@ -4646,5 +5237,5 @@
                                          )                                    &
                   ) * ddx                                                     &
-               +  ( v_comp(k)                                                 &
+               +  ( v_comp                                                    &
                                        * ( ibit21 + ibit22 + ibit23 )         &
                 - ( v(k,j,i)     + v(k,j-1,i) )                               &
@@ -4668,40 +5259,41 @@
 
                 tend(k,j,i) = tend(k,j,i) - (                                 &
-                       ( flux_r(k) + diss_r(k)                                &
-                     -   swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j)  &
-                       ) * ddx                                                &
-                     + ( flux_n(k) + diss_n(k)                                &
-                     -   swap_flux_y_local_v(k) - swap_diss_y_local_v(k)      &
-                       ) * ddy                                                &
-                     + ( ( flux_t(k) + diss_t(k) )                            &
-                     -   ( flux_d    + diss_d    )                            &
-                       ) * drho_air(k) * ddzw(k)                              &
+                       ( ( flux_r + diss_r )                                  &
+                     -   ( flux_l + diss_l ) ) * ddx                          &
+                     + ( ( flux_n + diss_n )                                  &
+                     -   ( flux_s + diss_s ) ) * ddy                          &
+                     + ( ( flux_t + diss_t )                                  &
+                     -   ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k)        &
                                             )  + v(k,j,i) * div
 
-                swap_flux_x_local_v(k,j) = flux_r(k)
-                swap_diss_x_local_v(k,j) = diss_r(k)
-                swap_flux_y_local_v(k)   = flux_n(k)
-                swap_diss_y_local_v(k)   = diss_n(k)
-                flux_d                   = flux_t(k)
-                diss_d                   = diss_t(k)
+#ifndef _OPENACC
+                swap_flux_x_local_v(k,j) = flux_r
+                swap_diss_x_local_v(k,j) = diss_r
+                swap_flux_y_local_v(k)   = flux_n
+                swap_diss_y_local_v(k)   = diss_n
+#endif
+                flux_d                   = flux_t
+                diss_d                   = diss_t
 
 !
 !--             Statistical Evaluation of v'v'. The factor has to be applied
 !--             for right evaluation when gallilei_trans = .T. .
+                !$ACC ATOMIC
                 sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                      &
-                + ( flux_n(k)                                                  &
-                    * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)                   )  &
-                    / ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) )  &
-               +   diss_n(k)                                                   &
-                    *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0)              )  &
-                    / ( ABS( v_comp(k) - gv ) + 1.0E-20_wp                  )  &
+                + ( flux_n                                                     &
+                    * ( v_comp - 2.0_wp * hom(k,1,2,0)                   )     &
+                    / ( v_comp - gv + SIGN( 1.0E-20_wp, v_comp - gv )    )     &
+               +   diss_n                                                      &
+                    *   ABS( v_comp - 2.0_wp * hom(k,1,2,0)              )     &
+                    / ( ABS( v_comp - gv ) + 1.0E-20_wp                  )     &
                   ) *   weight_substep(intermediate_timestep_count)
 !
 !--             Statistical Evaluation of w'u'.
+                !$ACC ATOMIC
                 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                    &
-                + ( flux_t(k)                                                  &
+                + ( flux_t                                                     &
                     * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
                     / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
-               +   diss_t(k)                                                   &
+               +   diss_t                                                      &
                     *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
                     / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
@@ -4712,26 +5304,60 @@
              DO  k = nzb_max+1, nzt
 
-                u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
-                flux_r(k) = u_comp * (                                        &
+                u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
+                flux_r = u_comp * (                                           &
                       37.0_wp * ( v(k,j,i+1) + v(k,j,i)   )                      &
                     -  8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) )                      &
                     +           ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
 
-                diss_r(k) = - ABS( u_comp ) * (                               &
+                diss_r = - ABS( u_comp ) * (                                  &
                       10.0_wp * ( v(k,j,i+1) - v(k,j,i) )                        &
                     -  5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) )                      &
                     +           ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
 
-
-                v_comp(k) = v(k,j+1,i) + v(k,j,i)
-                flux_n(k) = ( v_comp(k) - gv ) * (                            &
+#ifdef _OPENACC
+!
+!--             Recompute the left fluxes.
+                u_comp_l = u(k,j-1,i) + u(k,j,i) - gu
+                flux_l   = u_comp_l * (                                       &
+                             37.0_wp * ( v(k,j,i) + v(k,j,i-1)   )               &
+                           -  8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) )               &
+                           +           ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
+                diss_l   = - ABS( u_comp_l ) * (                              &
+                             10.0_wp * ( v(k,j,i) - v(k,j,i-1)   )               &
+                           -  5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) )               &
+                           +           ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
+#else
+                flux_l = swap_flux_x_local_v(k,j)
+                diss_l = swap_diss_x_local_v(k,j)
+#endif
+
+                v_comp = v(k,j+1,i) + v(k,j,i)
+                flux_n = ( v_comp - gv ) * (                                  &
                       37.0_wp * ( v(k,j+1,i) + v(k,j,i)   )                      &
                     -  8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) )                      &
                       +         ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
 
-                diss_n(k) = - ABS( v_comp(k) - gv ) * (                       &
+                diss_n = - ABS( v_comp - gv ) * (                             &
                       10.0_wp * ( v(k,j+1,i) - v(k,j,i)   )                      &
                     -  5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) )                      &
                     +           ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
+
+#ifdef _OPENACC
+!
+!--             Recompute the south fluxes.
+                v_comp_s = v(k,j,i) + v(k,j-1,i) - gv
+                flux_s   = v_comp_s * (                                       &
+                           37.0_wp * ( v(k,j,i) + v(k,j-1,i)   )                 &
+                         -  8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) )                 &
+                         +           ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
+                diss_s   = - ABS( v_comp_s ) * (                              &
+                           10.0_wp * ( v(k,j,i) - v(k,j-1,i)   )                 &
+                         -  5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) )                 &
+                         +           ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
+#else
+                flux_s = swap_flux_y_local_v(k)
+                diss_s = swap_diss_y_local_v(k)
+#endif
+
 !
 !--             k index has to be modified near bottom and top, else array
@@ -4745,6 +5371,6 @@
                 k_mm  = k - 2 * ibit26
 
-                w_comp    = w(k,j-1,i) + w(k,j,i)
-                flux_t(k) = w_comp * rho_air_zw(k) * (                       &
+                w_comp = w(k,j-1,i) + w(k,j,i)
+                flux_t = w_comp * rho_air_zw(k) * (                          &
                           ( 37.0_wp * ibit26 * adv_mom_5                        &
                        +     7.0_wp * ibit25 * adv_mom_3                        &
@@ -4761,5 +5387,5 @@
                                       )
 
-                diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * (              &
+                diss_t = - ABS( w_comp ) * rho_air_zw(k) * (                 &
                           ( 10.0_wp * ibit26 * adv_mom_5                        &
                        +     3.0_wp * ibit25 * adv_mom_3                        &
@@ -4780,5 +5406,5 @@
 !--             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(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
+                     +  ( v_comp      - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
                      +  (   w_comp                      * rho_air_zw(k) -     &
                           ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)     &
@@ -4787,40 +5413,41 @@
  
                 tend(k,j,i) = tend(k,j,i) - (                                 &
-                       ( flux_r(k) + diss_r(k)                                &
-                     -   swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j)  &
-                       ) * ddx                                                &
-                     + ( flux_n(k) + diss_n(k)                                &
-                     -   swap_flux_y_local_v(k) - swap_diss_y_local_v(k)      &
-                       ) * ddy                                                &
-                     + ( ( flux_t(k) + diss_t(k) )                            &
-                     -   ( flux_d    + diss_d    )                            &
-                       ) * drho_air(k) * ddzw(k)                              &
+                       ( ( flux_r + diss_r )                                  &
+                     -   ( flux_l + diss_l ) ) * ddx                          &
+                     + ( ( flux_n + diss_n )                                  &
+                     -   ( flux_s + diss_s ) ) * ddy                          &
+                     + ( ( flux_t + diss_t )                                  &
+                     -   ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k)        &
                                             )  + v(k,j,i) * div
 
-                swap_flux_x_local_v(k,j) = flux_r(k)
-                swap_diss_x_local_v(k,j) = diss_r(k)
-                swap_flux_y_local_v(k)   = flux_n(k)
-                swap_diss_y_local_v(k)   = diss_n(k)
-                flux_d                   = flux_t(k)
-                diss_d                   = diss_t(k)
+#ifndef _OPENACC
+                swap_flux_x_local_v(k,j) = flux_r
+                swap_diss_x_local_v(k,j) = diss_r
+                swap_flux_y_local_v(k)   = flux_n
+                swap_diss_y_local_v(k)   = diss_n
+#endif
+                flux_d                   = flux_t
+                diss_d                   = diss_t
 
 !
 !--             Statistical Evaluation of v'v'. The factor has to be applied
 !--             for right evaluation when gallilei_trans = .T. .
+                !$ACC ATOMIC
                 sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                      &
-                + ( flux_n(k)                                                  &
-                    * ( v_comp(k) - 2.0_wp * hom(k,1,2,0)                   )  &
-                    / ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) )  &
-               +   diss_n(k)                                                   &
-                    *   ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0)              )  &
-                    / ( ABS( v_comp(k) - gv ) + 1.0E-20_wp                  )  &
+                + ( flux_n                                                     &
+                    * ( v_comp - 2.0_wp * hom(k,1,2,0)                   )     &
+                    / ( v_comp - gv + SIGN( 1.0E-20_wp, v_comp - gv )    )     &
+               +   diss_n                                                      &
+                    *   ABS( v_comp - 2.0_wp * hom(k,1,2,0)              )     &
+                    / ( ABS( v_comp - gv ) + 1.0E-20_wp                  )     &
                   ) *   weight_substep(intermediate_timestep_count)
 !
 !--             Statistical Evaluation of w'u'.
+                !$ACC ATOMIC
                 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                    &
-                + ( flux_t(k)                                                  &
+                + ( flux_t                                                     &
                     * ( w_comp - 2.0_wp * hom(k,1,3,0)                   )     &
                     / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              )     &
-               +   diss_t(k)                                                   &
+               +   diss_t                                                      &
                     *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              )     &
                     / ( ABS( w_comp ) + 1.0E-20_wp                       )     &
@@ -4830,5 +5457,7 @@
           ENDDO
        ENDDO
+!$ACC UPDATE HOST(sums_vs2_ws_l(nzb+1,tn))
        sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
+!$ACC UPDATE DEVICE(sums_vs2_ws_l(nzb,tn))
 
 
@@ -4874,7 +5503,17 @@
        REAL(wp)    ::  ibit28 !< flag indicating 3rd-order scheme along x-direction
        REAL(wp)    ::  ibit29 !< flag indicating 5th-order scheme along x-direction
+#ifdef _OPENACC
+       REAL(wp)    ::  ibit27_l !< flag indicating 1st-order scheme along x-direction
+       REAL(wp)    ::  ibit28_l !< flag indicating 3rd-order scheme along x-direction
+       REAL(wp)    ::  ibit29_l !< flag indicating 5th-order scheme along x-direction
+#endif
        REAL(wp)    ::  ibit30 !< flag indicating 1st-order scheme along y-direction
        REAL(wp)    ::  ibit31 !< flag indicating 3rd-order scheme along y-direction
        REAL(wp)    ::  ibit32 !< flag indicating 5th-order scheme along y-direction
+#ifdef _OPENACC
+       REAL(wp)    ::  ibit30_s !< flag indicating 1st-order scheme along y-direction
+       REAL(wp)    ::  ibit31_s !< flag indicating 3rd-order scheme along y-direction
+       REAL(wp)    ::  ibit32_s !< flag indicating 5th-order scheme along y-direction
+#endif
        REAL(wp)    ::  ibit33 !< flag indicating 1st-order scheme along z-direction
        REAL(wp)    ::  ibit34 !< flag indicating 3rd-order scheme along z-direction
@@ -4886,22 +5525,39 @@
        REAL(wp)    ::  gv     !< Galilei-transformation velocity along y
        REAL(wp)    ::  u_comp !< advection velocity along x
+#ifdef _OPENACC
+       REAL(wp)    ::  u_comp_l !< advection velocity along x
+#endif
        REAL(wp)    ::  v_comp !< advection velocity along y
+#ifdef _OPENACC
+       REAL(wp)    ::  v_comp_s !< advection velocity along y
+#endif
        REAL(wp)    ::  w_comp !< advection velocity along z
        
-       REAL(wp), DIMENSION(nzb:nzt)    ::  diss_t !< discretized artificial dissipation at top of the grid box
-       REAL(wp), DIMENSION(nzb:nzt)    ::  flux_t !< discretized 6th-order flux at top of the grid box
+       REAL(wp)    ::  diss_t !< discretized artificial dissipation at top of the grid box
+       REAL(wp)    ::  flux_t !< discretized 6th-order flux at top of the grid box
        
-       REAL(wp), DIMENSION(nzb+1:nzt)  ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
-       REAL(wp), DIMENSION(nzb+1:nzt)  ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
-       REAL(wp), DIMENSION(nzb+1:nzt)  ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
-       REAL(wp), DIMENSION(nzb+1:nzt)  ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
+       REAL(wp)    ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
+       REAL(wp)    ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
+       REAL(wp)    ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
+       REAL(wp)    ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
+
+       REAL(wp)    ::  diss_s !< discretized artificial dissipation at southward-side of the grid box
+       REAL(wp)    ::  flux_s !< discretized 6th-order flux at southward-side of the grid box
+#ifndef _OPENACC
        REAL(wp), DIMENSION(nzb+1:nzt)  ::  swap_diss_y_local_w !< discretized artificial dissipation at southward-side of the grid box
        REAL(wp), DIMENSION(nzb+1:nzt)  ::  swap_flux_y_local_w !< discretized 6th-order flux at southward-side of the grid box
+#endif
        
+       REAL(wp)    ::  diss_l !< discretized artificial dissipation at leftward-side of the grid box
+       REAL(wp)    ::  flux_l !< discretized 6th-order flux at leftward-side of the grid box
+#ifndef _OPENACC
        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_diss_x_local_w !< discretized artificial dissipation at leftward-side of the grid box
        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_flux_x_local_w !< discretized 6th-order flux at leftward-side of the grid box
+#endif
  
        gu = 2.0_wp * u_gtrans
        gv = 2.0_wp * v_gtrans
+
+#ifndef _OPENACC
 !
 !--   compute the whole left boundary of the processor domain
@@ -4962,7 +5618,26 @@
 
        ENDDO
-
+#endif
+
+       !$ACC PARALLEL LOOP COLLAPSE(2) FIRSTPRIVATE(tn, gu, gv) &
+       !$ACC PRIVATE(i, j, k, k_mm, k_pp, k_ppp) &
+       !$ACC PRIVATE(ibit27, ibit28, ibit29, ibit30, ibit31, ibit32) &
+       !$ACC PRIVATE(ibit27_l, ibit28_l, ibit29_l) &
+       !$ACC PRIVATE(ibit30_s, ibit31_s, ibit32_s) &
+       !$ACC PRIVATE(ibit33, ibit34, ibit35) &
+       !$ACC PRIVATE(flux_r, diss_r, flux_l, diss_l) &
+       !$ACC PRIVATE(flux_n, diss_n, flux_s, diss_s) &
+       !$ACC PRIVATE(flux_t, diss_t, flux_d, diss_d) &
+       !$ACC PRIVATE(div, u_comp, u_comp_l, v_comp, v_comp_s, w_comp) &
+       !$ACC PRESENT(advc_flags_1, advc_flags_2) &
+       !$ACC PRESENT(u, v, w) &
+       !$ACC PRESENT(rho_air, drho_air_zw, ddzu) &
+       !$ACC PRESENT(tend) &
+       !$ACC PRESENT(hom(nzb+1:nzb_max,1,3,0)) &
+       !$ACC PRESENT(weight_substep(intermediate_timestep_count)) &
+       !$ACC PRESENT(sums_ws2_ws_l(nzb+1:nzb_max,0))
        DO i = nxl, nxr
 
+#ifndef _OPENACC
           j = nys
           DO  k = nzb+1, nzb_max
@@ -5018,4 +5693,5 @@
 
           ENDDO
+#endif
 
           DO  j = nys, nyn
@@ -5025,10 +5701,8 @@
 !--          at the first w-level. For topography wall this is done implicitely
 !--          by advc_flags_1.
-             k         = nzb + 1
-             w_comp    = w(k,j,i) + w(k-1,j,i)
-             flux_t(0) = w_comp       * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
-             diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
-             flux_d    = flux_t(0)
-             diss_d    = diss_t(0)
+             k      = nzb + 1
+             w_comp = w(k,j,i) + w(k-1,j,i)
+             flux_d = w_comp       * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
+             diss_d = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
 
              DO  k = nzb+1, nzb_max
@@ -5038,6 +5712,6 @@
                 ibit27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp )
 
-                u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
-                flux_r(k) = u_comp * (                                       &
+                u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
+                flux_r = u_comp * (                                          &
                           ( 37.0_wp * ibit29 * adv_mom_5                        &
                        +     7.0_wp * ibit28 * adv_mom_3                        &
@@ -5054,5 +5728,5 @@
                                      )
 
-                diss_r(k) = - ABS( u_comp ) * (                              &
+                diss_r = - ABS( u_comp ) * (                                 &
                           ( 10.0_wp * ibit29 * adv_mom_5                        &
                        +     3.0_wp * ibit28 * adv_mom_3                        &
@@ -5069,10 +5743,53 @@
                                               )
 
+#ifdef _OPENACC
+!
+!--             Recompute the left fluxes.
+                ibit29_l = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )
+                ibit28_l = REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp )
+                ibit27_l = REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp )
+
+                u_comp_l = u(k+1,j,i) + u(k,j,i) - gu
+                flux_l   = u_comp_l * (                                      &
+                                      ( 37.0_wp * ibit29_l * adv_mom_5          &
+                                   +     7.0_wp * ibit28_l * adv_mom_3          &
+                                   +              ibit27_l * adv_mom_1          &
+                                      ) *                                    &
+                                     ( w(k,j,i)   + w(k,j,i-1) )             &
+                               -      (  8.0_wp * ibit29_l * adv_mom_5          &
+                                   +              ibit28_l * adv_mom_3          &
+                                      ) *                                    &
+                                     ( w(k,j,i+1) + w(k,j,i-2) )             &
+                               +      (           ibit29_l * adv_mom_5          &
+                                      ) *                                    &
+                                     ( w(k,j,i+2) + w(k,j,i-3) )             &
+                                                 )
+
+                diss_l   = - ABS( u_comp_l ) * (                             &
+                                        ( 10.0_wp * ibit29_l * adv_mom_5        &
+                                     +     3.0_wp * ibit28_l * adv_mom_3        &
+                                     +              ibit27_l * adv_mom_1        &
+                                        ) *                                  &
+                                     ( w(k,j,i)   - w(k,j,i-1) )             &
+                                 -      (  5.0_wp * ibit29_l * adv_mom_5        &
+                                     +              ibit28_l * adv_mom_3        &
+                                        ) *                                  &
+                                     ( w(k,j,i+1) - w(k,j,i-2) )             &
+                                 +      (           ibit29_l * adv_mom_5        &
+                                        ) *                                  &
+                                     ( w(k,j,i+2) - w(k,j,i-3) )             &
+                                                            )
+#else
+                flux_l = swap_flux_x_local_w(k,j)
+                diss_l = swap_diss_x_local_w(k,j)
+#endif
+
+
                 ibit32 = REAL( IBITS(advc_flags_2(k,j,i),0,1),  KIND = wp )
                 ibit31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp )
                 ibit30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp )
 
-                v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
-                flux_n(k) = v_comp * (                                       &
+                v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
+                flux_n = v_comp * (                                          &
                           ( 37.0_wp * ibit32 * adv_mom_5                        &
                        +     7.0_wp * ibit31 * adv_mom_3                        &
@@ -5089,5 +5806,5 @@
                                      )
 
-                diss_n(k) = - ABS( v_comp ) * (                              &
+                diss_n = - ABS( v_comp ) * (                                 &
                           ( 10.0_wp * ibit32 * adv_mom_5                        &
                        +     3.0_wp * ibit31 * adv_mom_3                        &
@@ -5103,4 +5820,47 @@
                                  ( w(k,j+3,i) - w(k,j-2,i) )                 &
                                               )
+
+#ifdef _OPENACC
+!
+!--             Recompute the south fluxes.
+                ibit32_s = REAL( IBITS(advc_flags_2(k,j-1,i),0,1),  KIND = wp )
+                ibit31_s = REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )
+                ibit30_s = REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp )
+
+                v_comp_s = v(k+1,j,i) + v(k,j,i) - gv
+                flux_s   = v_comp_s * (                                      &
+                                    ( 37.0_wp * ibit32_s * adv_mom_5            &
+                                 +     7.0_wp * ibit31_s * adv_mom_3            &
+                                 +              ibit30_s * adv_mom_1            &
+                                    ) *                                      &
+                                     ( w(k,j,i)   + w(k,j-1,i) )             &
+                             -      (  8.0_wp * ibit32_s * adv_mom_5            &
+                                 +              ibit31_s * adv_mom_3            &
+                                    ) *                                      &
+                                     ( w(k,j+1,i) + w(k,j-2,i) )             &
+                             +      (           ibit32_s * adv_mom_5            &
+                                    ) *                                      &
+                                     ( w(k,j+2,i) + w(k,j-3,i) )             &
+                                               )
+
+                diss_s   = - ABS( v_comp_s ) * (                             &
+                                    ( 10.0_wp * ibit32_s * adv_mom_5            &
+                                 +     3.0_wp * ibit31_s * adv_mom_3            &
+                                 +              ibit30_s * adv_mom_1            &
+                                    ) *                                      &
+                                     ( w(k,j,i)   - w(k,j-1,i) )             &
+                             -      (  5.0_wp * ibit32_s * adv_mom_5            &
+                                 +              ibit31_s * adv_mom_3            &
+                                    ) *                                      &
+                                     ( w(k,j+1,i) - w(k,j-2,i) )             &
+                             +      (           ibit32_s * adv_mom_5            &
+                                    ) *                                      &
+                                     ( w(k,j+2,i) - w(k,j-3,i) )             &
+                                                        )
+#else
+                flux_s = swap_flux_y_local_w(k)
+                diss_s = swap_diss_y_local_w(k)
+#endif
+
 !
 !--             k index has to be modified near bottom and top, else array
@@ -5114,6 +5874,6 @@
                 k_mm  = k - 2 * ibit35
 
-                w_comp    = w(k+1,j,i) + w(k,j,i)
-                flux_t(k) = w_comp * rho_air(k+1) * (                        &
+                w_comp = w(k+1,j,i) + w(k,j,i)
+                flux_t = w_comp * rho_air(k+1) * (                           &
                           ( 37.0_wp * ibit35 * adv_mom_5                        &
                        +     7.0_wp * ibit34 * adv_mom_3                        &
@@ -5130,5 +5890,5 @@
                                        )
 
-                diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (               &
+                diss_t = - ABS( w_comp ) * rho_air(k+1) * (                  &
                           ( 10.0_wp * ibit35 * adv_mom_5                        &
                        +     3.0_wp * ibit34 * adv_mom_3                        &
@@ -5177,27 +5937,27 @@
 
                 tend(k,j,i) = tend(k,j,i) - (                                 &
-                      ( flux_r(k) + diss_r(k)                                 &
-                    -   swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j)   &
-                      ) * ddx                                                 &
-                    + ( flux_n(k) + diss_n(k)                                 &
-                    -   swap_flux_y_local_w(k)   - swap_diss_y_local_w(k)     &
-                      ) * ddy                                                 &
-                    + ( ( flux_t(k) + diss_t(k) )                             &
-                    -   ( flux_d    + diss_d    )                             &
-                      ) * drho_air_zw(k) * ddzu(k+1)                          &
+                      ( ( flux_r + diss_r )                                   &
+                    -   ( flux_l + diss_l ) ) * ddx                           &
+                    + ( ( flux_n + diss_n )                                   &
+                    -   ( flux_s + diss_s ) ) * ddy                           &
+                    + ( ( flux_t + diss_t )                                   &
+                    -   ( flux_d + diss_d ) ) * drho_air_zw(k) * ddzu(k+1)    &
                                             )  + div * w(k,j,i)
 
-                swap_flux_x_local_w(k,j) = flux_r(k)
-                swap_diss_x_local_w(k,j) = diss_r(k)
-                swap_flux_y_local_w(k)   = flux_n(k)
-                swap_diss_y_local_w(k)   = diss_n(k)
-                flux_d                   = flux_t(k)
-                diss_d                   = diss_t(k)
-
+#ifndef _OPENACC
+                swap_flux_x_local_w(k,j) = flux_r
+                swap_diss_x_local_w(k,j) = diss_r
+                swap_flux_y_local_w(k)   = flux_n
+                swap_diss_y_local_w(k)   = diss_n
+#endif
+                flux_d                   = flux_t
+                diss_d                   = diss_t
+
+                !$ACC ATOMIC
                 sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                    &
-                      + ( flux_t(k)                                           &
+                      + ( flux_t                                              &
                        * ( w_comp - 2.0_wp * hom(k,1,3,0)                   ) &
                        / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              ) &
-                        + diss_t(k)                                           &
+                        + diss_t                                              &
                        *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              ) &
                        / ( ABS( w_comp ) + 1.0E-20_wp                       ) &
@@ -5208,25 +5968,60 @@
              DO  k = nzb_max+1, nzt
 
-                u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
-                flux_r(k) = u_comp * (                                      &
+                u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
+                flux_r = u_comp * (                                         &
                       37.0_wp * ( w(k,j,i+1) + w(k,j,i)   )                    &
                     -  8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) )                    &
                     +           ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
 
-                diss_r(k) = - ABS( u_comp ) * (                             &
+                diss_r = - ABS( u_comp ) * (                                &
                       10.0_wp * ( w(k,j,i+1) - w(k,j,i)   )                    &
                     -  5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) )                    &
                     +           ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
 
-                v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
-                flux_n(k) = v_comp * (                                      &
+#ifdef _OPENACC
+!
+!--             Recompute the left fluxes.
+                u_comp_l = u(k+1,j,i) + u(k,j,i) - gu
+                flux_l   = u_comp_l * (                                     &
+                            37.0_wp * ( w(k,j,i) + w(k,j,i-1)   )              &
+                          -  8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) )              &
+                          +           ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
+                diss_l   = - ABS( u_comp_l ) * (                            &
+                            10.0_wp * ( w(k,j,i) - w(k,j,i-1)   )              &
+                          -  5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) )              &
+                          +           ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
+#else
+                flux_l = swap_flux_x_local_w(k,j)
+                diss_l = swap_diss_x_local_w(k,j)
+#endif
+
+                v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
+                flux_n = v_comp * (                                         &
                       37.0_wp * ( w(k,j+1,i) + w(k,j,i)   )                    &
                     -  8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) )                    &
                     +           ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
 
-                diss_n(k) = - ABS( v_comp ) * (                             &
+                diss_n = - ABS( v_comp ) * (                                &
                       10.0_wp * ( w(k,j+1,i) - w(k,j,i)   )                    &
                     -  5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) )                    &
                     +           ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
+
+#ifdef _OPENACC
+!
+!--             Recompute the south fluxes.
+                v_comp_s = v(k+1,j,i) + v(k,j,i) - gv
+                flux_s   = v_comp_s * (                                     &
+                           37.0_wp * ( w(k,j,i) + w(k,j-1,i)   )               &
+                         -  8.0_wp * ( w(k,j+1,i) +w(k,j-2,i)  )               &
+                         +           ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
+                diss_s   = - ABS( v_comp_s ) * (                            &
+                           10.0_wp * ( w(k,j,i) - w(k,j-1,i)   )               &
+                         -  5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) )               &
+                         +           ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
+#else
+                flux_s = swap_flux_y_local_w(k)
+                diss_s = swap_diss_y_local_w(k)
+#endif
+
 !
 !--             k index has to be modified near bottom and top, else array
@@ -5240,6 +6035,6 @@
                 k_mm  = k - 2 * ibit35
 
-                w_comp    = w(k+1,j,i) + w(k,j,i)
-                flux_t(k) = w_comp * rho_air(k+1) * (                        &
+                w_comp = w(k+1,j,i) + w(k,j,i)
+                flux_t = w_comp * rho_air(k+1) * (                           &
                           ( 37.0_wp * ibit35 * adv_mom_5                        &
                        +     7.0_wp * ibit34 * adv_mom_3                        &
@@ -5256,5 +6051,5 @@
                                        )
 
-                diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * (               &
+                diss_t = - ABS( w_comp ) * rho_air(k+1) * (                  &
                           ( 10.0_wp * ibit35 * adv_mom_5                        &
                        +     3.0_wp * ibit34 * adv_mom_3                        &
@@ -5282,27 +6077,27 @@
 
                 tend(k,j,i) = tend(k,j,i) - (                                 &
-                      ( flux_r(k) + diss_r(k)                                 &
-                    -   swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j)   &
-                      ) * ddx                                                 &
-                    + ( flux_n(k) + diss_n(k)                                 &
-                    -   swap_flux_y_local_w(k)   - swap_diss_y_local_w(k)     &
-                      ) * ddy                                                 &
-                    + ( ( flux_t(k) + diss_t(k) )                             &
-                    -   ( flux_d    + diss_d    )                             &
-                      ) * drho_air_zw(k) * ddzu(k+1)                          &
+                      ( ( flux_r + diss_r )                                   &
+                    -   ( flux_l + diss_l ) ) * ddx                           &
+                    + ( ( flux_n + diss_n )                                   &
+                    -   ( flux_s + diss_s ) ) * ddy                           &
+                    + ( ( flux_t + diss_t )                                   &
+                    -   ( flux_d + diss_d ) ) * drho_air_zw(k) * ddzu(k+1)    &
                                             )  + div * w(k,j,i)
 
-                swap_flux_x_local_w(k,j) = flux_r(k)
-                swap_diss_x_local_w(k,j) = diss_r(k)
-                swap_flux_y_local_w(k)   = flux_n(k)
-                swap_diss_y_local_w(k)   = diss_n(k)
-                flux_d                   = flux_t(k)
-                diss_d                   = diss_t(k)
-
+#ifndef _OPENACC
+                swap_flux_x_local_w(k,j) = flux_r
+                swap_diss_x_local_w(k,j) = diss_r
+                swap_flux_y_local_w(k)   = flux_n
+                swap_diss_y_local_w(k)   = diss_n
+#endif
+                flux_d                   = flux_t
+                diss_d                   = diss_t
+
+                !$ACC ATOMIC
                 sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                    &
-                      + ( flux_t(k)                                           &
+                      + ( flux_t                                              &
                        * ( w_comp - 2.0_wp * hom(k,1,3,0)                   ) &
                        / ( w_comp + SIGN( 1.0E-20_wp, w_comp )              ) &
-                        + diss_t(k)                                           &
+                        + diss_t                                              &
                        *   ABS( w_comp - 2.0_wp * hom(k,1,3,0)              ) &
                        / ( ABS( w_comp ) + 1.0E-20_wp                       ) &
Index: palm/trunk/SOURCE/basic_constants_and_equations_mod.f90
===================================================================
--- palm/trunk/SOURCE/basic_constants_and_equations_mod.f90	(revision 3633)
+++ palm/trunk/SOURCE/basic_constants_and_equations_mod.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3449 2018-10-29 19:36:56Z suehring
 ! +degc_to_k 
 ! 
@@ -58,4 +61,5 @@
     REAL(wp), PARAMETER ::  molecular_weight_of_water = 0.01801528_wp !< mol. m. H2O (kg mol-1)
     REAL(wp), PARAMETER ::  pi = 3.141592654_wp                       !< PI
+    !$ACC DECLARE COPYIN(pi)
     REAL(wp), PARAMETER ::  rho_l = 1.0E3_wp                          !< density of water (kg m-3)
     REAL(wp), PARAMETER ::  rho_nacl = 2165.0_wp                      !< density of NaCl (kg m-3)
Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 3633)
+++ palm/trunk/SOURCE/boundary_conds.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3589 2018-11-30 15:09:51Z suehring
 ! Move the control parameter "salsa" from salsa_mod to control_parameters 
 ! (M. Kurppa)
@@ -289,4 +292,6 @@
        kb = MERGE( -1, 1, l == 0 )
        !$OMP PARALLEL DO PRIVATE( i, j, k )
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+       !$ACC PRESENT(bc_h, w_p)
        DO  m = 1, bc_h(l)%ns
           i = bc_h(l)%i(m)            
@@ -300,6 +305,8 @@
 !-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
     IF ( ibc_uv_t == 0 )  THEN
+        !$ACC KERNELS PRESENT(u_p, v_p, u_init, v_init)
         u_p(nzt+1,:,:) = u_init(nzt+1)
         v_p(nzt+1,:,:) = v_init(nzt+1)
+        !$ACC END KERNELS
     ELSEIF ( ibc_uv_t == 1 )  THEN
         u_p(nzt+1,:,:) = u_p(nzt,:,:)
@@ -311,5 +318,7 @@
     IF (  .NOT.  child_domain  .AND.  .NOT.  nesting_offline  .AND.            &
                  TRIM(coupling_mode) /= 'vnested_fine' )  THEN
+       !$ACC KERNELS PRESENT(w_p)
        w_p(nzt:nzt+1,:,:) = 0.0_wp  !< nzt is not a prognostic level (but cf. pres)
+       !$ACC END KERNELS
     ENDIF
 
@@ -342,4 +351,6 @@
           kb = MERGE( -1, 1, l == 0 )
           !$OMP PARALLEL DO PRIVATE( i, j, k )
+          !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+          !$ACC PRESENT(bc_h, pt_p)
           DO  m = 1, bc_h(l)%ns
              i = bc_h(l)%i(m)            
@@ -364,5 +375,7 @@
         pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
     ELSEIF ( ibc_pt_t == 2 )  THEN
+        !$ACC KERNELS PRESENT(pt_p, dzu)
         pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
+        !$ACC END KERNELS
     ENDIF
 
@@ -379,4 +392,6 @@
              kb = MERGE( -1, 1, l == 0 )
              !$OMP PARALLEL DO PRIVATE( i, j, k )
+             !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+             !$ACC PRESENT(bc_h, e_p)
              DO  m = 1, bc_h(l)%ns
                 i = bc_h(l)%i(m)            
@@ -444,5 +459,7 @@
 
        IF ( .NOT. child_domain )  THEN
+          !$ACC KERNELS PRESENT(e_p)
           e_p(nzt+1,:,:) = e_p(nzt,:,:)
+          !$ACC END KERNELS
 !
 !--    Nesting case: if parent operates in RANS mode and child in LES mode,
Index: palm/trunk/SOURCE/buoyancy.f90
===================================================================
--- palm/trunk/SOURCE/buoyancy.f90	(revision 3633)
+++ palm/trunk/SOURCE/buoyancy.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3538 2018-11-20 10:55:41Z suehring
 ! Remove unnecessary double-masking of topography
 ! 
@@ -175,4 +178,8 @@
 !
 !--       Normal case: horizontal surface
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+          !$ACC PRESENT(var) &
+          !$ACC PRESENT(ref_state) &
+          !$ACC PRESENT(tend)
           DO  i = nxl, nxr
              DO  j = nys, nyn
Index: palm/trunk/SOURCE/coriolis.f90
===================================================================
--- palm/trunk/SOURCE/coriolis.f90	(revision 3633)
+++ palm/trunk/SOURCE/coriolis.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3538 2018-11-20 10:55:41Z suehring
 ! Note concerning topography masking added
 ! 
@@ -145,4 +148,8 @@
 !--       u-component
           CASE ( 1 )
+             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag) &
+             !$ACC PRESENT(wall_flags_0) &
+             !$ACC PRESENT(v, w, vg) &
+             !$ACC PRESENT(tend)
              DO  i = nxlu, nxr
                 DO  j = nys, nyn
@@ -167,4 +174,8 @@
 !--       v-component
           CASE ( 2 )
+             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag) &
+             !$ACC PRESENT(wall_flags_0) &
+             !$ACC PRESENT(u, ug) &
+             !$ACC PRESENT(tend)
              DO  i = nxl, nxr
                 DO  j = nysv, nyn
@@ -185,4 +196,8 @@
 !--       w-component
           CASE ( 3 )
+             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag) &
+             !$ACC PRESENT(wall_flags_0) &
+             !$ACC PRESENT(u) &
+             !$ACC PRESENT(tend)
              DO  i = nxl, nxr
                 DO  j = nys, nyn
Index: palm/trunk/SOURCE/cuda_fft_interfaces.f90
===================================================================
--- palm/trunk/SOURCE/cuda_fft_interfaces.f90	(revision 3634)
+++ palm/trunk/SOURCE/cuda_fft_interfaces.f90	(revision 3634)
@@ -0,0 +1,160 @@
+!> @file cuda_fft_interfaces.f90
+!--------------------------------------------------------------------------------!
+! This file is part of PALM.
+!
+! PALM is free software: you can redistribute it and/or modify it under the terms
+! of the GNU General Public License as published by the Free Software Foundation,
+! either version 3 of the License, or (at your option) any later version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2018 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------!
+!
+! Current revisions:
+! -----------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
+! Initial revision ( OpenACC port for SPEC )
+! 
+! 1683 2015-10-07 23:57:51Z knoop
+!
+!
+! Description:
+! ------------
+!> FORTRAN interfaces for the CUDA fft
+!> Routines for the fft along x and y (forward/backward) using the CUDA fft
+!--------------------------------------------------------------------------------!
+ MODULE cuda_fft_interfaces
+ 
+
+#if defined ( __cuda_fft )
+
+    USE kinds
+
+    INTEGER(iwp) ::  CUFFT_FORWARD = -1   !<
+    INTEGER(iwp) ::  CUFFT_INVERSE =  1   !<
+    INTEGER(iwp) ::  CUFFT_R2C = Z'2a'    !< Real to Complex (interleaved)
+    INTEGER(iwp) ::  CUFFT_C2R = Z'2c'    !< Complex (interleaved) to Real
+    INTEGER(iwp) ::  CUFFT_C2C = Z'29'    !< Complex to Complex, interleaved
+    INTEGER(iwp) ::  CUFFT_D2Z = Z'6a'    !< Double to Double-Complex
+    INTEGER(iwp) ::  CUFFT_Z2D = Z'6c'    !< Double-Complex to Double
+    INTEGER(iwp) ::  CUFFT_Z2Z = Z'69'    !< Double-Complex to Double-Complex
+
+    PUBLIC
+
+
+!
+!-- cufftPlan1d( cufftHandle *plan, int nx, cufftType type, int batch )
+    INTERFACE CUFFTPLAN1D
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> @todo Missing subroutine description.
+!------------------------------------------------------------------------------!
+       SUBROUTINE CUFFTPLAN1D( plan, nx, type, batch ) bind( C, name='cufftPlan1d' )
+
+          USE ISO_C_BINDING
+
+          INTEGER(C_INT)        ::  plan   !< 
+          INTEGER(C_INT), value ::  batch  !< 
+          INTEGER(C_INT), value ::  nx     !< 
+          INTEGER(C_INT), value ::  type   !< 
+       END SUBROUTINE CUFFTPLAN1D
+
+    END INTERFACE CUFFTPLAN1D
+
+!
+!-- cufftDestroy( cufftHandle plan )  !!! remove later if not really needed !!!
+    INTERFACE CUFFTDESTROY
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> @todo Missing subroutine description.
+!------------------------------------------------------------------------------!
+       SUBROUTINE CUFFTDESTROY( plan ) bind( C, name='cufftDestroy' )
+
+          USE ISO_C_BINDING
+
+          INTEGER(C_INT), VALUE ::  plan
+
+       END SUBROUTINE CUFFTDESTROY
+
+    END INTERFACE CUFFTDESTROY
+
+
+    INTERFACE CUFFTEXECZ2D
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> @todo Missing subroutine description.
+!------------------------------------------------------------------------------!
+       SUBROUTINE CUFFTEXECZ2D( plan, idata, odata ) bind( C, name='cufftExecZ2D' )
+
+          USE ISO_C_BINDING
+          USE kinds
+
+          INTEGER(C_INT), VALUE ::  plan          !<
+          COMPLEX(dp), DEVICE   ::  idata(:,:,:)  !<
+          REAL(dp), DEVICE      ::  odata(:,:,:)  !<
+
+       END SUBROUTINE CUFFTEXECZ2D
+
+    END INTERFACE CUFFTEXECZ2D
+
+
+    INTERFACE CUFFTEXECD2Z
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> @todo Missing subroutine description.
+!------------------------------------------------------------------------------!
+       SUBROUTINE CUFFTEXECD2Z( plan, idata, odata ) bind( C, name='cufftExecD2Z' )
+
+          USE ISO_C_BINDING
+          
+          USE kinds
+
+          INTEGER(C_INT), VALUE ::  plan          !<
+          REAL(dp), DEVICE      ::  idata(:,:,:)  !<
+          COMPLEX(dp), DEVICE   ::  odata(:,:,:)  !<
+
+       END SUBROUTINE CUFFTEXECD2Z
+
+    END INTERFACE CUFFTEXECD2Z
+
+#else
+
+    INTERFACE CUFFTdummy
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Dummy interface to avoid compiler warnings in case of no bublic objects
+!> declared.
+!------------------------------------------------------------------------------!
+       SUBROUTINE CUFFTdummy( dummy )
+       
+          USE kinds
+
+          REAL(wp) ::  dummy  !<
+
+       END SUBROUTINE CUFFTdummy
+
+    END INTERFACE CUFFTdummy
+
+#endif
+
+ END MODULE cuda_fft_interfaces
Index: palm/trunk/SOURCE/diffusion_s.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_s.f90	(revision 3633)
+++ palm/trunk/SOURCE/diffusion_s.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3547 2018-11-21 13:21:24Z suehring
 ! variables documented
 ! 
@@ -194,4 +197,24 @@
 #endif
 
+       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, m) &
+       !$ACC PRIVATE(surf_e, surf_s, flag, mask_top, mask_bottom) &
+       !$ACC PRIVATE(mask_north, mask_south, mask_west, mask_east) &
+       !$ACC PRESENT(wall_flags_0, kh) &
+       !$ACC PRESENT(s) &
+       !$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
+       !$ACC PRESENT(surf_def_h(0:2), surf_def_v(0:3)) &
+       !$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) &
+       !$ACC PRESENT(surf_usm_h, surf_usm_v(0:3)) &
+       !$ACC PRESENT(s_flux_def_h_up, s_flux_def_h_down) &
+       !$ACC PRESENT(s_flux_t) &
+       !$ACC PRESENT(s_flux_def_v_north, s_flux_def_v_south) &
+       !$ACC PRESENT(s_flux_def_v_east, s_flux_def_v_west) &
+       !$ACC PRESENT(s_flux_lsm_h_up) &
+       !$ACC PRESENT(s_flux_lsm_v_north, s_flux_lsm_v_south) &
+       !$ACC PRESENT(s_flux_lsm_v_east, s_flux_lsm_v_west) &
+       !$ACC PRESENT(s_flux_usm_h_up) &
+       !$ACC PRESENT(s_flux_usm_v_north, s_flux_usm_v_south) &
+       !$ACC PRESENT(s_flux_usm_v_east, s_flux_usm_v_west) &
+       !$ACC PRESENT(tend)
        DO  i = nxl, nxr
           DO  j = nys,nyn
Index: palm/trunk/SOURCE/diffusion_u.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_u.f90	(revision 3633)
+++ palm/trunk/SOURCE/diffusion_u.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3547 2018-11-21 13:21:24Z suehring
 ! variables documented
 ! 
@@ -172,4 +175,14 @@
 
 
+       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
+       !$ACC PRIVATE(surf_e, surf_s, flag, kmym, kmyp, kmzm, kmzp) &
+       !$ACC PRIVATE(mask_bottom, mask_north, mask_south, mask_top) &
+       !$ACC PRESENT(wall_flags_0, km) &
+       !$ACC PRESENT(u, v, w) &
+       !$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
+       !$ACC PRESENT(surf_def_h(0:2), surf_def_v(0:1)) &
+       !$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:1)) &
+       !$ACC PRESENT(surf_usm_h, surf_usm_v(0:1)) &
+       !$ACC PRESENT(tend)
        DO  i = nxlu, nxr
           DO  j = nys, nyn
Index: palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_v.f90	(revision 3633)
+++ palm/trunk/SOURCE/diffusion_v.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3547 2018-11-21 13:21:24Z suehring
 ! variables documented
 ! 
@@ -165,4 +168,14 @@
        REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface      
 
+       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
+       !$ACC PRIVATE(surf_e, surf_s, flag, kmxm, kmxp, kmzm, kmzp) &
+       !$ACC PRIVATE(mask_bottom, mask_east, mask_west, mask_top) &
+       !$ACC PRESENT(wall_flags_0, km) &
+       !$ACC PRESENT(u, v, w) &
+       !$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
+       !$ACC PRESENT(surf_def_h(0:2), surf_def_v(2:3)) &
+       !$ACC PRESENT(surf_lsm_h, surf_lsm_v(2:3)) &
+       !$ACC PRESENT(surf_usm_h, surf_usm_v(0:3)) &
+       !$ACC PRESENT(tend)
        DO  i = nxl, nxr
           DO  j = nysv, nyn
Index: palm/trunk/SOURCE/diffusion_w.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_w.f90	(revision 3633)
+++ palm/trunk/SOURCE/diffusion_w.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3547 2018-11-21 13:21:24Z suehring
 ! variables documented
 ! 
@@ -166,4 +169,14 @@
 
 
+       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
+       !$ACC PRIVATE(surf_e, surf_s, flag, kmxm, kmxp, kmym, kmyp) &
+       !$ACC PRIVATE(mask_west, mask_east, mask_south, mask_north) &
+       !$ACC PRESENT(wall_flags_0, km) &
+       !$ACC PRESENT(u, v, w) &
+       !$ACC PRESENT(ddzu, ddzw, rho_air, drho_air_zw) &
+       !$ACC PRESENT(surf_def_v(0:3)) &
+       !$ACC PRESENT(surf_lsm_v(0:3)) &
+       !$ACC PRESENT(surf_usm_v(0:3)) &
+       !$ACC PRESENT(tend)
        DO  i = nxl, nxr
           DO  j = nys, nyn
Index: palm/trunk/SOURCE/exchange_horiz.f90
===================================================================
--- palm/trunk/SOURCE/exchange_horiz.f90	(revision 3633)
+++ palm/trunk/SOURCE/exchange_horiz.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3241 2018-09-12 15:02:00Z raasch
 ! unused variables removed
 ! 
@@ -117,4 +120,7 @@
 
 
+#ifdef _OPENACC
+    INTEGER(iwp) ::  i           !<
+#endif
     INTEGER(iwp) ::  nbgp_local  !<
     
@@ -124,4 +130,15 @@
 
     CALL cpu_log( log_point_s(2), 'exchange_horiz', 'start' )
+
+#ifdef _OPENACC
+    !$ACC UPDATE IF_PRESENT &
+    !$ACC HOST(ar(:,:,nxr-nbgp_local+1:nxr)) &
+    !$ACC HOST(ar(:,:,nxl:nxl+nbgp_local-1))
+    DO i = nxl-nbgp_local, nxr+nbgp_local
+       !$ACC UPDATE IF_PRESENT &
+       !$ACC HOST(ar(:,nyn-nbgp_local+1:nyn,i)) &
+       !$ACC HOST(ar(:,nys:nys+nbgp_local-1,i))
+    ENDDO
+#endif
 
 #if defined( __parallel )
@@ -263,4 +280,16 @@
 
 #endif
+
+#ifdef _OPENACC
+    !$ACC UPDATE IF_PRESENT &
+    !$ACC DEVICE(ar(:,:,nxl-nbgp_local:nxl-1)) &
+    !$ACC DEVICE(ar(:,:,nxr+1:nxr+nbgp_local))
+    DO i = nxl-nbgp_local, nxr+nbgp_local
+       !$ACC UPDATE IF_PRESENT &
+       !$ACC DEVICE(ar(:,nys-nbgp_local:nys-1,i)) &
+       !$ACC DEVICE(ar(:,nyn+1:nyn+nbgp_local,i))
+    ENDDO
+#endif
+
     CALL cpu_log( log_point_s(2), 'exchange_horiz', 'stop' )
 
Index: palm/trunk/SOURCE/fft_xy_mod.f90
===================================================================
--- palm/trunk/SOURCE/fft_xy_mod.f90	(revision 3633)
+++ palm/trunk/SOURCE/fft_xy_mod.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3241 2018-09-12 15:02:00Z raasch
 ! preprocessor switches for variables that are required on NEC only
 ! 
@@ -156,8 +159,12 @@
         ONLY:  fft_method, message_string
         
+    USE cuda_fft_interfaces
+        
     USE indices,                                                               &
         ONLY:  nx, ny, nz
         
-#if defined( __fftw )
+#if defined( __cuda_fft )
+    USE ISO_C_BINDING
+#elif defined( __fftw )
     USE, INTRINSIC ::  ISO_C_BINDING
 #endif
@@ -210,4 +217,10 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trig_yf  !<
     
+#elif defined( __cuda_fft )
+    INTEGER(C_INT), SAVE ::  plan_xf  !<
+    INTEGER(C_INT), SAVE ::  plan_xi  !<
+    INTEGER(C_INT), SAVE ::  plan_yf  !<
+    INTEGER(C_INT), SAVE ::  plan_yi  !<
+
 #endif
 
@@ -298,4 +311,8 @@
        ENDIF
 
+#if defined( _OPENACC ) && defined( __cuda_fft )
+       fft_method = 'system-specific'
+#endif
+
        IF ( fft_method == 'system-specific' )  THEN
 
@@ -346,4 +363,9 @@
           CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4,      &
                        trig_yb, worky, 0 )
+#elif defined( __cuda_fft )
+          CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
+          CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
+          CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
+          CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
 #else
           message_string = 'no system-specific fft-call available'
@@ -425,4 +447,8 @@
 #elif defined( __nec )
        REAL(wp), DIMENSION(6*(nx+1)) ::  work2  !<
+#elif defined( __cuda_fft )
+       COMPLEX(dp), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) ::           &
+          ar_tmp  !<
+       !$ACC DECLARE CREATE(ar_tmp)
 #endif
 
@@ -726,4 +752,53 @@
              ENDDO
              !$OMP END PARALLEL
+
+          ENDIF
+
+#elif defined( __cuda_fft )
+
+          IF ( forward_fft )  THEN
+
+             !$ACC HOST_DATA USE_DEVICE(ar, ar_tmp)
+             CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
+             !$ACC END HOST_DATA
+
+             !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i,j,k) &
+             !$ACC PRESENT(ar, ar_tmp)
+             DO  k = nzb_x, nzt_x
+                DO  j = nys_x, nyn_x
+
+                   DO  i = 0, (nx+1)/2
+                      ar(i,j,k)      = REAL( ar_tmp(i,j,k), KIND=wp )  * dnx
+                   ENDDO
+
+                   DO  i = 1, (nx+1)/2 - 1
+                      ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
+                   ENDDO
+
+                ENDDO
+             ENDDO
+
+          ELSE
+
+             !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i,j,k) &
+             !$ACC PRESENT(ar, ar_tmp)
+             DO  k = nzb_x, nzt_x
+                DO  j = nys_x, nyn_x
+
+                   ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
+
+                   DO  i = 1, (nx+1)/2 - 1
+                      ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k),        &
+                                             KIND=wp )
+                   ENDDO
+                   ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp,     &
+                                                 KIND=wp )
+
+                ENDDO
+             ENDDO
+
+             !$ACC HOST_DATA USE_DEVICE(ar, ar_tmp)
+             CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
+             !$ACC END HOST_DATA
 
           ENDIF
@@ -1001,4 +1076,8 @@
 #elif defined( __nec )
        REAL(wp), DIMENSION(6*(ny+1)) ::  work2  !<
+#elif defined( __cuda_fft )
+       COMPLEX(dp), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) ::           &
+          ar_tmp  !<
+       !$ACC DECLARE CREATE(ar_tmp)
 #endif
 
@@ -1278,4 +1357,53 @@
 
           ENDIF
+#elif defined( __cuda_fft )
+
+          IF ( forward_fft )  THEN
+
+             !$ACC HOST_DATA USE_DEVICE(ar, ar_tmp)
+             CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
+             !$ACC END HOST_DATA
+
+             !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i,j,k) &
+             !$ACC PRESENT(ar, ar_tmp)
+             DO  k = nzb_y, nzt_y
+                DO  i = nxl_y, nxr_y
+
+                   DO  j = 0, (ny+1)/2
+                      ar(j,i,k)      = REAL( ar_tmp(j,i,k), KIND=wp )  * dny
+                   ENDDO
+
+                   DO  j = 1, (ny+1)/2 - 1
+                      ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
+                   ENDDO
+
+                ENDDO
+             ENDDO
+
+          ELSE
+
+             !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i,j,k) &
+             !$ACC PRESENT(ar, ar_tmp)
+             DO  k = nzb_y, nzt_y
+                DO  i = nxl_y, nxr_y
+
+                   ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0_wp, KIND=wp )
+
+                   DO  j = 1, (ny+1)/2 - 1
+                      ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k),        &
+                                             KIND=wp )
+                   ENDDO
+                   ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0_wp,     &
+                                                 KIND=wp )
+
+                ENDDO
+             ENDDO
+
+             !$ACC HOST_DATA USE_DEVICE(ar, ar_tmp)
+             CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
+             !$ACC END HOST_DATA
+
+          ENDIF
+
 #endif
 
Index: palm/trunk/SOURCE/poisfft_mod.f90
===================================================================
--- palm/trunk/SOURCE/poisfft_mod.f90	(revision 3633)
+++ palm/trunk/SOURCE/poisfft_mod.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3241 2018-09-12 15:02:00Z raasch
 ! unused variables removed,
 ! declarations of omp_get_thread_num now as omp-directive
@@ -250,4 +253,5 @@
        REAL(wp), DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar      !<
        REAL(wp), DIMENSION(nys:nyn,nxl:nxr,1:nz) ::  ar_inv  !<
+       !$ACC DECLARE CREATE(ar_inv)
 
        REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  ar1      !<
@@ -262,4 +266,5 @@
        IF ( .NOT. poisfft_initialized )  CALL poisfft_init
 
+#ifndef _OPENACC
 !
 !--    Two-dimensional Fourier Transformation in x- and y-direction.
@@ -295,4 +300,5 @@
 
        ELSEIF ( .NOT. transpose_compute_overlap )  THEN
+#endif
 
 !
@@ -366,4 +372,5 @@
           CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
 
+#ifndef _OPENACC
        ELSE
 
@@ -698,4 +705,5 @@
 
        ENDIF
+#endif
 
        CALL cpu_log( log_point_s(3), 'poisfft', 'stop' )
Index: palm/trunk/SOURCE/pres.f90
===================================================================
--- palm/trunk/SOURCE/pres.f90	(revision 3633)
+++ palm/trunk/SOURCE/pres.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3347 2018-10-15 14:21:08Z suehring
 ! Bugfixes in offline nesting.
 ! Add comment. 
@@ -232,4 +235,17 @@
 
     CALL cpu_log( log_point(8), 'pres', 'start' )
+
+!$ACC DATA COPYOUT(d, tend) &
+!$ACC COPY(u, v, w, p) &
+!$ACC COPYIN(rho_air, rho_air_zw, ddzu, ddzw, wall_flags_0)
+
+!$ACC DATA &
+!$ACC COPYIN(bc_h(0:1)) &
+!$ACC COPYIN(bc_h(0)%i(1:bc_h(0)%ns)) &
+!$ACC COPYIN(bc_h(0)%j(1:bc_h(0)%ns)) &
+!$ACC COPYIN(bc_h(0)%k(1:bc_h(0)%ns)) &
+!$ACC COPYIN(bc_h(1)%i(1:bc_h(1)%ns)) &
+!$ACC COPYIN(bc_h(1)%j(1:bc_h(1)%ns)) &
+!$ACC COPYIN(bc_h(1)%k(1:bc_h(1)%ns))
 
 
@@ -442,4 +458,6 @@
     ELSE
        !$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+       !$ACC PRESENT(d)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -490,4 +508,7 @@
     !$OMP PARALLEL PRIVATE (i,j,k)
     !$OMP DO SCHEDULE( STATIC )
+    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+    !$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_0) &
+    !$ACC PRESENT(d)
     DO  i = nxl, nxr
        DO  j = nys, nyn
@@ -513,4 +534,7 @@
        !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
        !$OMP DO SCHEDULE( STATIC )
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+       !$ACC REDUCTION(+:threadsum) COPY(threadsum) &
+       !$ACC PRESENT(d)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -547,4 +571,6 @@
 !--    z-direction
        !$OMP PARALLEL DO PRIVATE (i,j,k)
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+       !$ACC PRESENT(d, tend)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -565,4 +591,6 @@
 !--       Upward facing
           !$OMP PARALLEL DO PRIVATE( i, j, k )
+          !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+          !$ACC PRESENT(bc_h, tend)
           DO  m = 1, bc_h(0)%ns
              i = bc_h(0)%i(m)
@@ -574,4 +602,6 @@
 !--       Downward facing
           !$OMP PARALLEL DO PRIVATE( i, j, k )
+          !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+          !$ACC PRESENT(bc_h, tend)
           DO  m = 1, bc_h(1)%ns
              i = bc_h(1)%i(m)
@@ -621,4 +651,6 @@
 !--       Dirichlet
           !$OMP PARALLEL DO PRIVATE (i,j,k)
+          !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
+          !$ACC PRESENT(tend)
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
@@ -688,4 +720,6 @@
        !$OMP PARALLEL PRIVATE (i,j,k)
        !$OMP DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+       !$ACC PRESENT(p, tend)
        DO  i = nxl-1, nxr+1
           DO  j = nys-1, nyn+1
@@ -701,4 +735,6 @@
        !$OMP PARALLEL PRIVATE (i,j,k)
        !$OMP DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+       !$ACC PRESENT(p, tend)
        DO  i = nxl-1, nxr+1
           DO  j = nys-1, nyn+1
@@ -733,4 +769,6 @@
     !$OMP PARALLEL PRIVATE (i,j,k)
     !$OMP DO
+    !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) &
+    !$ACC PRESENT(u, v, w, tend, ddzu, wall_flags_0)
     DO  i = nxl, nxr   
        DO  j = nys, nyn
@@ -900,4 +938,7 @@
 #else
        !$OMP DO SCHEDULE( STATIC )
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+       !$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_0) &
+       !$ACC PRESENT(d)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -916,4 +957,7 @@
 !--    Compute possible PE-sum of divergences for flow_statistics
        !$OMP DO SCHEDULE( STATIC )
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+       !$ACC REDUCTION(+:threadsum) COPY(threadsum) &
+       !$ACC PRESENT(d)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -939,4 +983,6 @@
     CALL cpu_log( log_point(8), 'pres', 'stop' )
 
+!$ACC END DATA
+!$ACC END DATA
 
  END SUBROUTINE pres
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 3633)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3589 2018-11-30 15:09:51Z suehring
 ! Move the control parameter "salsa" from salsa_mod to control_parameters 
 ! (M. Kurppa)
@@ -1566,5 +1569,7 @@
     CALL cpu_log( log_point(5), 'u-equation', 'start' )
 
+    !$ACC KERNELS PRESENT(tend)
     tend = 0.0_wp
+    !$ACC END KERNELS
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( ws_scheme_mom )  THEN
@@ -1614,4 +1619,8 @@
 !
 !-- Prognostic equation for u-velocity component
+    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+    !$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
+    !$ACC PRESENT(tsc(2:5)) &
+    !$ACC PRESENT(u_p)
     DO  i = nxlu, nxr
        DO  j = nys, nyn
@@ -1640,4 +1649,6 @@
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( intermediate_timestep_count == 1 )  THEN
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+          !$ACC PRESENT(tend, tu_m)
           DO  i = nxlu, nxr
              DO  j = nys, nyn
@@ -1649,4 +1660,6 @@
        ELSEIF ( intermediate_timestep_count < &
                 intermediate_timestep_count_max )  THEN
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+          !$ACC PRESENT(tend, tu_m)
           DO  i = nxlu, nxr
              DO  j = nys, nyn
@@ -1666,5 +1679,7 @@
     CALL cpu_log( log_point(6), 'v-equation', 'start' )
 
+    !$ACC KERNELS PRESENT(tend)
     tend = 0.0_wp
+    !$ACC END KERNELS
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( ws_scheme_mom )  THEN
@@ -1711,4 +1726,8 @@
 !
 !-- Prognostic equation for v-velocity component
+    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+    !$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
+    !$ACC PRESENT(tsc(2:5)) &
+    !$ACC PRESENT(v_p)
     DO  i = nxl, nxr
        DO  j = nysv, nyn
@@ -1737,4 +1756,6 @@
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( intermediate_timestep_count == 1 )  THEN
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+          !$ACC PRESENT(tend, tv_m)
           DO  i = nxl, nxr
              DO  j = nysv, nyn
@@ -1746,4 +1767,6 @@
        ELSEIF ( intermediate_timestep_count < &
                 intermediate_timestep_count_max )  THEN
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+          !$ACC PRESENT(tend, tv_m)
           DO  i = nxl, nxr
              DO  j = nysv, nyn
@@ -1763,5 +1786,7 @@
     CALL cpu_log( log_point(7), 'w-equation', 'start' )
 
+    !$ACC KERNELS PRESENT(tend)
     tend = 0.0_wp
+    !$ACC END KERNELS
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( ws_scheme_mom )  THEN
@@ -1804,4 +1829,8 @@
 !
 !-- Prognostic equation for w-velocity component
+    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+    !$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
+    !$ACC PRESENT(tsc(2:5)) &
+    !$ACC PRESENT(w_p)
     DO  i = nxl, nxr
        DO  j = nys, nyn
@@ -1821,4 +1850,6 @@
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( intermediate_timestep_count == 1 )  THEN
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+          !$ACC PRESENT(tend, tw_m)
           DO  i = nxl, nxr
              DO  j = nys, nyn
@@ -1830,4 +1861,6 @@
        ELSEIF ( intermediate_timestep_count < &
                 intermediate_timestep_count_max )  THEN
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+          !$ACC PRESENT(tend, tw_m)
           DO  i = nxl, nxr
              DO  j = nys, nyn
@@ -1868,5 +1901,7 @@
 !--    pt-tendency terms with no communication
        IF ( scalar_advec /= 'bc-scheme' )  THEN
+          !$ACC KERNELS PRESENT(tend)
           tend = 0.0_wp
+          !$ACC END KERNELS
           IF ( timestep_scheme(1:5) == 'runge' )  THEN
              IF ( ws_scheme_sca )  THEN
@@ -1926,4 +1961,9 @@
 !
 !--    Prognostic equation for potential temperature
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+       !$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
+       !$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
+       !$ACC PRESENT(tsc(3:5)) &
+       !$ACC PRESENT(pt_p)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -1945,4 +1985,6 @@
        IF ( timestep_scheme(1:5) == 'runge' )  THEN
           IF ( intermediate_timestep_count == 1 )  THEN
+             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+             !$ACC PRESENT(tend, tpt_m)
              DO  i = nxl, nxr
                 DO  j = nys, nyn
@@ -1954,4 +1996,6 @@
           ELSEIF ( intermediate_timestep_count < &
                    intermediate_timestep_count_max )  THEN
+             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+             !$ACC PRESENT(tend, tpt_m)
              DO  i = nxl, nxr
                 DO  j = nys, nyn
Index: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 3633)
+++ palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 3634)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3597 2018-12-04 08:40:18Z maronga
 ! Added routine for calculating near surface air potential temperature (moved
 ! from urban_surface_mod)
@@ -1009,4 +1012,6 @@
        ibit = MERGE( 1, 0, .NOT. downward )
 
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) &
+       !$ACC PRESENT(surf, u, v)
        DO  m = 1, surf%ns
 
@@ -1297,4 +1302,6 @@
              ELSE
                 !$OMP PARALLEL DO PRIVATE( k, z_mo )
+                !$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
+                !$ACC PRESENT(surf, drho_air_zw)
                 DO  m = 1, surf%ns
 
@@ -1319,4 +1326,7 @@
        IF ( TRIM( most_method ) == 'newton' )  THEN
 
+          !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
+          !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
+          !$ACC PRESENT(surf)
           DO  m = 1, surf%ns
 
@@ -1575,4 +1585,6 @@
           IF ( .NOT. downward )  THEN
              !$OMP PARALLEL  DO PRIVATE( z_mo )
+             !$ACC PARALLEL LOOP PRIVATE(z_mo) &
+             !$ACC PRESENT(surf)
              DO  m = 1, surf%ns
 
@@ -1591,4 +1603,6 @@
           ELSE
              !$OMP PARALLEL  DO PRIVATE( z_mo )
+             !$ACC PARALLEL LOOP PRIVATE(z_mo) &
+             !$ACC PRESENT(surf)
              DO  m = 1, surf%ns
 
@@ -1605,4 +1619,6 @@
        ELSE
           !$OMP PARALLEL DO PRIVATE( z_mo )
+          !$ACC PARALLEL LOOP PRIVATE(z_mo) &
+          !$ACC PRESENT(surf)
           DO  m = 1, surf%ns
              z_mo = surf%z_mo(m)
@@ -1624,4 +1640,6 @@
 
        !$OMP PARALLEL DO PRIVATE( i, j, k )
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+       !$ACC PRESENT(surf, pt)
        DO  m = 1, surf%ns 
 
@@ -1630,4 +1648,5 @@
           k   = surf%k(m)
 
+#ifndef _OPENACC
           IF ( bulk_cloud_model ) THEN
              surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
@@ -1637,9 +1656,13 @@
              surf%qv1(m) = q(k,j,i) 
           ELSE
+#endif
              surf%pt1(m) = pt(k,j,i)
+#ifndef _OPENACC
              IF ( humidity )  THEN
                 surf%qv1(m) = q(k,j,i)
              ELSE
+#endif
                 surf%qv1(m) = 0.0_wp
+#ifndef _OPENACC
              ENDIF
           ENDIF 
@@ -1648,4 +1671,5 @@
              surf%vpt1(m) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
           ENDIF
+#endif
           
        ENDDO
@@ -1658,4 +1682,53 @@
 !-- ( only for upward-facing surfs )
     SUBROUTINE calc_pt_surface
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 
+       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements 
+       
+       k_off = surf%koff
+       !$OMP PARALLEL DO PRIVATE( i, j, k )
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+       !$ACC PRESENT(surf, pt)
+       DO  m = 1, surf%ns 
+
+          i   = surf%i(m)            
+          j   = surf%j(m)
+          k   = surf%k(m)
+
+          surf%pt_surface(m) = pt(k+k_off,j,i)
+
+       ENDDO
+
+    END SUBROUTINE calc_pt_surface
+
+!
+!-- Set mixing ratio at surface grid level. ( Only for upward-facing surfs. )
+    SUBROUTINE calc_q_surface
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  k_off   !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 
+       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements 
+       
+       k_off = surf%koff
+       !$OMP PARALLEL DO PRIVATE( i, j, k )
+       DO  m = 1, surf%ns 
+
+          i   = surf%i(m)            
+          j   = surf%j(m)
+          k   = surf%k(m)
+
+          surf%q_surface(m) = q(k+k_off,j,i)
+
+       ENDDO
+
+    END SUBROUTINE calc_q_surface
+    
+!
+!-- Set virtual potential temperature at surface grid level.
+!-- ( only for upward-facing surfs )
+    SUBROUTINE calc_vpt_surface
 
        IMPLICIT NONE
@@ -1672,51 +1745,4 @@
           k   = surf%k(m)
 
-          surf%pt_surface(m) = pt(k+k_off,j,i)
-
-       ENDDO
-
-    END SUBROUTINE calc_pt_surface
-
-!
-!-- Set mixing ratio at surface grid level. ( Only for upward-facing surfs. )
-    SUBROUTINE calc_q_surface
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  k_off   !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 
-       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements 
-       
-       k_off = surf%koff
-       !$OMP PARALLEL DO PRIVATE( i, j, k )
-       DO  m = 1, surf%ns 
-
-          i   = surf%i(m)            
-          j   = surf%j(m)
-          k   = surf%k(m)
-
-          surf%q_surface(m) = q(k+k_off,j,i)
-
-       ENDDO
-
-    END SUBROUTINE calc_q_surface
-    
-!
-!-- Set virtual potential temperature at surface grid level.
-!-- ( only for upward-facing surfs )
-    SUBROUTINE calc_vpt_surface
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  k_off    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 
-       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements 
-       
-       k_off = surf%koff
-       !$OMP PARALLEL DO PRIVATE( i, j, k )
-       DO  m = 1, surf%ns 
-
-          i   = surf%i(m)            
-          j   = surf%j(m)
-          k   = surf%k(m)
-
           surf%vpt_surface(m) = vpt(k+k_off,j,i)
 
@@ -1741,4 +1767,6 @@
 
           !$OMP PARALLEL DO PRIVATE( i, j, k )
+          !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+          !$ACC PRESENT(surf, drho_air_zw)
           DO  m = 1, surf%ns 
 
@@ -2039,4 +2067,6 @@
           IF ( .NOT. downward )  THEN
              !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
+             !$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) &
+             !$ACC PRESENT(surf, u, rho_air_zw)
              DO  m = 1, surf%ns  
    
@@ -2086,4 +2116,6 @@
           IF ( .NOT. downward )  THEN
              !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
+             !$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) &
+             !$ACC PRESENT(surf, v, rho_air_zw)
              DO  m = 1, surf%ns  
                 i = surf%i(m)            
@@ -2392,4 +2424,5 @@
 !-- Integrated stability function for momentum
     FUNCTION psi_m( zeta ) 
+       !$ACC ROUTINE SEQ
        
        USE kinds
@@ -2429,4 +2462,5 @@
 !-- Integrated stability function for heat and moisture
     FUNCTION psi_h( zeta ) 
+       !$ACC ROUTINE SEQ
        
        USE kinds
@@ -2469,4 +2503,5 @@
 !------------------------------------------------------------------------------!
     FUNCTION phi_m( zeta ) 
+       !$ACC ROUTINE SEQ
    
        IMPLICIT NONE 
Index: palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_mod.f90	(revision 3633)
+++ palm/trunk/SOURCE/surface_mod.f90	(revision 3634)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3597 2018-12-04 08:40:18Z maronga
 ! Added pt_2m and renamed t_surf_10cm to pt_10cm. Removed some _eb variables as
 ! they are no longer used.
@@ -608,4 +611,12 @@
     END INTERFACE init_surface_arrays
 
+    INTERFACE enter_surface_arrays
+       MODULE PROCEDURE enter_surface_arrays
+    END INTERFACE
+
+    INTERFACE exit_surface_arrays
+       MODULE PROCEDURE exit_surface_arrays
+    END INTERFACE
+
     INTERFACE surface_rrd_local
        MODULE PROCEDURE surface_rrd_local
@@ -634,6 +645,7 @@
 !-- Public subroutines and functions
     PUBLIC get_topography_top_index, get_topography_top_index_ji, init_bc,     &
-           init_surfaces, init_surface_arrays, surface_rrd_local,              & 
-           surface_restore_elements, surface_wrd_local, surface_last_actions
+           init_surfaces, init_surface_arrays, enter_surface_arrays,           &
+           exit_surface_arrays, surface_rrd_local, surface_restore_elements,   &
+           surface_wrd_local, surface_last_actions
 
 
@@ -1084,4 +1096,92 @@
 
     END SUBROUTINE init_surface_arrays
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Enter horizontal and vertical surfaces.
+!------------------------------------------------------------------------------!
+    SUBROUTINE enter_surface_arrays
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  l     !<
+       
+       !$ACC ENTER DATA &
+       !$ACC COPYIN(surf_def_h(0:2)) &
+       !$ACC COPYIN(surf_def_v(0:3)) &
+       !$ACC COPYIN(surf_lsm_h) &
+       !$ACC COPYIN(surf_lsm_v(0:3)) &
+       !$ACC COPYIN(surf_usm_h) &
+       !$ACC COPYIN(surf_usm_v(0:3))
+
+       ! Copy data in surf_def_h(0:2)
+       DO  l = 0, 1
+          CALL enter_surface_attributes_h(surf_def_h(l))
+       ENDDO
+       CALL enter_surface_attributes_h_top(surf_def_h(2))
+       ! Copy data in surf_def_v(0:3)
+       DO  l = 0, 3
+          CALL enter_surface_attributes_v(surf_def_v(l))
+       ENDDO
+       ! Copy data in surf_lsm_h
+       CALL enter_surface_attributes_h(surf_lsm_h)
+       ! Copy data in surf_lsm_v(0:3)
+       DO  l = 0, 3
+          CALL enter_surface_attributes_v(surf_lsm_v(l))
+       ENDDO
+       ! Copy data in surf_usm_h
+       CALL enter_surface_attributes_h(surf_usm_h)
+       ! Copy data in surf_usm_v(0:3)
+       DO  l = 0, 3
+          CALL enter_surface_attributes_v(surf_usm_v(l))
+       ENDDO
+
+    END SUBROUTINE enter_surface_arrays
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Exit horizontal and vertical surfaces.
+!------------------------------------------------------------------------------!
+    SUBROUTINE exit_surface_arrays
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  l     !<
+       
+       ! Delete data in surf_def_h(0:2)
+       DO  l = 0, 1
+          CALL exit_surface_attributes_h(surf_def_h(l))
+       ENDDO
+       CALL exit_surface_attributes_h(surf_def_h(2))
+       ! Delete data in surf_def_v(0:3)
+       DO  l = 0, 3
+          CALL exit_surface_attributes_v(surf_def_v(l))
+       ENDDO
+       ! Delete data in surf_lsm_h
+       CALL exit_surface_attributes_h(surf_lsm_h)
+       ! Delete data in surf_lsm_v(0:3)
+       DO  l = 0, 3
+          CALL exit_surface_attributes_v(surf_lsm_v(l))
+       ENDDO
+       ! Delete data in surf_usm_h
+       CALL exit_surface_attributes_h(surf_usm_h)
+       ! Delete data in surf_usm_v(0:3)
+       DO  l = 0, 3
+          CALL exit_surface_attributes_v(surf_usm_v(l))
+       ENDDO
+
+       !$ACC EXIT DATA &
+       !$ACC DELETE(surf_def_h(0:2)) &
+       !$ACC DELETE(surf_def_v(0:3)) &
+       !$ACC DELETE(surf_lsm_h) &
+       !$ACC DELETE(surf_lsm_v(0:3)) &
+       !$ACC DELETE(surf_usm_h) &
+       !$ACC DELETE(surf_usm_v(0:3))
+
+    END SUBROUTINE exit_surface_arrays
 
 
@@ -1328,4 +1428,86 @@
 ! Description:
 ! ------------
+!> Exit memory for upward and downward-facing horizontal surface types,
+!> except for top fluxes.
+!------------------------------------------------------------------------------!
+    SUBROUTINE exit_surface_attributes_h( surfaces )
+
+       IMPLICIT NONE
+   
+       TYPE(surf_type) ::  surfaces  !< respective surface type
+   
+       !$ACC EXIT DATA &
+       !$ACC DELETE(surfaces%start_index(nys:nyn,nxl:nxr)) &
+       !$ACC DELETE(surfaces%end_index(nys:nyn,nxl:nxr)) &
+       !$ACC DELETE(surfaces%i(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%j(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%k(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%z_mo(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%uvw_abs(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%z0(1:surfaces%ns)) &
+       !$ACC COPYOUT(surfaces%us(1:surfaces%ns)) &
+       !$ACC COPYOUT(surfaces%ol(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%rib(1:surfaces%ns)) &
+       !$ACC COPYOUT(surfaces%usws(1:surfaces%ns)) &
+       !$ACC COPYOUT(surfaces%vsws(1:surfaces%ns)) &
+       !$ACC COPYOUT(surfaces%ts(1:surfaces%ns)) &
+       !$ACC COPYOUT(surfaces%shf(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%pt_surface(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%pt1(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%qv1(1:surfaces%ns))
+
+       IF ( .NOT. constant_diffusion )  THEN
+          !$ACC EXIT DATA &
+          !$ACC DELETE(surfaces%u_0(1:surfaces%ns)) &
+          !$ACC DELETE(surfaces%v_0(1:surfaces%ns))
+       ENDIF
+   
+    END SUBROUTINE exit_surface_attributes_h
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Enter memory for upward and downward-facing horizontal surface types,
+!> except for top fluxes.
+!------------------------------------------------------------------------------!
+    SUBROUTINE enter_surface_attributes_h( surfaces )
+
+       IMPLICIT NONE
+
+       TYPE(surf_type) ::  surfaces  !< respective surface type
+
+       !$ACC ENTER DATA &
+       !$ACC COPYIN(surfaces%start_index(nys:nyn,nxl:nxr)) &
+       !$ACC COPYIN(surfaces%end_index(nys:nyn,nxl:nxr)) &
+       !$ACC COPYIN(surfaces%i(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%j(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%k(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%z_mo(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%uvw_abs(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%z0(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%us(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%ol(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%rib(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%usws(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%vsws(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%ts(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%shf(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%pt1(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%qv1(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%pt_surface(1:surfaces%ns))
+
+       IF ( .NOT. constant_diffusion )  THEN
+          !$ACC ENTER DATA &
+          !$ACC COPYIN(surfaces%u_0(1:surfaces%ns)) &
+          !$ACC COPYIN(surfaces%v_0(1:surfaces%ns))
+       ENDIF
+
+    END SUBROUTINE enter_surface_attributes_h
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Deallocating memory for model-top fluxes  
 !------------------------------------------------------------------------------!
@@ -1460,4 +1642,64 @@
 
     END SUBROUTINE allocate_surface_attributes_h_top
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Exit memory for model-top fluxes.
+!------------------------------------------------------------------------------!
+    SUBROUTINE exit_surface_attributes_h_top( surfaces )
+
+       IMPLICIT NONE
+   
+       TYPE(surf_type) ::  surfaces  !< respective surface type
+   
+       !$ACC EXIT DATA &
+       !$ACC DELETE(surfaces%start_index(nys:nyn,nxl:nxr)) &
+       !$ACC DELETE(surfaces%end_index(nys:nyn,nxl:nxr)) &
+       !$ACC DELETE(surfaces%i(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%j(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%k(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%usws(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%vsws(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%shf(1:surfaces%ns))
+
+       IF ( .NOT. constant_diffusion )  THEN
+          !$ACC EXIT DATA &
+          !$ACC DELETE(surfaces%u_0(1:surfaces%ns)) &
+          !$ACC DELETE(surfaces%v_0(1:surfaces%ns))
+       ENDIF
+   
+    END SUBROUTINE exit_surface_attributes_h_top
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Enter memory for model-top fluxes.
+!------------------------------------------------------------------------------!
+    SUBROUTINE enter_surface_attributes_h_top( surfaces )
+
+       IMPLICIT NONE
+
+       TYPE(surf_type) ::  surfaces  !< respective surface type
+
+       !$ACC ENTER DATA &
+       !$ACC COPYIN(surfaces%start_index(nys:nyn,nxl:nxr)) &
+       !$ACC COPYIN(surfaces%end_index(nys:nyn,nxl:nxr)) &
+       !$ACC COPYIN(surfaces%i(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%j(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%k(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%usws(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%vsws(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%shf(1:surfaces%ns))
+
+       IF ( .NOT. constant_diffusion )  THEN
+          !$ACC ENTER DATA &
+          !$ACC COPYIN(surfaces%u_0(1:surfaces%ns)) &
+          !$ACC COPYIN(surfaces%v_0(1:surfaces%ns))
+       ENDIF
+
+    END SUBROUTINE enter_surface_attributes_h_top
 
 
@@ -1703,4 +1945,66 @@
 
     END SUBROUTINE allocate_surface_attributes_v
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Exit memory for vertical surface types. 
+!------------------------------------------------------------------------------!
+    SUBROUTINE exit_surface_attributes_v( surfaces )
+
+       IMPLICIT NONE
+
+       TYPE(surf_type) ::  surfaces  !< respective surface type
+
+       !$ACC EXIT DATA &
+       !$ACC DELETE(surfaces%start_index(nys:nyn,nxl:nxr)) &
+       !$ACC DELETE(surfaces%end_index(nys:nyn,nxl:nxr)) &
+       !$ACC DELETE(surfaces%i(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%j(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%k(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%uvw_abs(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%z0(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%rib(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%mom_flux_uv(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%mom_flux_w(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%mom_flux_tke(0:1,1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%ts(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%shf(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%pt1(1:surfaces%ns)) &
+       !$ACC DELETE(surfaces%qv1(1:surfaces%ns))
+
+    END SUBROUTINE exit_surface_attributes_v
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Enter memory for vertical surface types. 
+!------------------------------------------------------------------------------!
+    SUBROUTINE enter_surface_attributes_v( surfaces )
+   
+       IMPLICIT NONE
+   
+       TYPE(surf_type) ::  surfaces  !< respective surface type
+   
+       !$ACC ENTER DATA &
+       !$ACC COPYIN(surfaces%start_index(nys:nyn,nxl:nxr)) &
+       !$ACC COPYIN(surfaces%end_index(nys:nyn,nxl:nxr)) &
+       !$ACC COPYIN(surfaces%i(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%j(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%k(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%uvw_abs(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%z0(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%rib(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%mom_flux_uv(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%mom_flux_w(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%mom_flux_tke(0:1,1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%ts(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%shf(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%pt1(1:surfaces%ns)) &
+       !$ACC COPYIN(surfaces%qv1(1:surfaces%ns))
+   
+    END SUBROUTINE enter_surface_attributes_v
 
 
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 3633)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3597 2018-12-04 08:40:18Z maronga
 ! Removed call to calculation of near air (10 cm) potential temperature (now in
 ! surface layer fluxes)
@@ -421,8 +424,11 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  diss, diss_p, dzu, e, e_p, nc, nc_p, nr, nr_p, prho, pt, pt_p, pt_init, &
-               q_init, q, qc, qc_p, ql, ql_c, ql_v, ql_vp, qr, qr_p, q_p,      &
-               ref_state, rho_ocean, s, s_p, sa_p, tend, u, u_p, v, vpt,       &
-               v_p, w, w_p
+        ONLY:  d, diss, diss_p, ddzu, dd2zu, ddzw, drho_air, drho_air_zw, dzu, &
+               dzw, e, e_p, kh, km, nc, nc_p, nr, nr_p, p, prho, pt, pt_p,     &
+               pt_init, ptdf_x, ptdf_y, q_init, q, qc, qc_p, ql, ql_c, ql_v,   &
+               ql_vp, qr, qr_p, q_p, rdf, rdf_sc, ref_state, rho_air,          &
+               rho_air_zw, rho_ocean, s, s_p, sa_p, te_m, tend, tpt_m, tu_m,   &
+               tv_m, tw_m, u, ug, u_init, u_p, u_stokes_zu, v, vg, v_init, vpt,&
+               v_p, v_stokes_zu, w, w_p, zu
 
     USE biometeorology_mod,                                                    &
@@ -479,5 +485,5 @@
                time_dopr_listing, time_dopts, time_dosp, time_dosp_av,         &
                time_dots, time_do_av, time_do_sla, time_disturb, time_dvrp,    &
-               time_run_control, time_since_reference_point,                   &
+               time_run_control, time_since_reference_point, tsc,              &
                turbulent_inflow, turbulent_outflow, urban_surface,             &
                use_initial_profile_as_reference,                               &
@@ -499,5 +505,6 @@
 
     USE indices,                                                               &
-        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, nzb, nzt
+        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, &
+               nz, nzb_max, advc_flags_1, advc_flags_2, wall_flags_0
 
     USE indoor_model_mod,                                                      &
@@ -556,4 +563,7 @@
               salsa_boundary_conds, salsa_gas, salsa_gases_from_chem,          &
               skip_time_do_salsa                      
+    
+    USE salsa_util_mod,                                                        &
+           ONLY:  sums_salsa_ws_l      
 
     USE spectra_mod,                                                           &
@@ -562,5 +572,10 @@
 
     USE statistics,                                                            &
-        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l
+        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l,                &
+               rmask, statistic_regions, weight_substep, sums_l_l,             &
+               sums_us2_ws_l, sums_wsus_ws_l, sums_vs2_ws_l, sums_wsvs_ws_l,   &
+               sums_ws2_ws_l, sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,&
+               sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsncs_ws_l,              &
+               sums_wsnrs_ws_l, sums_wsss_ws_l
 
     USE surface_layer_fluxes_mod,                                              &
@@ -568,5 +583,6 @@
 
     USE surface_mod,                                                           &
-        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
+        ONLY:  bc_h, surf_def_h, surf_lsm_h, surf_usm_h,                       &
+               enter_surface_arrays, exit_surface_arrays
 
     USE surface_output_mod,                                                    &
@@ -623,4 +639,86 @@
     REAL(wp) ::  time_since_reference_point_save  !< original value of
                                                   !< time_since_reference_point
+
+
+! Copy data from arrays_3d
+!$ACC DATA & 
+!$ACC COPY(d(nzb+1:nzt,nys:nyn,nxl:nxr)) &
+!$ACC COPY(e(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(u(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(v(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(w(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(km(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
+
+!$ACC DATA &
+!$ACC COPY(e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(v_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(w_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(tv_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(tw_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPY(tpt_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
+
+!$ACC DATA &
+!$ACC COPYIN(rho_air(nzb:nzt+1), drho_air(nzb:nzt+1)) &
+!$ACC COPYIN(rho_air_zw(nzb:nzt+1), drho_air_zw(nzb:nzt+1)) &
+!$ACC COPYIN(zu(nzb:nzt+1)) &
+!$ACC COPYIN(dzu(1:nzt+1), dzw(1:nzt+1)) &
+!$ACC COPYIN(ddzu(1:nzt+1), dd2zu(1:nzt)) &
+!$ACC COPYIN(ddzw(1:nzt+1)) &
+!$ACC COPYIN(rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt)) &
+!$ACC COPYIN(ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng)) &
+!$ACC COPYIN(ref_state(0:nz+1)) &
+!$ACC COPYIN(u_init(0:nz+1), v_init(0:nz+1)) &
+!$ACC COPYIN(u_stokes_zu(nzb:nzt+1), v_stokes_zu(nzb:nzt+1)) &
+!$ACC COPYIN(pt_init(0:nz+1)) &
+!$ACC COPYIN(ug(0:nz+1), vg(0:nz+1))
+
+! Copy data from control_parameters
+!$ACC DATA &
+!$ACC COPYIN(tsc(1:5))
+
+! Copy data from indices
+!$ACC DATA &
+!$ACC COPYIN(advc_flags_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPYIN(advc_flags_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+!$ACC COPYIN(wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
+
+! Copy data from surface_mod
+!$ACC DATA &
+!$ACC COPYIN(bc_h(0:1)) &
+!$ACC COPYIN(bc_h(0)%i(1:bc_h(0)%ns)) &
+!$ACC COPYIN(bc_h(0)%j(1:bc_h(0)%ns)) &
+!$ACC COPYIN(bc_h(0)%k(1:bc_h(0)%ns)) &
+!$ACC COPYIN(bc_h(1)%i(1:bc_h(1)%ns)) &
+!$ACC COPYIN(bc_h(1)%j(1:bc_h(1)%ns)) &
+!$ACC COPYIN(bc_h(1)%k(1:bc_h(1)%ns))
+
+! Copy data from statistics
+!$ACC DATA &
+!$ACC COPYIN(hom(nzb+1:nzb_max,1,1:3,0)) &
+!$ACC COPYIN(rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions)) &
+!$ACC COPYIN(weight_substep(1:intermediate_timestep_count_max)) &
+!$ACC COPY(sums_l_l(nzb:nzt+1,0:statistic_regions,0)) &
+!$ACC COPY(sums_us2_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_wsus_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_vs2_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_wsvs_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_ws2_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_wspts_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_wssas_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_wsqs_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_wsqcs_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_wsqrs_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_wsncs_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_wsnrs_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_wsss_ws_l(nzb:nzt+1,0)) &
+!$ACC COPY(sums_salsa_ws_l(nzb:nzt+1,0))
 
 !
@@ -739,4 +837,8 @@
        ENDIF   
        
+#ifdef _OPENACC
+       CALL enter_surface_arrays
+#endif
+        
 !
 !--    Start of intermediate step loop
@@ -1252,4 +1354,29 @@
 
        ENDDO   ! Intermediate step loop
+
+!
+!--    Will be used at some point by flow_statistics.
+       !$ACC UPDATE HOST(e, u, v, w, pt) &
+       !$ACC HOST(kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+       !$ACC HOST(km(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+       !$ACC HOST(sums_l_l(nzb:nzt+1,0:statistic_regions,0)) &
+       !$ACC HOST(sums_us2_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_wsus_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_vs2_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_wsvs_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_ws2_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_wspts_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_wssas_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_wsqs_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_wsqcs_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_wsqrs_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_wsncs_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_wsnrs_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_wsss_ws_l(nzb:nzt+1,0)) &
+       !$ACC HOST(sums_salsa_ws_l(nzb:nzt+1,0))
+#ifdef _OPENACC
+       CALL exit_surface_arrays
+#endif
+
 !
 !--    If required, consider chemical emissions
@@ -1655,4 +1782,12 @@
     ENDDO   ! time loop
 
+!$ACC END DATA
+!$ACC END DATA
+!$ACC END DATA
+!$ACC END DATA
+!$ACC END DATA
+!$ACC END DATA
+!$ACC END DATA
+
 !
 !-- Vertical nesting: Deallocate variables initialized for vertical nesting    
Index: palm/trunk/SOURCE/timestep.f90
===================================================================
--- palm/trunk/SOURCE/timestep.f90	(revision 3633)
+++ palm/trunk/SOURCE/timestep.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3311 2018-10-05 12:34:56Z raasch
 ! Stokes drift is regarded in timestep calculation
 ! 
@@ -204,4 +207,5 @@
     REAL(wp), DIMENSION(3)         ::  reduce_l    !< 
     REAL(wp), DIMENSION(nzb+1:nzt) ::  dxyz2_min   !<  
+    !$ACC DECLARE CREATE(dxyz2_min)
 
 
@@ -277,4 +281,8 @@
        dt_v_l = 999999.9_wp
        dt_w_l = 999999.9_wp
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+       !$ACC COPY(dt_u_l, dt_v_l, dt_w_l, u_stokes_zu, v_stokes_zu) &
+       !$ACC REDUCTION(MIN: dt_u_l, dt_v_l, dt_w_l) &
+       !$ACC PRESENT(u, v, w, dzu)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -313,4 +321,5 @@
        dt_diff_l = 999999.0_wp
 
+       !$ACC PARALLEL LOOP PRESENT(dxyz2_min, dzw)
        DO  k = nzb+1, nzt
            dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125_wp
@@ -319,4 +328,7 @@
        !$OMP PARALLEL private(i,j,k) reduction(MIN: dt_diff_l)
        !$OMP DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+       !$ACC COPY(dt_diff_l) REDUCTION(MIN: dt_diff_l) &
+       !$ACC PRESENT(dxyz2_min, kh, km)
        DO  i = nxl, nxr
           DO  j = nys, nyn
Index: palm/trunk/SOURCE/timestep_scheme_steering.f90
===================================================================
--- palm/trunk/SOURCE/timestep_scheme_steering.f90	(revision 3633)
+++ palm/trunk/SOURCE/timestep_scheme_steering.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 2718 2018-01-02 08:49:38Z maronga
 ! Corrected "Former revisions" section
 ! 
@@ -90,4 +93,6 @@
        ENDIF
 
+       !$ACC UPDATE DEVICE(tsc(1:5))
+
     ELSEIF ( timestep_scheme == 'euler' )  THEN
 !
Index: palm/trunk/SOURCE/transpose.f90
===================================================================
--- palm/trunk/SOURCE/transpose.f90	(revision 3633)
+++ palm/trunk/SOURCE/transpose.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3241 2018-09-12 15:02:00Z raasch
 ! unused variables removed
 ! 
@@ -119,4 +122,6 @@
     !$OMP  PARALLEL PRIVATE ( i, j, k )
     !$OMP  DO
+     !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+     !$ACC PRESENT(f_inv, f_in)
      DO  i = 0, nx
          DO  k = nzb_x, nzt_x
@@ -166,4 +171,5 @@
 
     REAL(wp), DIMENSION(nyn_x-nys_x+1,nzb_y:nzt_y,nxl_y:nxr_y,0:pdims(2)-1) ::  work !< 
+    !$ACC DECLARE CREATE(work)
 
 
@@ -174,8 +180,10 @@
 !--    Transpose array
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start', cpu_log_nowait )
+       !$ACC UPDATE HOST(f_inv)
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
        CALL MPI_ALLTOALL( f_inv(nys_x,nzb_x,0),  sendrecvcount_xy, MPI_REAL, &
                           work(1,nzb_y,nxl_y,0), sendrecvcount_xy, MPI_REAL, &
                           comm1dy, ierr )
+       !$ACC UPDATE DEVICE(work)
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
 
@@ -186,4 +194,6 @@
        DO  l = 0, pdims(2) - 1
           ys = 0 + l * ( nyn_x - nys_x + 1 )
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+          !$ACC PRESENT(f_out, work)
           DO  i = nxl_y, nxr_y
              DO  k = nzb_y, nzt_y
@@ -203,4 +213,6 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+       !$ACC PRESENT(f_out, f_inv)
        DO  k = nzb_y, nzt_y
           DO  i = nxl_y, nxr_y
@@ -246,4 +258,6 @@
     !$OMP  PARALLEL PRIVATE ( i, j, k )
     !$OMP  DO
+     !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+     !$ACC PRESENT(f_out, f_inv)
      DO  k = 1, nz
          DO  i = nxl, nxr
@@ -293,4 +307,5 @@
 
     REAL(wp), DIMENSION(nys_x:nyn_x,nnx,nzb_x:nzt_x,0:pdims(1)-1) ::  work !< 
+    !$ACC DECLARE CREATE(work)
 
 
@@ -307,4 +322,6 @@
        DO  l = 0, pdims(1) - 1
           xs = 0 + l * nnx
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+          !$ACC PRESENT(work, f_in)
           DO  k = nzb_x, nzt_x
              DO  i = xs, xs + nnx - 1
@@ -320,8 +337,10 @@
 !--    Transpose array
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start', cpu_log_nowait )
+       !$ACC UPDATE HOST(work)
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
        CALL MPI_ALLTOALL( work(nys_x,1,nzb_x,0), sendrecvcount_zx, MPI_REAL, &
                           f_inv(nys,nxl,1),      sendrecvcount_zx, MPI_REAL, &
                           comm1dx, ierr )
+       !$ACC UPDATE DEVICE(f_inv)
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
 #endif
@@ -333,4 +352,6 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+       !$ACC PRESENT(f_inv, f_in)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -378,4 +399,6 @@
     !$OMP  PARALLEL PRIVATE ( i, j, k )
     !$OMP  DO
+     !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+     !$ACC PRESENT(f_out, f_inv)
      DO  i = 0, nx
          DO  k = nzb_x, nzt_x
@@ -425,4 +448,5 @@
 
     REAL(wp), DIMENSION(nyn_x-nys_x+1,nzb_y:nzt_y,nxl_y:nxr_y,0:pdims(2)-1) ::  work !< 
+    !$ACC DECLARE CREATE(work)
 
 
@@ -436,4 +460,6 @@
        DO  l = 0, pdims(2) - 1
           ys = 0 + l * ( nyn_x - nys_x + 1 )
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+          !$ACC PRESENT(work, f_in)
           DO  i = nxl_y, nxr_y
              DO  k = nzb_y, nzt_y
@@ -449,8 +475,10 @@
 !--    Transpose array
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start', cpu_log_nowait )
+       !$ACC UPDATE HOST(work)
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
        CALL MPI_ALLTOALL( work(1,nzb_y,nxl_y,0), sendrecvcount_xy, MPI_REAL, &
                           f_inv(nys_x,nzb_x,0),  sendrecvcount_xy, MPI_REAL, &
                           comm1dy, ierr )
+       !$ACC UPDATE DEVICE(f_inv)
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
 #endif
@@ -462,4 +490,6 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+       !$ACC PRESENT(f_inv, f_in)
        DO  i = nxl_y, nxr_y
           DO  k = nzb_y, nzt_y
@@ -587,4 +617,6 @@
     !$OMP  PARALLEL PRIVATE ( i, j, k )
     !$OMP  DO
+     !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+     !$ACC PRESENT(f_inv, f_in)
      DO  j = 0, ny
          DO  k = nzb_y, nzt_y
@@ -634,4 +666,5 @@
 
     REAL(wp), DIMENSION(nxl_z:nxr_z,nzt_y-nzb_y+1,nys_z:nyn_z,0:pdims(1)-1) ::  work !< 
+    !$ACC DECLARE CREATE(work)
 
 
@@ -643,4 +676,6 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+       !$ACC PRESENT(f_out, f_inv)
        DO  j = 0, ny
           DO  k = nzb_y, nzt_y
@@ -658,8 +693,10 @@
 !--    Transpose array
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start', cpu_log_nowait )
+       !$ACC UPDATE HOST(f_inv)
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
        CALL MPI_ALLTOALL( f_inv(nxl_y,nzb_y,0),  sendrecvcount_yz, MPI_REAL, &
                           work(nxl_z,1,nys_z,0), sendrecvcount_yz, MPI_REAL, &
                           comm1dx, ierr )
+       !$ACC UPDATE DEVICE(work)
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
 
@@ -670,4 +707,6 @@
        DO  l = 0, pdims(1) - 1
           zs = 1 + l * ( nzt_y - nzb_y + 1 )
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+          !$ACC PRESENT(f_out, work)
           DO  j = nys_z, nyn_z
              DO  k = zs, zs + nzt_y - nzb_y
@@ -714,4 +753,6 @@
     !$OMP  PARALLEL PRIVATE ( i, j, k )
     !$OMP  DO
+    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+    !$ACC PRESENT(f_in, f_inv)
      DO  k = 1,nz
          DO  i = nxl, nxr
@@ -761,4 +802,5 @@
 
     REAL(wp), DIMENSION(nys_x:nyn_x,nnx,nzb_x:nzt_x,0:pdims(1)-1) ::  work !< 
+    !$ACC DECLARE CREATE(work)
 
 
@@ -770,4 +812,6 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+       !$ACC PRESENT(f_out, f_inv)
        DO  k = 1, nz
           DO  i = nxl, nxr
@@ -785,8 +829,10 @@
 !--    Transpose array
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start', cpu_log_nowait )
+       !$ACC UPDATE HOST(f_inv)
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
        CALL MPI_ALLTOALL( f_inv(nys,nxl,1),      sendrecvcount_zx, MPI_REAL, &
                           work(nys_x,1,nzb_x,0), sendrecvcount_zx, MPI_REAL, &
                           comm1dx, ierr )
+       !$ACC UPDATE DEVICE(work)
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
 
@@ -797,4 +843,6 @@
        DO  l = 0, pdims(1) - 1
           xs = 0 + l * nnx
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+          !$ACC PRESENT(f_out, work)
           DO  k = nzb_x, nzt_x
              DO  i = xs, xs + nnx - 1
@@ -845,4 +893,6 @@
     !$OMP  PARALLEL PRIVATE ( i, j, k )
     !$OMP  DO
+    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+    !$ACC PRESENT(f_out, f_inv)
      DO  k = nzb_y, nzt_y
          DO  j = 0, ny
@@ -892,4 +942,5 @@
 
     REAL(wp), DIMENSION(nxl_z:nxr_z,nzt_y-nzb_y+1,nys_z:nyn_z,0:pdims(1)-1) ::  work !<
+    !$ACC DECLARE CREATE(work)
 
 !
@@ -905,4 +956,6 @@
        DO  l = 0, pdims(1) - 1
           zs = 1 + l * ( nzt_y - nzb_y + 1 )
+          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+          !$ACC PRESENT(work, f_in)
           DO  j = nys_z, nyn_z
              DO  k = zs, zs + nzt_y - nzb_y
@@ -918,8 +971,10 @@
 !--    Transpose array
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start', cpu_log_nowait )
+       !$ACC UPDATE HOST(work)
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
        CALL MPI_ALLTOALL( work(nxl_z,1,nys_z,0), sendrecvcount_yz, MPI_REAL, &
                           f_inv(nxl_y,nzb_y,0),  sendrecvcount_yz, MPI_REAL, &
                           comm1dx, ierr )
+       !$ACC UPDATE DEVICE(f_inv)
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
 #endif
@@ -930,4 +985,6 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+       !$ACC PRESENT(f_inv, f_in)
        DO  k = nzb_y, nzt_y
           DO  j = 0, ny
Index: palm/trunk/SOURCE/tridia_solver_mod.f90
===================================================================
--- palm/trunk/SOURCE/tridia_solver_mod.f90	(revision 3633)
+++ palm/trunk/SOURCE/tridia_solver_mod.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -150,5 +153,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu_pres, ddzw, rho_air_zw
+           ONLY:  ddzu_pres, ddzw, rho_air_zw, tri
 
        IMPLICIT NONE
@@ -169,4 +172,8 @@
        CALL maketri
        CALL split
+
+       !$ACC ENTER DATA &
+       !$ACC COPYIN(ddzuw(0:nz-1,1:3)) &
+       !$ACC COPYIN(tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,1:2))
 
     END SUBROUTINE tridia_init
@@ -290,8 +297,11 @@
 
           REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1)   ::  ar1 !< 
+          !$ACC DECLARE CREATE(ar1)
 
 !
 !--       Forward substitution
+          !$ACC PARALLEL PRESENT(ar, ar1, tri) PRIVATE(i,j,k)
           DO  k = 0, nz - 1
+             !$ACC LOOP COLLAPSE(2)
              DO  j = nys_z, nyn_z
                 DO  i = nxl_z, nxr_z
@@ -306,4 +316,5 @@
              ENDDO
           ENDDO
+          !$ACC END PARALLEL
 
 !
@@ -312,5 +323,7 @@
 !--       by zero appearing if the pressure bc is set to neumann at the top of
 !--       the model domain.
+          !$ACC PARALLEL PRESENT(ar, ar1, ddzuw, tri) PRIVATE(i,j,k)
           DO  k = nz-1, 0, -1
+             !$ACC LOOP COLLAPSE(2)
              DO  j = nys_z, nyn_z
                 DO  i = nxl_z, nxr_z
@@ -325,4 +338,5 @@
              ENDDO
           ENDDO
+          !$ACC END PARALLEL
 
 !
@@ -332,4 +346,5 @@
           IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
              IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
+                !$ACC PARALLEL LOOP PRESENT(ar)
                 DO  k = 1, nz
                    ar(nxl_z,nys_z,k) = 0.0_wp
Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3633)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3634)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! OpenACC port for SPEC
+! 
+! 3550 2018-11-21 16:01:01Z gronemeier
 ! - calculate diss production same in vector and cache optimization
 ! - move boundary condition for e and diss to boundary_conds
@@ -1619,4 +1622,6 @@
        ENDDO  !k loop
 
+       !$ACC ENTER DATA COPYIN(l_black(nzb:nzt+1))
+
     ENDIF  !LES or RANS mode
 
@@ -1624,4 +1629,7 @@
 !-- Set lateral boundary conditions for l_wall
     CALL exchange_horiz( l_wall, nbgp )
+
+    !$ACC ENTER DATA COPYIN(l_grid(nzb:nzt+1)) &
+    !$ACC COPYIN(l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
 
     CONTAINS
@@ -1777,4 +1785,6 @@
 !--    Default surfaces, upward-facing
        !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
+       !$ACC PRESENT(surf_def_h(0), u, v, drho_air_zw, zu)
        DO  m = 1, surf_def_h(0)%ns
 
@@ -1811,4 +1821,6 @@
 !--    Default surfaces, downward-facing surfaces
        !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
+       !$ACC PRESENT(surf_def_h(1), u, v, drho_air_zw, zu, km)
        DO  m = 1, surf_def_h(1)%ns
 
@@ -1842,4 +1854,6 @@
 !--    Natural surfaces, upward-facing
        !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
+       !$ACC PRESENT(surf_lsm_h, u, v, drho_air_zw, zu)
        DO  m = 1, surf_lsm_h%ns
 
@@ -1876,4 +1890,6 @@
 !--    Urban surfaces, upward-facing
        !$OMP PARALLEL DO PRIVATE(i,j,k,m)
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
+       !$ACC PRESENT(surf_usm_h, u, v, drho_air_zw, zu)
        DO  m = 1, surf_usm_h%ns
 
@@ -1970,5 +1986,7 @@
              CALL advec_s_up( e )
           ELSE
+             !$ACC KERNELS PRESENT(tend)
              tend = 0.0_wp
+             !$ACC END KERNELS
              IF ( timestep_scheme(1:5) == 'runge' )  THEN
                 IF ( ws_scheme_sca )  THEN
@@ -2006,4 +2024,8 @@
 !--    reasons in the course of the integration. In such cases the old TKE
 !--    value is reduced by 90%.
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+       !$ACC PRESENT(e, tend, te_m, wall_flags_0) &
+       !$ACC PRESENT(tsc(3:3)) &
+       !$ACC PRESENT(e_p)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -2024,4 +2046,6 @@
        IF ( timestep_scheme(1:5) == 'runge' )  THEN
           IF ( intermediate_timestep_count == 1 )  THEN
+             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+             !$ACC PRESENT(tend, te_m)
              DO  i = nxl, nxr
                 DO  j = nys, nyn
@@ -2033,4 +2057,6 @@
           ELSEIF ( intermediate_timestep_count < &
                    intermediate_timestep_count_max )  THEN
+             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+             !$ACC PRESENT(tend, te_m)
              DO  i = nxl, nxr
                 DO  j = nys, nyn
@@ -2380,6 +2406,15 @@
 !-- points first, gradients at surface-bounded grid points will be
 !-- overwritten further below.
+    !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, l) &
+    !$ACC PRIVATE(surf_s, surf_e) &
+    !$ACC PRIVATE(dudx(:), dudy(:), dudz(:), dvdx(:), dvdy(:), dvdz(:), dwdx(:), dwdy(:), dwdz(:)) &
+    !$ACC PRESENT(e, u, v, w, diss, dd2zu, ddzw, km, wall_flags_0) &
+    !$ACC PRESENT(tend) &
+    !$ACC PRESENT(surf_def_h(0:1), surf_def_v(0:3)) &
+    !$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) &
+    !$ACC PRESENT(surf_usm_h, surf_usm_v(0:3))
     DO  i = nxl, nxr
        DO  j = nys, nyn
+          !$ACC LOOP PRIVATE(k)
           DO  k = nzb+1, nzt
 
@@ -2423,4 +2458,5 @@
                 surf_s = surf_def_v(l)%start_index(j,i)
                 surf_e = surf_def_v(l)%end_index(j,i)
+                !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
                 DO  m = surf_s, surf_e
                    k           = surf_def_v(l)%k(m)
@@ -2441,4 +2477,5 @@
                 surf_s = surf_lsm_v(l)%start_index(j,i)
                 surf_e = surf_lsm_v(l)%end_index(j,i)
+                !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
                 DO  m = surf_s, surf_e
                    k           = surf_lsm_v(l)%k(m)
@@ -2459,4 +2496,5 @@
                 surf_s = surf_usm_v(l)%start_index(j,i)
                 surf_e = surf_usm_v(l)%end_index(j,i)
+                !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
                 DO  m = surf_s, surf_e
                    k           = surf_usm_v(l)%k(m)
@@ -2479,4 +2517,5 @@
                 surf_s = surf_def_v(l)%start_index(j,i)
                 surf_e = surf_def_v(l)%end_index(j,i)
+                !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
                 DO  m = surf_s, surf_e
                    k     = surf_def_v(l)%k(m)
@@ -2497,4 +2536,5 @@
                 surf_s = surf_lsm_v(l)%start_index(j,i)
                 surf_e = surf_lsm_v(l)%end_index(j,i)
+                !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
                 DO  m = surf_s, surf_e
                    k     = surf_lsm_v(l)%k(m)
@@ -2515,4 +2555,5 @@
                 surf_s = surf_usm_v(l)%start_index(j,i)
                 surf_e = surf_usm_v(l)%end_index(j,i)
+                !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
                 DO  m = surf_s, surf_e
                    k     = surf_usm_v(l)%k(m)
@@ -2534,4 +2575,5 @@
              surf_s = surf_def_h(0)%start_index(j,i)
              surf_e = surf_def_h(0)%end_index(j,i)
+             !$ACC LOOP PRIVATE(m, k)
              DO  m = surf_s, surf_e
                 k = surf_def_h(0)%k(m)
@@ -2550,4 +2592,5 @@
              surf_s = surf_lsm_h%start_index(j,i)
              surf_e = surf_lsm_h%end_index(j,i)
+             !$ACC LOOP PRIVATE(m, k)
              DO  m = surf_s, surf_e
                 k = surf_lsm_h%k(m)
@@ -2561,4 +2604,5 @@
              surf_s = surf_usm_h%start_index(j,i)
              surf_e = surf_usm_h%end_index(j,i)
+             !$ACC LOOP PRIVATE(m, k)
              DO  m = surf_s, surf_e
                 k = surf_usm_h%k(m)
@@ -2573,4 +2617,5 @@
              surf_s = surf_def_h(1)%start_index(j,i)
              surf_e = surf_def_h(1)%end_index(j,i)
+             !$ACC LOOP PRIVATE(m, k)
              DO  m = surf_s, surf_e
                 k = surf_def_h(1)%k(m)
@@ -2585,4 +2630,5 @@
 
 
+          !$ACC LOOP PRIVATE(k, def, flag)
           DO  k = nzb+1, nzt
 
@@ -2684,7 +2730,16 @@
           ELSE ! or IF ( .NOT. ocean_mode )  THEN
 
+             !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
+             !$ACC PRIVATE(surf_s, surf_e) &
+             !$ACC PRIVATE(tmp_flux(nzb+1:nzt)) &
+             !$ACC PRESENT(e, diss, kh, pt, dd2zu, drho_air_zw, wall_flags_0) &
+             !$ACC PRESENT(tend) &
+             !$ACC PRESENT(surf_def_h(0:2)) &
+             !$ACC PRESENT(surf_lsm_h) &
+             !$ACC PRESENT(surf_usm_h)
              DO  i = nxl, nxr
                 DO  j = nys, nyn
 
+                   !$ACC LOOP PRIVATE(k)
                    DO  k = nzb+1, nzt
                       tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
@@ -2697,4 +2752,5 @@
                          surf_s = surf_def_h(l)%start_index(j,i)
                          surf_e = surf_def_h(l)%end_index(j,i)
+                         !$ACC LOOP PRIVATE(m, k)
                          DO  m = surf_s, surf_e
                             k = surf_def_h(l)%k(m)
@@ -2706,4 +2762,5 @@
                       surf_s = surf_lsm_h%start_index(j,i)
                       surf_e = surf_lsm_h%end_index(j,i)
+                      !$ACC LOOP PRIVATE(m, k)
                       DO  m = surf_s, surf_e
                          k = surf_lsm_h%k(m)
@@ -2714,4 +2771,5 @@
                       surf_s = surf_usm_h%start_index(j,i)
                       surf_e = surf_usm_h%end_index(j,i)
+                      !$ACC LOOP PRIVATE(m, k)
                       DO  m = surf_s, surf_e
                          k = surf_usm_h%k(m)
@@ -2723,4 +2781,5 @@
                       surf_s = surf_def_h(2)%start_index(j,i)
                       surf_e = surf_def_h(2)%end_index(j,i)
+                      !$ACC LOOP PRIVATE(m, k)
                       DO  m = surf_s, surf_e
                          k = surf_def_h(2)%k(m)
@@ -2732,4 +2791,5 @@
 
 !--                   Compute tendency for TKE-production from shear
+                     !$ACC LOOP PRIVATE(k, flag)
                       DO  k = nzb+1, nzt
                          flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
@@ -3628,5 +3688,9 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  ddzu, ddzw, drho_air, rho_air_zw
+        ONLY:  ddzu, dd2zu, ddzw, drho_air, rho_air_zw
+
+    USE control_parameters,                                                    &
+        ONLY:  atmos_ocean_sign, use_single_reference_value,                   &
+               wall_adjustment, wall_adjustment_factor
 
     USE grid_variables,                                                        &
@@ -3653,4 +3717,17 @@
     REAL(wp)     ::  ll             !< adjusted l
     REAL(wp)     ::  var_reference  !< reference temperature
+#ifdef _OPENACC
+!
+!-- From mixing_length_les:
+    REAL(wp)     :: l_stable        !< mixing length according to stratification
+!
+!-- From mixing_length_rans:
+    REAL(wp)     :: duv2_dz2        !< squared vertical gradient of wind vector
+    REAL(wp)     :: l_diss          !< mixing length for dissipation
+    REAL(wp)     :: rif             !< Richardson flux number
+!
+!-- From both:
+    REAL(wp)     :: dvar_dz         !< vertical gradient of var
+#endif
 
 #if defined( __nopointer )
@@ -3659,9 +3736,16 @@
     REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
 #endif
-    REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation  !< TKE dissipation
+    REAL(wp)     ::  dissipation  !< TKE dissipation
 
 
 !
 !-- Calculate the tendency terms
+    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
+    !$ACC PRIVATE(flag, l, ll, dissipation) &
+    !$ACC PRIVATE(l_stable, duv2_dz2, l_diss, rif, dvar_dz) &
+    !$ACC PRESENT(e, u, v, km, var, wall_flags_0) &
+    !$ACC PRESENT(ddzu, dd2zu, ddzw, rho_air_zw, drho_air) &
+    !$ACC PRESENT(l_black, l_grid, l_wall) &
+    !$ACC PRESENT(diss, tend)
     DO  i = nxl, nxr
        DO  j = nys, nyn
@@ -3675,20 +3759,80 @@
              IF ( les_dynamic .OR. les_mw )  THEN
 
+#ifdef _OPENACC
+!
+!--             Cannot call subroutine mixing_length_les because after adding all required
+!--             OpenACC directives the execution crashed reliably when inside the called
+!--             subroutine. I'm not sure why that is...
+                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                IF ( dvar_dz > 0.0_wp ) THEN
+                   IF ( use_single_reference_value )  THEN
+                      l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
+                                         / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
+                   ELSE
+                      l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
+                                         / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
+                   ENDIF
+                ELSE
+                   l_stable = l_grid(k)
+                ENDIF
+!
+!--             Adjustment of the mixing length
+                IF ( wall_adjustment )  THEN
+                   l  = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
+                   ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
+                ELSE
+                   l  = MIN( l_grid(k), l_stable )
+                   ll = l_grid(k)
+                ENDIF
+#else
                 CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
-
-                dissipation(k,j) = ( 0.19_wp + 0.74_wp * l / ll )              &
-                                   * e(k,j,i) * SQRT( e(k,j,i) ) / l
+#endif
+
+                dissipation = ( 0.19_wp + 0.74_wp * l / ll )                   &
+                              * e(k,j,i) * SQRT( e(k,j,i) ) / l
 
              ELSEIF ( rans_tke_l )  THEN
 
+#ifdef _OPENACC
+!
+!--             Same reason as above...
+                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+
+                duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2                 &
+                           + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2                 &
+                           + 1E-30_wp
+
+                IF ( use_single_reference_value )  THEN
+                   rif = g / var_reference * dvar_dz / duv2_dz2
+                ELSE
+                   rif = g / var(k,j,i) * dvar_dz / duv2_dz2
+                ENDIF
+
+                rif = MAX( rif, -5.0_wp )
+                rif = MIN( rif,  1.0_wp )
+
+!
+!--             Calculate diabatic mixing length using Dyer-profile functions
+                IF ( rif >= 0.0_wp )  THEN
+                   l      = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif ), l_wall(k,j,i) )
+                   l_diss = l
+                ELSE
+!
+!--                In case of unstable stratification, use mixing length of neutral case
+!--                for l, but consider profile functions for l_diss
+                   l      = l_wall(k,j,i)
+                   l_diss = l * SQRT( 1.0_wp - 16.0_wp * rif )
+                ENDIF
+#else
                 CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
-
-                dissipation(k,j) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
-
-                diss(k,j,i) = dissipation(k,j) * flag
+#endif
+
+                dissipation = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
+
+                diss(k,j,i) = dissipation * flag
 
              ELSEIF ( rans_tke_e )  THEN
 
-                dissipation(k,j) = diss(k,j,i)
+                dissipation = diss(k,j,i)
 
              ENDIF
@@ -3710,22 +3854,15 @@
                                            ) * ddzw(k) * drho_air(k)           &
                                          ) * flag * dsig_e                     &
-                          - dissipation(k,j) * flag
+                          - dissipation * flag
+!
+!--          Store dissipation if needed for calculating the sgs particle
+!--          velocities
+             IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles  .OR.         &
+                  wang_kernel  .OR.  collision_turbulence  ) )  THEN
+                diss(k,j,i) = dissipation * flag
+             ENDIF
 
           ENDDO
        ENDDO
-
-!
-!--    Store dissipation if needed for calculating the sgs particle
-!--    velocities
-       IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles  .OR.               &
-            wang_kernel  .OR.  collision_turbulence  ) )  THEN
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-                diss(k,j,i) = dissipation(k,j) * MERGE( 1.0_wp, 0.0_wp,        &
-                                               BTEST( wall_flags_0(k,j,i), 0 ) )
-             ENDDO
-          ENDDO
-       ENDIF
-
     ENDDO
 
@@ -4257,4 +4394,6 @@
 !-- Downward facing surfaces
     !$OMP PARALLEL DO PRIVATE(i,j,k)
+    !$ACC PARALLEL LOOP PRIVATE(i,j,k) &
+    !$ACC PRESENT(bc_h(1), kh, km)
     DO  m = 1, bc_h(1)%ns
        i = bc_h(1)%i(m)
@@ -4267,4 +4406,6 @@
 !-- Downward facing surfaces
     !$OMP PARALLEL DO PRIVATE(i,j,k)
+    !$ACC PARALLEL LOOP PRIVATE(i,j,k) &
+    !$ACC PRESENT(bc_h(0), kh, km)
     DO  m = 1, bc_h(0)%ns
        i = bc_h(0)%i(m)
@@ -4277,4 +4418,6 @@
 !-- Model top
     !$OMP PARALLEL DO
+    !$ACC PARALLEL LOOP COLLAPSE(2) &
+    !$ACC PRESENT(kh, km)
     DO  i = nxlg, nxrg
        DO  j = nysg, nyng
@@ -4315,4 +4458,11 @@
  SUBROUTINE tcm_diffusivities_default( var, var_reference )
  
+    USE arrays_3d,                                                             &
+        ONLY:  dd2zu
+
+    USE control_parameters,                                                    &
+        ONLY:  atmos_ocean_sign, use_single_reference_value,                   &
+               wall_adjustment, wall_adjustment_factor
+
     USE statistics,                                                            &
         ONLY :  rmask, sums_l_l
@@ -4331,4 +4481,10 @@
     REAL(wp)     ::  ll                  !< adjusted mixing length
     REAL(wp)     ::  var_reference       !< reference temperature
+#ifdef _OPENACC
+!
+!-- From mixing_length_les:
+    REAL(wp)     :: l_stable        !< mixing length according to stratification
+    REAL(wp)     :: dvar_dz         !< vertical gradient of var
+#endif
 
 #if defined( __nopointer )
@@ -4344,5 +4500,7 @@
 !
 !-- Initialization for calculation of the mixing length profile
+    !$ACC KERNELS PRESENT(sums_l_l)
     sums_l_l = 0.0_wp
+    !$ACC END KERNELS
 
 !
@@ -4353,4 +4511,8 @@
     IF ( les_dynamic .OR. les_mw )  THEN
        !$OMP DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(sr) &
+       !$ACC PRIVATE(flag, dvar_dz, l_stable, l, ll) &
+       !$ACC PRESENT(wall_flags_0, var, dd2zu, e, l_wall, l_grid, rmask) &
+       !$ACC PRESENT(kh, km, sums_l_l)
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -4361,5 +4523,34 @@
 !
 !--             Determine the mixing length for LES closure
+#ifdef _OPENACC
+!
+!--             Cannot call subroutine mixing_length_les because after adding all required
+!--             OpenACC directives the execution crashed reliably when inside the called
+!--             subroutine. I'm not sure why that is...
+                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                IF ( dvar_dz > 0.0_wp ) THEN
+                   IF ( use_single_reference_value )  THEN
+                      l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
+                                         / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
+                   ELSE
+                      l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
+                                         / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
+                   ENDIF
+                ELSE
+                   l_stable = l_grid(k)
+                ENDIF
+!
+!--             Adjustment of the mixing length
+                IF ( wall_adjustment )  THEN
+                   l  = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
+                   ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
+                ELSE
+                   l  = MIN( l_grid(k), l_stable )
+                   ll = l_grid(k)
+                ENDIF
+#else
                 CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
+#endif
+
 !
 !--             Compute diffusion coefficients for momentum and heat
@@ -4429,6 +4620,8 @@
     ENDIF
 
+    !$ACC KERNELS PRESENT(sums_l_l)
     sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
                                                 ! data output
+    !$ACC END KERNELS
 !$OMP END PARALLEL
 
