Index: /palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- /palm/trunk/SOURCE/flow_statistics.f90	(revision 3657)
+++ /palm/trunk/SOURCE/flow_statistics.f90	(revision 3658)
@@ -427,4 +427,9 @@
 !--    Initialize (local) summation array
        sums_l = 0.0_wp
+#ifdef _OPENACC
+       !$ACC KERNELS PRESENT(sums_l)
+       sums_l = 0.0_wp
+       !$ACC END KERNELS
+#endif
 
 !
@@ -500,12 +505,17 @@
        !$ tn = omp_get_thread_num()
        !$OMP DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag) &
+       !$ACC PRESENT(wall_flags_0, u, v, pt, rmask, sums_l)
        DO  i = nxl, nxr
           DO  j =  nys, nyn
              DO  k = nzb, nzt+1
                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
+                !$ACC ATOMIC
                 sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)  &
                                                               * flag
+                !$ACC ATOMIC
                 sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)  &
                                                               * flag
+                !$ACC ATOMIC
                 sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)  &
                                                               * flag
@@ -513,4 +523,5 @@
           ENDDO
        ENDDO
+       !$ACC UPDATE HOST(sums_l(:,1,tn), sums_l(:,2,tn), sums_l(:,4,tn))
 
 !
@@ -670,4 +681,5 @@
        hom(:,1,2,sr) = sums(:,2)             ! v
        hom(:,1,4,sr) = sums(:,4)             ! pt
+       !$ACC UPDATE DEVICE(hom(:,1,1,sr), hom(:,1,2,sr), hom(:,1,4,sr))
 
 
@@ -714,4 +726,11 @@
        !$ tn = omp_get_thread_num()
        !$OMP DO
+       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, m) &
+       !$ACC PRIVATE(sums_l_etot, flag) &
+       !$ACC PRESENT(wall_flags_0, rmask, momentumflux_output_conversion) &
+       !$ACC PRESENT(hom(:,1,4,sr)) &
+       !$ACC PRESENT(e, u, v, w, km, kh, p, pt) &
+       !$ACC PRESENT(surf_def_h(0), surf_lsm_h, surf_usm_h) &
+       !$ACC PRESENT(sums_l)
        DO  i = nxl, nxr
           DO  j =  nys, nyn
@@ -721,20 +740,26 @@
 !
 !--             Prognostic and diagnostic variables
+                !$ACC ATOMIC
                 sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)  &
                                                               * flag
+                !$ACC ATOMIC
                 sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)  &
                                                               * flag
+                !$ACC ATOMIC
                 sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)  &
                                                               * flag
+                !$ACC ATOMIC
                 sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)  &
                                                               * flag
+                !$ACC ATOMIC
                 sums_l(k,40,tn) = sums_l(k,40,tn) + ( p(k,j,i)                 &
                                          / momentumflux_output_conversion(k) ) &
                                                               * flag
 
+                !$ACC ATOMIC
                 sums_l(k,33,tn) = sums_l(k,33,tn) + &
                                   ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)&
                                                                  * flag
-
+#ifndef _OPENACC
                 IF ( humidity )  THEN
                    sums_l(k,70,tn) = sums_l(k,70,tn) + &
@@ -747,7 +772,9 @@
                                                                   * flag
                 ENDIF
+#endif
 !
 !--             Higher moments
 !--             (Computation of the skewness of w further below)
+                !$ACC ATOMIC
                 sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr) &
                                                                 * flag
@@ -763,4 +790,5 @@
 !--          quantities is seperated from the previous loop in order to
 !--          allow vectorization of that loop.
+             !$ACC ATOMIC
              sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
 !
@@ -769,12 +797,17 @@
                   surf_def_h(0)%start_index(j,i) )  THEN
                 m = surf_def_h(0)%start_index(j,i)
+                !$ACC ATOMIC
                 sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +            &
                                         surf_def_h(0)%us(m)   * rmask(j,i,sr)
+                !$ACC ATOMIC
                 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +          &
                                         surf_def_h(0)%usws(m) * rmask(j,i,sr)
+                !$ACC ATOMIC
                 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +          &
                                         surf_def_h(0)%vsws(m) * rmask(j,i,sr)
+                !$ACC ATOMIC
                 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +          &
                                         surf_def_h(0)%ts(m)   * rmask(j,i,sr)
+#ifndef _OPENACC
                 IF ( humidity )  THEN
                    sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +     &
@@ -785,6 +818,8 @@
                                             surf_def_h(0)%ss(m)   * rmask(j,i,sr)
                 ENDIF
+#endif
 !
 !--             Summation of surface temperature.
+                !$ACC ATOMIC
                 sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn)   +      &
                                             surf_def_h(0)%pt_surface(m) *      &
@@ -793,12 +828,17 @@
              IF ( surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i) )  THEN
                 m = surf_lsm_h%start_index(j,i)
+                !$ACC ATOMIC
                 sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +            &
                                         surf_lsm_h%us(m)   * rmask(j,i,sr)
+                !$ACC ATOMIC
                 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +          &
                                         surf_lsm_h%usws(m) * rmask(j,i,sr)
+                !$ACC ATOMIC
                 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +          &
                                         surf_lsm_h%vsws(m) * rmask(j,i,sr)
+                !$ACC ATOMIC
                 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +          &
                                         surf_lsm_h%ts(m)   * rmask(j,i,sr)
+#ifndef _OPENACC
                 IF ( humidity )  THEN
                    sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +     &
@@ -809,6 +849,8 @@
                                             surf_lsm_h%ss(m)   * rmask(j,i,sr)
                 ENDIF
+#endif
 !
 !--             Summation of surface temperature.
+                !$ACC ATOMIC
                 sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn)   +      &
                                             surf_lsm_h%pt_surface(m)    *      &
@@ -817,12 +859,17 @@
              IF ( surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i) )  THEN
                 m = surf_usm_h%start_index(j,i)
+                !$ACC ATOMIC
                 sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +            &
                                         surf_usm_h%us(m)   * rmask(j,i,sr)
+                !$ACC ATOMIC
                 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +          &
                                         surf_usm_h%usws(m) * rmask(j,i,sr)
+                !$ACC ATOMIC
                 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +          &
                                         surf_usm_h%vsws(m) * rmask(j,i,sr)
+                !$ACC ATOMIC
                 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +          &
                                         surf_usm_h%ts(m)   * rmask(j,i,sr)
+#ifndef _OPENACC
                 IF ( humidity )  THEN
                    sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +     &
@@ -833,6 +880,8 @@
                                             surf_usm_h%ss(m)   * rmask(j,i,sr)
                 ENDIF
+#endif
 !
 !--             Summation of surface temperature.
+                !$ACC ATOMIC
                 sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn)   +      &
                                             surf_usm_h%pt_surface(m)    *      &
@@ -841,4 +890,9 @@
           ENDDO
        ENDDO
+       !$ACC UPDATE &
+       !$ACC HOST(sums_l(:,3,tn), sums_l(:,8,tn), sums_l(:,9,tn)) &
+       !$ACC HOST(sums_l(:,10,tn), sums_l(:,40,tn), sums_l(:,33,tn)) &
+       !$ACC HOST(sums_l(:,38,tn)) &
+       !$ACC HOST(sums_l(nzb:nzb+4,pr_palm,tn), sums_l(nzb+14:nzb+14,pr_palm,tn))
 
 !
@@ -882,4 +936,7 @@
 !--    term. 
        !$OMP DO
+       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag, w2, ust2, vst2) &
+       !$ACC PRESENT(wall_flags_0, u, v, w, rmask, hom(:,1,1:2,sr)) &
+       !$ACC PRESENT(sums_l)
        DO  i = nxl, nxr
           DO  j =  nys, nyn
@@ -892,4 +949,5 @@
                 w2   = w(k,j,i)**2
 
+                !$ACC ATOMIC
                 sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)            &
                                  + 0.5_wp * ( ust2 + vst2 + w2 )               &
@@ -899,4 +957,5 @@
           ENDDO
        ENDDO
+       !$ACC UPDATE HOST(sums_l(nzb+5:nzb+5,pr_palm,tn))
 
 !
@@ -904,4 +963,11 @@
 
        !$OMP DO
+       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
+       !$ACC PRIVATE(ki, flag, ust, vst, pts) &
+       !$ACC PRESENT(kh, km, u, v, w, pt) &
+       !$ACC PRESENT(wall_flags_0, rmask, ddzu, rho_air_zw, hom(:,1,1:4,sr)) &
+       !$ACC PRESENT(heatflux_output_conversion, momentumflux_output_conversion) &
+       !$ACC PRESENT(surf_def_h(0:2), surf_lsm_h, surf_usm_h) &
+       !$ACC PRESENT(sums_l)
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -923,4 +989,5 @@
 !
 !--             Momentum flux w"u"
+                !$ACC ATOMIC
                 sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * (                &
                                km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
@@ -934,4 +1001,5 @@
 !
 !--             Momentum flux w"v"
+                !$ACC ATOMIC
                 sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * (                &
                                km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
@@ -945,4 +1013,5 @@
 !
 !--             Heat flux w"pt"
+                !$ACC ATOMIC
                 sums_l(k,16,tn) = sums_l(k,16,tn)                              &
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
@@ -953,7 +1022,7 @@
                                                * flag
 
-
 !
 !--             Salinity flux w"sa"
+#ifndef _OPENACC
                 IF ( ocean_mode )  THEN
                    sums_l(k,65,tn) = sums_l(k,65,tn)                           &
@@ -1000,4 +1069,5 @@
                                                               * flag
                 ENDIF
+#endif
 
              ENDDO
@@ -1007,5 +1077,8 @@
              IF ( use_surface_fluxes )  THEN
                 DO  l = 0, 1
-                   ki = MERGE( -1, 0, l == 0 )
+                   ! The original code using MERGE doesn't work with the PGI
+                   ! compiler when running on the GPU.
+                   ! ki = MERGE( -1, 0, l == 0 )
+                   ki = -1 + l
                    IF ( surf_def_h(l)%ns >= 1 )  THEN
                       DO  m = surf_def_h(l)%start_index(j,i),                  &
@@ -1013,17 +1086,23 @@
                          k = surf_def_h(l)%k(m)
 
+                         !$ACC ATOMIC
                          sums_l(k+ki,12,tn) = sums_l(k+ki,12,tn) + &
                                     momentumflux_output_conversion(k+ki) * &
                                     surf_def_h(l)%usws(m) * rmask(j,i,sr)     ! w"u"
+                         !$ACC ATOMIC
                          sums_l(k+ki,14,tn) = sums_l(k+ki,14,tn) + &
                                     momentumflux_output_conversion(k+ki) * &
                                     surf_def_h(l)%vsws(m) * rmask(j,i,sr)     ! w"v"
+                         !$ACC ATOMIC
                          sums_l(k+ki,16,tn) = sums_l(k+ki,16,tn) + &
                                     heatflux_output_conversion(k+ki) * &
                                     surf_def_h(l)%shf(m)  * rmask(j,i,sr)     ! w"pt"
+#if 0
                          sums_l(k+ki,58,tn) = sums_l(k+ki,58,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! u"pt"
                          sums_l(k+ki,61,tn) = sums_l(k+ki,61,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! v"pt"
+#endif
+#ifndef _OPENACC
                          IF ( ocean_mode )  THEN
                             sums_l(k+ki,65,tn) = sums_l(k+ki,65,tn) + &
@@ -1058,4 +1137,5 @@
                                         surf_def_h(l)%ssws(m) * rmask(j,i,sr) ! w"s"
                          ENDIF
+#endif
 
                       ENDDO
@@ -1066,17 +1146,23 @@
                      surf_lsm_h%start_index(j,i) )  THEN
                    m = surf_lsm_h%start_index(j,i)
+                   !$ACC ATOMIC
                    sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
                                     momentumflux_output_conversion(nzb) * &
                                     surf_lsm_h%usws(m) * rmask(j,i,sr)     ! w"u"
+                   !$ACC ATOMIC
                    sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
                                     momentumflux_output_conversion(nzb) * &
                                     surf_lsm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
+                   !$ACC ATOMIC
                    sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
                                     heatflux_output_conversion(nzb) * &
                                     surf_lsm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
+#if 0
                    sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! u"pt"
                    sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! v"pt"
+#endif
+#ifndef _OPENACC
                    IF ( ocean_mode )  THEN
                       sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
@@ -1111,5 +1197,5 @@
                                         surf_lsm_h%ssws(m) * rmask(j,i,sr) ! w"s"
                    ENDIF
-
+#endif
 
                 ENDIF
@@ -1117,17 +1203,23 @@
                      surf_usm_h%start_index(j,i) )  THEN
                    m = surf_usm_h%start_index(j,i)
+                   !$ACC ATOMIC
                    sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
                                     momentumflux_output_conversion(nzb) * &
                                     surf_usm_h%usws(m) * rmask(j,i,sr)     ! w"u"
+                   !$ACC ATOMIC
                    sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
                                     momentumflux_output_conversion(nzb) * &
                                     surf_usm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
+                   !$ACC ATOMIC
                    sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
                                     heatflux_output_conversion(nzb) * &
                                     surf_usm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
+#if 0
                    sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! u"pt"
                    sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! v"pt"
+#endif
+#ifndef _OPENACC
                    IF ( ocean_mode )  THEN
                       sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
@@ -1162,5 +1254,5 @@
                                         surf_usm_h%ssws(m) * rmask(j,i,sr) ! w"s"
                    ENDIF
-
+#endif
 
                 ENDIF
@@ -1168,4 +1260,5 @@
              ENDIF
 
+#ifndef _OPENACC
              IF ( .NOT. neutral )  THEN
                 IF ( surf_def_h(0)%end_index(j,i) >=                           &
@@ -1284,22 +1377,40 @@
 #endif
              ENDIF
+#endif
 !
 !--          Subgridscale fluxes at the top surface
              IF ( use_top_fluxes )  THEN
                 m = surf_def_h(2)%start_index(j,i)
-                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
-                                    momentumflux_output_conversion(nzt:nzt+1) * &
+                !$ACC ATOMIC
+                sums_l(nzt,12,tn) = sums_l(nzt,12,tn) + &
+                                    momentumflux_output_conversion(nzt) * &
                                     surf_def_h(2)%usws(m) * rmask(j,i,sr)    ! w"u"
-                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
-                                    momentumflux_output_conversion(nzt:nzt+1) * &
+                !$ACC ATOMIC
+                sums_l(nzt+1,12,tn) = sums_l(nzt+1,12,tn) + &
+                                    momentumflux_output_conversion(nzt+1) * &
+                                    surf_def_h(2)%usws(m) * rmask(j,i,sr)    ! w"u"
+                !$ACC ATOMIC
+                sums_l(nzt,14,tn) = sums_l(nzt,14,tn) + &
+                                    momentumflux_output_conversion(nzt) * &
                                     surf_def_h(2)%vsws(m) * rmask(j,i,sr)    ! w"v"
-                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
-                                    heatflux_output_conversion(nzt:nzt+1) * &
+                !$ACC ATOMIC
+                sums_l(nzt+1,14,tn) = sums_l(nzt+1,14,tn) + &
+                                    momentumflux_output_conversion(nzt+1) * &
+                                    surf_def_h(2)%vsws(m) * rmask(j,i,sr)    ! w"v"
+                !$ACC ATOMIC
+                sums_l(nzt,16,tn) = sums_l(nzt,16,tn) + &
+                                    heatflux_output_conversion(nzt) * &
                                     surf_def_h(2)%shf(m)  * rmask(j,i,sr)   ! w"pt"
+                !$ACC ATOMIC
+                sums_l(nzt+1,16,tn) = sums_l(nzt+1,16,tn) + &
+                                    heatflux_output_conversion(nzt+1) * &
+                                    surf_def_h(2)%shf(m)  * rmask(j,i,sr)   ! w"pt"
+#if 0
                 sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! u"pt"
                 sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! v"pt" 
-
+#endif
+#ifndef _OPENACC
                 IF ( ocean_mode )  THEN
                    sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
@@ -1337,4 +1448,5 @@
                                         surf_def_h(2)%ssws(m) * rmask(j,i,sr) ! w"s"
                 ENDIF
+#endif
              ENDIF
 
@@ -1354,6 +1466,8 @@
 
 !--             Higher moments
+                !$ACC ATOMIC
                 sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
                                                     rmask(j,i,sr) * flag
+                !$ACC ATOMIC
                 sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
                                                     rmask(j,i,sr) * flag
@@ -1362,4 +1476,5 @@
 !--             Salinity flux and density (density does not belong to here,
 !--             but so far there is no other suitable place to calculate)
+#ifndef _OPENACC
                 IF ( ocean_mode )  THEN
                    IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
@@ -1447,8 +1562,10 @@
                                                        rmask(j,i,sr) * flag
                 ENDIF
+#endif
 
 !
 !--             Energy flux w*e*
 !--             has to be adjusted
+                !$ACC ATOMIC
                 sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
                                              ( ust**2 + vst**2 + w(k,j,i)**2 ) &
@@ -1460,4 +1577,8 @@
        ENDDO
        !$OMP END PARALLEL
+
+       !$ACC UPDATE &
+       !$ACC HOST(sums_l(:,12,tn), sums_l(:,14,tn), sums_l(:,16,tn)) &
+       !$ACC HOST(sums_l(:,35,tn), sums_l(:,36,tn), sums_l(:,37,tn))
 !
 !--    Treat land-surface quantities according to new wall model structure. 
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 3657)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 3658)
@@ -442,5 +442,6 @@
                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
+               v_p, v_stokes_zu, w, w_p, zu, heatflux_output_conversion,       &
+               momentumflux_output_conversion
 
     USE biometeorology_mod,                                                    &
@@ -585,5 +586,5 @@
     USE statistics,                                                            &
         ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l,                &
-               rmask, statistic_regions, weight_substep, sums_l_l,             &
+               rmask, statistic_regions, weight_substep, sums_l, 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,&
@@ -685,4 +686,6 @@
 !$ACC COPYIN(ddzu(1:nzt+1), dd2zu(1:nzt)) &
 !$ACC COPYIN(ddzw(1:nzt+1)) &
+!$ACC COPYIN(heatflux_output_conversion(nzb:nzt+1)) &
+!$ACC COPYIN(momentumflux_output_conversion(nzb:nzt+1)) &
 !$ACC COPYIN(rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt)) &
 !$ACC COPYIN(ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng)) &
@@ -715,7 +718,8 @@
 ! Copy data from statistics
 !$ACC DATA &
-!$ACC COPYIN(hom(nzb+1:nzb_max,1,1:3,0)) &
+!$ACC COPYIN(hom(0:nz+1,1:2,1:4,0)) &
 !$ACC COPYIN(rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions)) &
 !$ACC COPYIN(weight_substep(1:intermediate_timestep_count_max)) &
+!$ACC COPY(sums_l(nzb:nzt+1,1:pr_palm,0)) &
 !$ACC COPY(sums_l_l(nzb:nzt+1,0:statistic_regions,0)) &
 !$ACC COPY(sums_us2_ws_l(nzb:nzt+1,0)) &
@@ -734,4 +738,8 @@
 !$ACC COPY(sums_salsa_ws_l(nzb:nzt+1,0))
 
+#ifdef _OPENACC
+    CALL enter_surface_arrays
+#endif
+
 !
 !-- At beginning determine the first time step
@@ -849,8 +857,4 @@
        ENDIF   
        
-#ifdef _OPENACC
-       CALL enter_surface_arrays
-#endif
-        
 !
 !--    Start of intermediate step loop
@@ -1369,7 +1373,5 @@
 !
 !--    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 UPDATE &
        !$ACC HOST(sums_l_l(nzb:nzt+1,0:statistic_regions,0)) &
        !$ACC HOST(sums_us2_ws_l(nzb:nzt+1,0)) &
@@ -1387,7 +1389,4 @@
        !$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
 
 !
@@ -1810,4 +1809,7 @@
     ENDDO   ! time loop
 
+#ifdef _OPENACC
+    CALL exit_surface_arrays
+#endif
 !$ACC END DATA
 !$ACC END DATA
Index: /palm/trunk/SOURCE/timestep.f90
===================================================================
--- /palm/trunk/SOURCE/timestep.f90	(revision 3657)
+++ /palm/trunk/SOURCE/timestep.f90	(revision 3658)
@@ -211,4 +211,12 @@
 
     CALL cpu_log( log_point(12), 'calculate_timestep', 'start' )
+
+    !$ACC UPDATE &
+    !$ACC HOST(u(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+    !$ACC HOST(v(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+    !$ACC HOST(w(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+    !$ACC HOST(kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
+    !$ACC HOST(km(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
+
 !
 !--    Save former time step as reference
