Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 1220)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 1221)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! wall_flags_00 introduced, which holds bits 32-...
 !
 ! Former revisions:
@@ -1648,5 +1648,5 @@
 
           DO  k = nzb+1, nzb_max
-             ibit32 = IBITS(wall_flags_0(k,j,i),32,1)
+             ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
              ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
              ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
@@ -1808,5 +1808,5 @@
                                         )
 
-          ibit32 = IBITS(wall_flags_0(k,j,i),32,1)
+          ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
           ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
           ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
@@ -1845,7 +1845,7 @@
 !--       k index has to be modified near bottom and top, else array
 !--       subscripts will be exceeded.
-          ibit35 = IBITS(wall_flags_0(k,j,i),35,1)
-          ibit34 = IBITS(wall_flags_0(k,j,i),34,1)
-          ibit33 = IBITS(wall_flags_0(k,j,i),33,1)
+          ibit35 = IBITS(wall_flags_00(k,j,i),3,1)
+          ibit34 = IBITS(wall_flags_00(k,j,i),2,1)
+          ibit33 = IBITS(wall_flags_00(k,j,i),1,1)
 
           k_ppp = k + 3 * ibit35
@@ -1942,7 +1942,7 @@
 !--       k index has to be modified near bottom and top, else array
 !--       subscripts will be exceeded.
-          ibit35 = IBITS(wall_flags_0(k,j,i),35,1)
-          ibit34 = IBITS(wall_flags_0(k,j,i),34,1)
-          ibit33 = IBITS(wall_flags_0(k,j,i),33,1)
+          ibit35 = IBITS(wall_flags_00(k,j,i),3,1)
+          ibit34 = IBITS(wall_flags_00(k,j,i),2,1)
+          ibit33 = IBITS(wall_flags_00(k,j,i),1,1)
 
           k_ppp = k + 3 * ibit35
@@ -2454,5 +2454,5 @@
 !
 !--    Computation of fluxes and tendency terms
-       !$acc kernels present( ddzw, sk, tend, u, v, w, wall_flags_0 )
+       !$acc kernels present( ddzw, sk, tend, u, v, w, wall_flags_0, wall_flags_00 )
        !$acc loop
        DO  i = i_left, i_right
@@ -3155,5 +3155,5 @@
 !
 !--    Computation of fluxes and tendency terms
-       !$acc  kernels present( ddzw, tend, u, v, w, wall_flags_0 )
+       !$acc  kernels present( ddzw, tend, u, v, w, wall_flags_0, wall_flags_00 )
        !$acc  loop
        DO i = i_left, i_right
@@ -3872,5 +3872,5 @@
 !
 !--    Computation of fluxes and tendency terms
-       !$acc kernels present( ddzw, tend, u, v, w, wall_flags_0 )
+       !$acc kernels present( ddzw, tend, u, v, w, wall_flags_0, wall_flags_00 )
        !$acc loop
        DO  i = i_left, i_right
@@ -4226,5 +4226,5 @@
           DO  k = nzb+1, nzb_max
 
-             ibit32 = IBITS(wall_flags_0(k,j,i),32,1)
+             ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
              ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
              ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
@@ -4327,5 +4327,5 @@
                                               )
 
-                ibit32 = IBITS(wall_flags_0(k,j,i),32,1)
+                ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
                 ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
                 ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
@@ -4364,7 +4364,7 @@
 !--             k index has to be modified near bottom and top, else array
 !--             subscripts will be exceeded.
-                ibit35 = IBITS(wall_flags_0(k,j,i),35,1)
-                ibit34 = IBITS(wall_flags_0(k,j,i),34,1)
-                ibit33 = IBITS(wall_flags_0(k,j,i),33,1)
+                ibit35 = IBITS(wall_flags_00(k,j,i),3,1)
+                ibit34 = IBITS(wall_flags_00(k,j,i),2,1)
+                ibit33 = IBITS(wall_flags_00(k,j,i),1,1)
 
                 k_ppp = k + 3 * ibit35
@@ -4463,7 +4463,7 @@
 !--             k index has to be modified near bottom and top, else array
 !--             subscripts will be exceeded.
-                ibit35 = IBITS(wall_flags_0(k,j,i),35,1)
-                ibit34 = IBITS(wall_flags_0(k,j,i),34,1)
-                ibit33 = IBITS(wall_flags_0(k,j,i),33,1)
+                ibit35 = IBITS(wall_flags_00(k,j,i),3,1)
+                ibit34 = IBITS(wall_flags_00(k,j,i),2,1)
+                ibit33 = IBITS(wall_flags_00(k,j,i),1,1)
 
                 k_ppp = k + 3 * ibit35
@@ -4565,7 +4565,8 @@
        gv = 2.0 * v_gtrans
 
+
 !
 !--    Computation of fluxes and tendency terms
-       !$acc kernels present( ddzu, tend, u, v, w, wall_flags_0 )
+       !$acc kernels present( ddzu, tend, u, v, w, wall_flags_0, wall_flags_00 )
        !$acc loop
        DO i = i_left, i_right
@@ -4574,8 +4575,7 @@
              DO  k = nzb+1, nzt
 
+                ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
+                ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
                 ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
-                ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
-                ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
-
 
                 u_comp_l                 = u(k+1,j,i) + u(k,j,i) - gu
@@ -4640,6 +4640,5 @@
                                  ( w(k,j,i+3) - w(k,j,i-2) )                 &
                                               )
-
-                ibit32 = IBITS(wall_flags_0(k,j,i),32,1)
+                ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
                 ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
                 ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
@@ -4707,8 +4706,7 @@
                                               )
 
-
-                ibit35 = IBITS(wall_flags_0(k-1,j,i),35,1)
-                ibit34 = IBITS(wall_flags_0(k-1,j,i),34,1)
-                ibit33 = IBITS(wall_flags_0(k-1,j,i),33,1)
+                ibit35 = IBITS(wall_flags_00(k-1,j,i),3,1)
+                ibit34 = IBITS(wall_flags_00(k-1,j,i),2,1)
+                ibit33 = IBITS(wall_flags_00(k-1,j,i),1,1)
 
                 k_pp  = k + 2 * ibit35
@@ -4750,7 +4748,7 @@
 !--             k index has to be modified near bottom and top, else array
 !--             subscripts will be exceeded.
-                ibit35 = IBITS(wall_flags_0(k,j,i),35,1)
-                ibit34 = IBITS(wall_flags_0(k,j,i),34,1)
-                ibit33 = IBITS(wall_flags_0(k,j,i),33,1)
+                ibit35 = IBITS(wall_flags_00(k,j,i),3,1)
+                ibit34 = IBITS(wall_flags_00(k,j,i),2,1)
+                ibit33 = IBITS(wall_flags_00(k,j,i),1,1)
 
                 k_ppp = k + 3 * ibit35
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 1220)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 1221)
@@ -1,2 +1,3 @@
+#if ! defined( __openacc )
  SUBROUTINE flow_statistics
 
@@ -20,5 +21,5 @@
 ! Current revisions:
 ! -----------------
-!
+! ported for openACC in separate branch
 !
 ! Former revisions:
@@ -64,5 +65,5 @@
 ! scheme. Furthermore the calculation will be the same for all advection
 ! schemes.
-!
+!, tend
 ! 696 2011-03-18 07:03:49Z raasch
 ! Bugfix: Summation of Wicker-Skamarock scheme fluxes and variances for all
@@ -182,4 +183,6 @@
     CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
 
+    !$acc update host( km, kh, e, pt, qs, qsws, rif, shf, ts, u, usws, v, vsws, w )
+
 !
 !-- To be on the safe side, check whether flow_statistics has already been
@@ -192,6 +195,4 @@
 
     ENDIF
-
-    !$acc update host( km, kh, e, pt, qs, qsws, rif, shf, ts, u, v, w )
 
 !
@@ -1321,2 +1322,1623 @@
 
  END SUBROUTINE flow_statistics
+
+
+#else
+
+
+!------------------------------------------------------------------------------!
+! flow statistics - accelerator version
+!------------------------------------------------------------------------------!
+ SUBROUTINE flow_statistics
+
+    USE arrays_3d
+    USE cloud_parameters
+    USE control_parameters
+    USE cpulog
+    USE grid_variables
+    USE indices
+    USE interfaces
+    USE pegrid
+    USE statistics
+
+    IMPLICIT NONE
+
+    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
+    LOGICAL ::  first
+    REAL    ::  dptdz_threshold, height, pts, sums_l_eper, sums_l_etot, ust, &
+                ust2, u2, vst, vst2, v2, w2, z_i(2)
+    REAL    ::  s1, s2, s3, s4, s5, s6, s7
+    REAL    ::  dptdz(nzb+1:nzt+1)
+    REAL    ::  sums_ll(nzb:nzt+1,2)
+
+    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
+
+!
+!-- To be on the safe side, check whether flow_statistics has already been
+!-- called once after the current time step
+    IF ( flow_statistics_called )  THEN
+
+       message_string = 'flow_statistics is called two times within one ' // &
+                        'timestep'
+       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
+
+    ENDIF
+
+    !$acc data copyin( hom ) create( sums, sums_l )
+
+!
+!-- Compute statistics for each (sub-)region
+    DO  sr = 0, statistic_regions
+
+!
+!--    Initialize (local) summation array
+       sums_l = 0.0
+
+!
+!--    Store sums that have been computed in other subroutines in summation
+!--    array
+       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
+!--    WARNING: next line still has to be adjusted for OpenMP
+       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
+       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
+       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
+
+!
+!--    Copy the turbulent quantities, evaluated in the advection routines to
+!--    the local array sums_l() for further computations
+       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
+
+!
+!--       According to the Neumann bc for the horizontal velocity components,
+!--       the corresponding fluxes has to satisfiy the same bc.
+          IF ( ocean )  THEN
+             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
+             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
+          ENDIF
+
+          DO  i = 0, threads_per_task-1
+!
+!--          Swap the turbulent quantities evaluated in advec_ws.
+             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
+             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
+             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
+             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
+             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
+             sums_l(:,34,i) = sums_l(:,34,i) + 0.5 *                        &
+                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +   &
+                                sums_ws2_ws_l(:,i) )    ! e*
+             DO  k = nzb, nzt
+                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * ( &
+                                                      sums_us2_ws_l(k,i) +  &
+                                                      sums_vs2_ws_l(k,i) +  &
+                                                      sums_ws2_ws_l(k,i) )
+             ENDDO
+          ENDDO
+
+       ENDIF
+
+       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
+
+          DO  i = 0, threads_per_task-1
+             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
+             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
+             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
+                                                   sums_wsqs_ws_l(:,i) !w*q*
+          ENDDO
+
+       ENDIF
+!
+!--    Horizontally averaged profiles of horizontal velocities and temperature.
+!--    They must have been computed before, because they are already required
+!--    for other horizontal averages.
+       tn = 0
+
+       !$OMP PARALLEL PRIVATE( i, j, k, tn )
+#if defined( __intel_openmp_bug )
+       tn = omp_get_thread_num()
+#else
+!$     tn = omp_get_thread_num()
+#endif
+
+       !$acc update device( sums_l )
+
+       !$OMP DO
+       !$acc parallel loop gang present( pt, rflags_invers, rmask, sums_l, u, v ) create( s1, s2, s3 )
+       DO  k = nzb, nzt+1
+          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+!
+!--             k+1 is used in rflags since rflags is set 0 at surface points
+                s1 = s1 + u(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                s2 = s2 + v(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                s3 = s3 + pt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+             ENDDO
+          ENDDO
+          sums_l(k,1,tn) = s1
+          sums_l(k,2,tn) = s2
+          sums_l(k,4,tn) = s3
+       ENDDO
+       !$acc end parallel
+
+!
+!--    Horizontally averaged profile of salinity
+       IF ( ocean )  THEN
+          !$OMP DO
+          !$acc parallel loop gang present( rflags_invers, rmask, sums_l, sa ) create( s1 )
+          DO  k = nzb, nzt+1
+             !$acc loop vector collapse( 2 ) reduction( +: s1 )
+             DO  i = nxl, nxr
+                DO  j =  nys, nyn
+                   s1 = s1 + sa(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,23,tn) = s1
+          ENDDO
+          !$acc end parallel
+       ENDIF
+
+!
+!--    Horizontally averaged profiles of virtual potential temperature,
+!--    total water content, specific humidity and liquid water potential
+!--    temperature
+       IF ( humidity )  THEN
+
+          !$OMP DO
+          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l, vpt ) create( s1, s2 )
+          DO  k = nzb, nzt+1
+             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
+             DO  i = nxl, nxr
+                DO  j =  nys, nyn
+                   s1 = s1 + q(k,j,i)   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s2 = s2 + vpt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,41,tn) = s1
+             sums_l(k,44,tn) = s2
+          ENDDO
+          !$acc end parallel
+
+          IF ( cloud_physics )  THEN
+             !$OMP DO
+             !$acc parallel loop gang present( pt, q, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
+             DO  k = nzb, nzt+1
+                !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
+                DO  i = nxl, nxr
+                   DO  j =  nys, nyn
+                      s1 = s1 + ( q(k,j,i) - ql(k,j,i) ) * &
+                                rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                      s2 = s2 + ( pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) ) * &
+                                rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   ENDDO
+                ENDDO
+                sums_l(k,42,tn) = s1
+                sums_l(k,43,tn) = s2
+             ENDDO
+             !$acc end parallel
+          ENDIF
+       ENDIF
+
+!
+!--    Horizontally averaged profiles of passive scalar
+       IF ( passive_scalar )  THEN
+          !$OMP DO
+          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 )
+          DO  k = nzb, nzt+1
+             !$acc loop vector collapse( 2 ) reduction( +: s1 )
+             DO  i = nxl, nxr
+                DO  j =  nys, nyn
+                   s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,41,tn) = s1
+          ENDDO
+          !$acc end parallel
+       ENDIF
+       !$OMP END PARALLEL
+
+!
+!--    Summation of thread sums
+       IF ( threads_per_task > 1 )  THEN
+          DO  i = 1, threads_per_task-1
+             !$acc parallel present( sums_l )
+             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
+             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
+             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
+             !$acc end parallel
+             IF ( ocean )  THEN
+                !$acc parallel present( sums_l )
+                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
+                !$acc end parallel
+             ENDIF
+             IF ( humidity )  THEN
+                !$acc parallel present( sums_l )
+                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
+                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
+                !$acc end parallel
+                IF ( cloud_physics )  THEN
+                   !$acc parallel present( sums_l )
+                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
+                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
+                   !$acc end parallel
+                ENDIF
+             ENDIF
+             IF ( passive_scalar )  THEN
+                !$acc parallel present( sums_l )
+                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
+                !$acc end parallel
+             ENDIF
+          ENDDO
+       ENDIF
+
+#if defined( __parallel )
+!
+!--    Compute total sum from local sums
+       !$acc update host( sums_l )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
+                           MPI_SUM, comm2d, ierr )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
+                           MPI_SUM, comm2d, ierr )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
+                           MPI_SUM, comm2d, ierr )
+       IF ( ocean )  THEN
+          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
+                              MPI_REAL, MPI_SUM, comm2d, ierr )
+       ENDIF
+       IF ( humidity ) THEN
+          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
+                              MPI_REAL, MPI_SUM, comm2d, ierr )
+          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
+                              MPI_REAL, MPI_SUM, comm2d, ierr )
+          IF ( cloud_physics ) THEN
+             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
+                                 MPI_REAL, MPI_SUM, comm2d, ierr )
+             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
+                                 MPI_REAL, MPI_SUM, comm2d, ierr )
+          ENDIF
+       ENDIF
+
+       IF ( passive_scalar )  THEN
+          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
+                              MPI_REAL, MPI_SUM, comm2d, ierr )
+       ENDIF
+       !$acc update device( sums )
+#else
+       !$acc parallel present( sums, sums_l )
+       sums(:,1) = sums_l(:,1,0)
+       sums(:,2) = sums_l(:,2,0)
+       sums(:,4) = sums_l(:,4,0)
+       !$acc end parallel
+       IF ( ocean )  THEN
+          !$acc parallel present( sums, sums_l )
+          sums(:,23) = sums_l(:,23,0)
+          !$acc end parallel
+       ENDIF
+       IF ( humidity )  THEN
+          !$acc parallel present( sums, sums_l )
+          sums(:,44) = sums_l(:,44,0)
+          sums(:,41) = sums_l(:,41,0)
+          !$acc end parallel
+          IF ( cloud_physics )  THEN
+             !$acc parallel present( sums, sums_l )
+             sums(:,42) = sums_l(:,42,0)
+             sums(:,43) = sums_l(:,43,0)
+             !$acc end parallel
+          ENDIF
+       ENDIF
+       IF ( passive_scalar )  THEN
+          !$acc parallel present( sums, sums_l )
+          sums(:,41) = sums_l(:,41,0)
+          !$acc end parallel
+       ENDIF
+#endif
+
+!
+!--    Final values are obtained by division by the total number of grid points
+!--    used for summation. After that store profiles.
+       !$acc parallel present( hom, ngp_2dh, ngp_2dh_s_inner, sums )
+       sums(:,1) = sums(:,1) / ngp_2dh(sr)
+       sums(:,2) = sums(:,2) / ngp_2dh(sr)
+       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
+       hom(:,1,1,sr) = sums(:,1)             ! u
+       hom(:,1,2,sr) = sums(:,2)             ! v
+       hom(:,1,4,sr) = sums(:,4)             ! pt
+       !$acc end parallel
+
+!
+!--    Salinity
+       IF ( ocean )  THEN
+          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
+          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
+          hom(:,1,23,sr) = sums(:,23)             ! sa
+          !$acc end parallel
+       ENDIF
+
+!
+!--    Humidity and cloud parameters
+       IF ( humidity ) THEN
+          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
+          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
+          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
+          hom(:,1,44,sr) = sums(:,44)                ! vpt
+          hom(:,1,41,sr) = sums(:,41)                ! qv (q)
+          !$acc end parallel
+          IF ( cloud_physics ) THEN
+             !$acc parallel present( hom, ngp_2dh_s_inner, sums )
+             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
+             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
+             hom(:,1,42,sr) = sums(:,42)             ! qv
+             hom(:,1,43,sr) = sums(:,43)             ! pt
+             !$acc end parallel
+          ENDIF
+       ENDIF
+
+!
+!--    Passive scalar
+       IF ( passive_scalar )  THEN
+          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
+          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
+          hom(:,1,41,sr) = sums(:,41)                ! s (q)
+          !$acc end parallel
+       ENDIF
+
+!
+!--    Horizontally averaged profiles of the remaining prognostic variables,
+!--    variances, the total and the perturbation energy (single values in last
+!--    column of sums_l) and some diagnostic quantities.
+!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
+!--    ----  speaking the following k-loop would have to be split up and
+!--          rearranged according to the staggered grid.
+!--          However, this implies no error since staggered velocity components
+!--          are zero at the walls and inside buildings.
+       tn = 0
+#if defined( __intel_openmp_bug )
+       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
+       !$OMP                    tn, ust, ust2, u2, vst, vst2, v2, w2 )
+       tn = omp_get_thread_num()
+#else
+       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
+!$     tn = omp_get_thread_num()
+#endif
+       !$OMP DO
+       !$acc parallel loop gang present( e, hom, kh, km, p, pt, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, s5, s6, s7 )
+       DO  k = nzb, nzt+1
+          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5, s6, s7 )
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+!
+!--             Prognostic and diagnostic variables
+                s1 = s1 + w(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                s2 = s2 + e(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                s3 = s3 + km(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                s4 = s4 + kh(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                s5 = s5 + p(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                s6 = s6 + ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr) * &
+                          rflags_invers(j,i,k+1)
+!
+!--             Higher moments
+!--             (Computation of the skewness of w further below)
+                s7 = s7 + w(k,j,i)**3 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+             ENDDO
+          ENDDO
+          sums_l(k,3,tn)  = s1
+          sums_l(k,8,tn)  = s2
+          sums_l(k,9,tn)  = s3
+          sums_l(k,10,tn) = s4
+          sums_l(k,40,tn) = s5
+          sums_l(k,33,tn) = s6
+          sums_l(k,38,tn) = s7
+       ENDDO
+       !$acc end parallel
+
+       IF ( humidity )  THEN
+          !$OMP DO
+          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l ) create( s1 )
+          DO  k = nzb, nzt+1
+             !$acc loop vector collapse( 2 ) reduction( +: s1 )
+             DO  i = nxl, nxr
+                DO  j =  nys, nyn
+                   s1 = s1 + ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr) * &
+                             rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,70,tn) = s1
+          ENDDO
+          !$acc end parallel
+       ENDIF
+
+!
+!--    Total and perturbation energy for the total domain (being
+!--    collected in the last column of sums_l).
+       !$OMP DO
+       !$acc parallel loop collapse(3) present( rflags_invers, rmask, u, v, w ) reduction(+:s1)
+       DO  i = nxl, nxr
+          DO  j =  nys, nyn
+             DO  k = nzb, nzt+1
+                s1 = s1 + 0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) * &
+                          rmask(j,i,sr) * rflags_invers(j,i,k+1)
+             ENDDO
+          ENDDO
+       ENDDO
+       !$acc end parallel
+       !$acc parallel present( sums_l )
+       sums_l(nzb+4,pr_palm,tn) = s1
+       !$acc end parallel
+
+       !$OMP DO
+       !$acc parallel present( rmask, sums_l, us, usws, vsws, ts ) create( s1, s2, s3, s4 )
+       !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
+       DO  i = nxl, nxr
+          DO  j =  nys, nyn
+!
+!--          2D-arrays (being collected in the last column of sums_l)
+             s1 = s1 + us(j,i)   * rmask(j,i,sr)
+             s2 = s2 + usws(j,i) * rmask(j,i,sr)
+             s3 = s3 + vsws(j,i) * rmask(j,i,sr)
+             s4 = s4 + ts(j,i)   * rmask(j,i,sr)
+          ENDDO
+       ENDDO
+       sums_l(nzb,pr_palm,tn)   = s1
+       sums_l(nzb+1,pr_palm,tn) = s2
+       sums_l(nzb+2,pr_palm,tn) = s3
+       sums_l(nzb+3,pr_palm,tn) = s4
+       !$acc end parallel
+
+       IF ( humidity )  THEN
+          !$acc parallel present( qs, rmask, sums_l ) create( s1 )
+          !$acc loop vector collapse( 2 ) reduction( +: s1 )
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                s1 = s1 + qs(j,i) * rmask(j,i,sr)
+             ENDDO
+          ENDDO
+          sums_l(nzb+12,pr_palm,tn) = s1
+          !$acc end parallel
+       ENDIF
+
+!
+!--    Computation of statistics when ws-scheme is not used. Else these
+!--    quantities are evaluated in the advection routines.
+       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0 )  THEN
+
+          !$OMP DO
+          !$acc parallel loop gang present( u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, ust2, vst2, w2 )
+          DO  k = nzb, nzt+1
+             !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
+             DO  i = nxl, nxr
+                DO  j =  nys, nyn
+                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
+                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
+                   w2   = w(k,j,i)**2
+
+                   s1 = s1 + ust2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s2 = s2 + vst2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s3 = s3 + w2   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+!
+!--                Perturbation energy
+                   s4 = s4 + 0.5 * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) * &
+                             rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,30,tn) = s1
+             sums_l(k,31,tn) = s2
+             sums_l(k,32,tn) = s3
+             sums_l(k,34,tn) = s4
+          ENDDO
+          !$acc end parallel
+!
+!--       Total perturbation TKE
+          !$OMP DO
+          !$acc parallel present( sums_l ) create( s1 )
+          !$acc loop reduction( +: s1 )
+          DO  k = nzb, nzt+1
+             s1 = s1 + sums_l(k,34,tn)
+          ENDDO
+          sums_l(nzb+5,pr_palm,tn) = s1
+          !$acc end parallel
+
+       ENDIF
+
+!
+!--    Horizontally averaged profiles of the vertical fluxes
+
+!
+!--    Subgridscale fluxes.
+!--    WARNING: If a Prandtl-layer is used (k=nzb for flat terrain), the fluxes
+!--    -------  should be calculated there in a different way. This is done
+!--             in the next loop further below, where results from this loop are
+!--             overwritten. However, THIS WORKS IN CASE OF FLAT TERRAIN ONLY!
+!--             The non-flat case still has to be handled.
+!--    NOTE: for simplicity, nzb_s_inner is used below, although
+!--    ----  strictly speaking the following k-loop would have to be
+!--          split up according to the staggered grid.
+!--          However, this implies no error since staggered velocity
+!--          components are zero at the walls and inside buildings.
+       !$OMP DO
+       !$acc parallel loop gang present( ddzu, kh, km, pt, u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
+       DO  k = nzb, nzt_diff
+          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+
+!
+!--             Momentum flux w"u"
+                s1 = s1 - 0.25 * (                   &
+                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
+                                                           ) * (               &
+                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
+                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
+                                                               )               &
+                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+!
+!--             Momentum flux w"v"
+                s2 = s2 - 0.25 * (                   &
+                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
+                                                           ) * (               &
+                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
+                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
+                                                               )               &
+                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+!
+!--             Heat flux w"pt"
+                s3 = s3 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
+                              * ( pt(k+1,j,i) - pt(k,j,i) )   &
+                              * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+             ENDDO
+          ENDDO
+          sums_l(k,12,tn) = s1
+          sums_l(k,14,tn) = s2
+          sums_l(k,16,tn) = s3
+       ENDDO
+       !$acc end parallel
+
+!
+!--    Salinity flux w"sa"
+       IF ( ocean )  THEN
+          !$acc parallel loop gang present( ddzu, kh, sa, rflags_invers, rmask, sums_l ) create( s1 )
+          DO  k = nzb, nzt_diff
+             !$acc loop vector collapse( 2 ) reduction( +: s1 )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
+                                 * ( sa(k+1,j,i) - sa(k,j,i) )   &
+                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,65,tn) = s1
+          ENDDO
+          !$acc end parallel
+       ENDIF
+
+!
+!--    Buoyancy flux, water flux (humidity flux) w"q"
+       IF ( humidity ) THEN
+
+          !$acc parallel loop gang present( ddzu, kh, q, vpt, rflags_invers, rmask, sums_l ) create( s1, s2 )
+          DO  k = nzb, nzt_diff
+             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
+                                 * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
+                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s2 = s2 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
+                                 * ( q(k+1,j,i) - q(k,j,i) )     &
+                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,45,tn) = s1
+             sums_l(k,48,tn) = s2
+          ENDDO
+          !$acc end parallel
+
+          IF ( cloud_physics ) THEN
+
+             !$acc parallel loop gang present( ddzu, kh, q, ql, rflags_invers, rmask, sums_l ) create( s1 )
+             DO  k = nzb, nzt_diff
+                !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )    &
+                                    * ( ( q(k+1,j,i) - ql(k+1,j,i) ) &
+                                      - ( q(k,j,i) - ql(k,j,i) ) )   &
+                                    * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   ENDDO
+                ENDDO
+                sums_l(k,51,tn) = s1
+             ENDDO
+             !$acc end parallel
+
+          ENDIF
+
+       ENDIF
+!
+!--    Passive scalar flux
+       IF ( passive_scalar )  THEN
+
+          !$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 )
+          DO  k = nzb, nzt_diff
+             !$acc loop vector collapse( 2 ) reduction( +: s1 )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
+                                 * ( q(k+1,j,i) - q(k,j,i) )     &
+                                 * ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,48,tn) = s1
+          ENDDO
+          !$acc end parallel
+
+       ENDIF
+
+       IF ( use_surface_fluxes )  THEN
+
+          !$OMP DO
+          !$acc parallel present( rmask, shf, sums_l, usws, vsws ) create( s1, s2, s3, s4, s5 )
+          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+!
+!--             Subgridscale fluxes in the Prandtl layer
+                s1 = s1 + usws(j,i) * rmask(j,i,sr)     ! w"u"
+                s2 = s2 + vsws(j,i) * rmask(j,i,sr)     ! w"v"
+                s3 = s3 + shf(j,i)  * rmask(j,i,sr)     ! w"pt"
+                s4 = s4 + 0.0 * rmask(j,i,sr)           ! u"pt"
+                s5 = s5 + 0.0 * rmask(j,i,sr)           ! v"pt"
+             ENDDO
+          ENDDO
+          sums_l(nzb,12,tn) = s1
+          sums_l(nzb,14,tn) = s2
+          sums_l(nzb,16,tn) = s3
+          sums_l(nzb,58,tn) = s4
+          sums_l(nzb,61,tn) = s5
+          !$acc end parallel
+
+          IF ( ocean )  THEN
+
+             !$OMP DO
+             !$acc parallel present( rmask, saswsb, sums_l ) create( s1 )
+             !$acc loop vector collapse( 2 ) reduction( +: s1 )
+             DO  i = nxl, nxr
+                DO  j =  nys, nyn
+                   s1 = s1 + saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
+                ENDDO
+             ENDDO
+             sums_l(nzb,65,tn) = s1
+             !$acc end parallel
+
+          ENDIF
+
+          IF ( humidity )  THEN
+
+             !$OMP DO
+             !$acc parallel present( pt, q, qsws, rmask, shf, sums_l ) create( s1, s2 )
+             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
+             DO  i = nxl, nxr
+                DO  j =  nys, nyn
+                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                   s2 = s2 + ( ( 1.0 + 0.61 * q(nzb,j,i) ) * shf(j,i) &
+                               + 0.61 * pt(nzb,j,i) * qsws(j,i) )
+                ENDDO
+             ENDDO
+             sums_l(nzb,48,tn) = s1
+             sums_l(nzb,45,tn) = s2
+             !$acc end parallel
+
+             IF ( cloud_droplets )  THEN
+
+                !$OMP DO
+                !$acc parallel present( pt, q, ql, qsws, rmask, shf, sums_l ) create( s1 )
+                !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                DO  i = nxl, nxr
+                   DO  j =  nys, nyn
+                      s1 = s1 + ( ( 1.0 + 0.61 * q(nzb,j,i) - ql(nzb,j,i) ) * &
+                                  shf(j,i) + 0.61 * pt(nzb,j,i) * qsws(j,i) )
+                   ENDDO
+                ENDDO
+                sums_l(nzb,45,tn) = s1
+                !$acc end parallel
+
+             ENDIF
+
+             IF ( cloud_physics )  THEN
+
+                !$OMP DO
+                !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
+                !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                DO  i = nxl, nxr
+                   DO  j =  nys, nyn
+!
+!--                   Formula does not work if ql(nzb) /= 0.0
+                      s1 = s1 + qsws(j,i) * rmask(j,i,sr)   ! w"q" (w"qv")
+                   ENDDO
+                ENDDO
+                sums_l(nzb,51,tn) = s1
+                !$acc end parallel
+
+             ENDIF
+
+          ENDIF
+
+          IF ( passive_scalar )  THEN
+
+             !$OMP DO
+             !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
+             !$acc loop vector collapse( 2 ) reduction( +: s1 )
+             DO  i = nxl, nxr
+                DO  j =  nys, nyn
+                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                ENDDO
+             ENDDO
+             sums_l(nzb,48,tn) = s1
+             !$acc end parallel
+
+          ENDIF
+
+       ENDIF
+
+!
+!--    Subgridscale fluxes at the top surface
+       IF ( use_top_fluxes )  THEN
+
+          !$OMP DO
+          !$acc parallel present( rmask, sums_l, tswst, uswst, vswst ) create( s1, s2, s3, s4, s5 )
+          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                s1 = s1 + uswst(j,i) * rmask(j,i,sr)    ! w"u"
+                s2 = s2 + vswst(j,i) * rmask(j,i,sr)    ! w"v"
+                s3 = s3 + tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
+                s4 = s4 + 0.0 * rmask(j,i,sr)           ! u"pt"
+                s5 = s5 + 0.0 * rmask(j,i,sr)           ! v"pt"
+             ENDDO
+          ENDDO
+          sums_l(nzt:nzt+1,12,tn) = s1
+          sums_l(nzt:nzt+1,14,tn) = s2
+          sums_l(nzt:nzt+1,16,tn) = s3
+          sums_l(nzt:nzt+1,58,tn) = s4
+          sums_l(nzt:nzt+1,61,tn) = s5
+          !$acc end parallel
+
+          IF ( ocean )  THEN
+
+             !$OMP DO
+             !$acc parallel present( rmask, saswst, sums_l ) create( s1 )
+             !$acc loop vector collapse( 2 ) reduction( +: s1 )
+             DO  i = nxl, nxr
+                DO  j =  nys, nyn
+                   s1 = s1 + saswst(j,i) * rmask(j,i,sr)  ! w"sa"
+                ENDDO
+             ENDDO
+             sums_l(nzt,65,tn) = s1
+             !$acc end parallel
+
+          ENDIF
+
+          IF ( humidity )  THEN
+
+             !$OMP DO
+             !$acc parallel present( pt, q, qswst, rmask, tswst, sums_l ) create( s1, s2 )
+             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
+             DO  i = nxl, nxr
+                DO  j =  nys, nyn
+                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
+                   s2 = s2 + ( ( 1.0 + 0.61 * q(nzt,j,i) ) * tswst(j,i) + &
+                               0.61 * pt(nzt,j,i) * qswst(j,i) )
+                ENDDO
+             ENDDO
+             sums_l(nzt,48,tn) = s1
+             sums_l(nzt,45,tn) = s2
+             !$acc end parallel
+
+             IF ( cloud_droplets )  THEN
+
+                !$OMP DO
+                !$acc parallel present( pt, q, ql, qswst, rmask, tswst, sums_l ) create( s1 )
+                !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                DO  i = nxl, nxr
+                   DO  j =  nys, nyn
+                      s1 = s1 + ( ( 1.0 + 0.61 * q(nzt,j,i) - ql(nzt,j,i) ) * &
+                                  tswst(j,i) + 0.61 * pt(nzt,j,i) * qswst(j,i) )
+                   ENDDO
+                ENDDO
+                sums_l(nzt,45,tn) = s1
+                !$acc end parallel
+
+             ENDIF
+
+             IF ( cloud_physics )  THEN
+
+                !$OMP DO
+                !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
+                !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                DO  i = nxl, nxr
+                   DO  j =  nys, nyn
+!
+!--                   Formula does not work if ql(nzb) /= 0.0
+                      s1 = s1 + qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                   ENDDO
+                ENDDO
+                sums_l(nzt,51,tn) = s1
+                !$acc end parallel
+
+             ENDIF
+
+          ENDIF
+
+          IF ( passive_scalar )  THEN
+
+             !$OMP DO
+             !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
+             !$acc loop vector collapse( 2 ) reduction( +: s1 )
+             DO  i = nxl, nxr
+                DO  j =  nys, nyn
+                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
+                ENDDO
+             ENDDO
+             sums_l(nzt,48,tn) = s1
+             !$acc end parallel
+
+          ENDIF
+
+       ENDIF
+
+!
+!--    Resolved fluxes (can be computed for all horizontal points)
+!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
+!--    ----  speaking the following k-loop would have to be split up and
+!--          rearranged according to the staggered grid.
+       !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2, s3 )
+       DO  k = nzb, nzt_diff
+          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
+                              u(k+1,j,i) - hom(k+1,1,1,sr) )
+                vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
+                              v(k+1,j,i) - hom(k+1,1,2,sr) )
+                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
+                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
+!
+!--             Higher moments
+                s1 = s1 + pts * w(k,j,i)**2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                s2 = s2 + pts**2 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+!
+!--             Energy flux w*e* (has to be adjusted?)
+                s3 = s3 + w(k,j,i) * 0.5 * ( ust**2 + vst**2 + w(k,j,i)**2 ) &
+                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+             ENDDO
+          ENDDO
+          sums_l(k,35,tn) = s1
+          sums_l(k,36,tn) = s2
+          sums_l(k,37,tn) = s3
+       ENDDO
+       !$acc end parallel
+
+!
+!--    Salinity flux and density (density does not belong to here,
+!--    but so far there is no other suitable place to calculate)
+       IF ( ocean )  THEN
+
+          IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
+
+             !$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 )
+             DO  k = nzb, nzt_diff
+                !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      s1 = s1 + 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) +   &
+                                        sa(k+1,j,i) - hom(k+1,1,23,sr) ) &
+                                    * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   ENDDO
+                ENDDO
+                sums_l(k,66,tn) = s1
+             ENDDO
+             !$acc end parallel
+
+          ENDIF
+
+          !$acc parallel loop gang present( rflags_invers, rho, prho, rmask, sums_l ) create( s1, s2 )
+          DO  k = nzb, nzt_diff
+             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   s1 = s1 + rho(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s2 = s2 + prho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,64,tn) = s1
+             sums_l(k,71,tn) = s2
+          ENDDO
+          !$acc end parallel
+
+       ENDIF
+
+!
+!--    Buoyancy flux, water flux, humidity flux, liquid water
+!--    content, rain drop concentration and rain water content
+       IF ( humidity )  THEN
+
+          IF ( cloud_physics  .OR.  cloud_droplets )  THEN
+
+             !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
+             DO  k = nzb, nzt_diff
+                !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      s1 = s1 + 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
+                                        vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
+                                      w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   ENDDO
+                ENDDO
+                sums_l(k,46,tn) = s1
+             ENDDO
+             !$acc end parallel
+
+             IF ( .NOT. cloud_droplets )  THEN
+
+                !$acc parallel loop gang present( hom, q, ql, rflags_invers, rmask, sums_l, w ) create( s1 )
+                DO  k = nzb, nzt_diff
+                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                   DO  i = nxl, nxr
+                      DO  j = nys, nyn
+                         s1 = s1 + 0.5 * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
+                                           ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
+                                       * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                      ENDDO
+                   ENDDO
+                   sums_l(k,52,tn) = s1
+                ENDDO
+                !$acc end parallel
+
+                IF ( icloud_scheme == 0  )  THEN
+
+                   !$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
+                   DO  k = nzb, nzt_diff
+                      !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
+                      DO  i = nxl, nxr
+                         DO  j = nys, nyn
+                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                            s2 = s2 + qc(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                         ENDDO
+                      ENDDO
+                      sums_l(k,54,tn) = s1
+                      sums_l(k,75,tn) = s2
+                   ENDDO
+                   !$acc end parallel
+
+                   IF ( precipitation )  THEN
+
+                      !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
+                      DO  k = nzb, nzt_diff
+                         !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
+                         DO  i = nxl, nxr
+                            DO  j = nys, nyn
+                               s1 = s1 + nr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                               s2 = s2 + qr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                               s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                            ENDDO
+                         ENDDO
+                         sums_l(k,73,tn) = s1
+                         sums_l(k,74,tn) = s2
+                         sums_l(k,76,tn) = s3
+                      ENDDO
+                      !$acc end parallel
+
+                   ENDIF
+
+                ELSE
+
+                   !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
+                   DO  k = nzb, nzt_diff
+                      !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                      DO  i = nxl, nxr
+                         DO  j = nys, nyn
+                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                         ENDDO
+                      ENDDO
+                      sums_l(k,54,tn) = s1
+                   ENDDO
+                   !$acc end parallel
+
+                ENDIF
+
+             ELSE
+
+                !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
+                DO  k = nzb, nzt_diff
+                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                   DO  i = nxl, nxr
+                      DO  j = nys, nyn
+                         s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                      ENDDO
+                   ENDDO
+                   sums_l(k,54,tn) = s1
+                ENDDO
+                !$acc end parallel
+
+             ENDIF
+
+          ELSE
+
+             IF( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
+
+                !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
+                DO  k = nzb, nzt_diff
+                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                   DO  i = nxl, nxr
+                      DO  j = nys, nyn
+                         s1 = s1 + 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
+                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
+                                       * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                      ENDDO
+                   ENDDO
+                   sums_l(k,46,tn) = s1
+                ENDDO
+                !$acc end parallel
+
+             ELSEIF ( ws_scheme_sca  .AND.  sr == 0 )  THEN
+
+                !$acc parallel loop present( hom, sums_l )
+                DO  k = nzb, nzt_diff
+                   sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
+                                             0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
+                ENDDO
+                !$acc end parallel
+
+             ENDIF
+
+          ENDIF
+
+       ENDIF
+!
+!--    Passive scalar flux
+       IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
+
+          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
+          DO  k = nzb, nzt_diff
+             !$acc loop vector collapse( 2 ) reduction( +: s1 )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   s1 = s1 + 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) +   &
+                                     q(k+1,j,i) - hom(k+1,1,41,sr) ) &
+                                 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,49,tn) = s1
+          ENDDO
+          !$acc end parallel
+
+       ENDIF
+
+!
+!--    For speed optimization fluxes which have been computed in part directly
+!--    inside the WS advection routines are treated seperatly
+!--    Momentum fluxes first:
+       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0  )  THEN
+
+          !$OMP DO
+          !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2 )
+          DO  k = nzb, nzt_diff
+             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
+                               u(k+1,j,i) - hom(k+1,1,1,sr) )
+                   vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
+                               v(k+1,j,i) - hom(k+1,1,2,sr) )
+!
+!--                Momentum flux w*u*
+                   s1 = s1 + 0.5 * ( w(k,j,i-1) + w(k,j,i) ) &
+                                 * ust * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+!
+!--                Momentum flux w*v*
+                   s2 = s2 + 0.5 * ( w(k,j-1,i) + w(k,j,i) ) &
+                                 * vst * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,13,tn) = s1
+             sums_l(k,15,tn) = s1
+          ENDDO
+          !$acc end parallel
+
+       ENDIF
+
+       IF ( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
+
+          !$OMP DO
+          !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, w ) create( s1 )
+          DO  k = nzb, nzt_diff
+             !$acc loop vector collapse( 2 ) reduction( +: s1 )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+!
+!--                Vertical heat flux
+                   s1 = s1 + 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
+                                     pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
+                                 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                ENDDO
+             ENDDO
+             sums_l(k,17,tn) = s1
+          ENDDO
+          !$acc end parallel
+
+          IF ( humidity )  THEN
+
+             !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
+             DO  k = nzb, nzt_diff
+                !$acc loop vector collapse( 2 ) reduction( +: s1 )
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      s1 = s1 + 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) +   &
+                                        q(k+1,j,i) - hom(k+1,1,41,sr) ) &
+                                    * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   ENDDO
+                ENDDO
+                sums_l(k,49,tn) = s1
+             ENDDO
+             !$acc end parallel
+
+          ENDIF
+
+       ENDIF
+
+
+!
+!--    Density at top follows Neumann condition
+       IF ( ocean )  THEN
+          !$acc parallel present( sums_l )
+          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
+          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
+          !$acc end parallel
+       ENDIF
+
+!
+!--    Divergence of vertical flux of resolved scale energy and pressure
+!--    fluctuations as well as flux of pressure fluctuation itself (68).
+!--    First calculate the products, then the divergence.
+!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
+       IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
+
+          STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
+          sums_ll = 0.0  ! local array
+
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb_s_inner(j,i)+1, nzt
+
+                   sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
+                  ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
+                           - 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
+                           ) )**2                                          &
+                + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
+                           - 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
+                           ) )**2                                          &
+                   + w(k,j,i)**2                                  )
+
+                   sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
+                                               * ( p(k,j,i) + p(k+1,j,i) )
+
+                ENDDO
+             ENDDO
+          ENDDO
+          sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
+          sums_ll(nzt+1,1) = 0.0
+          sums_ll(0,2)     = 0.0
+          sums_ll(nzt+1,2) = 0.0
+
+          DO  k = nzb+1, nzt
+             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
+             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
+             sums_l(k,68,tn) = sums_ll(k,2)
+          ENDDO
+          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
+          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
+          sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
+
+       ENDIF
+
+!
+!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
+       IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
+
+          STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb_s_inner(j,i)+1, nzt
+
+                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
+                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
+                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
+                                                             ) * ddzw(k)
+
+                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
+                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
+                                                              )
+
+                ENDDO
+             ENDDO
+          ENDDO
+          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
+          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
+
+       ENDIF
+
+!
+!--    Horizontal heat fluxes (subgrid, resolved, total).
+!--    Do it only, if profiles shall be plotted.
+       IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
+
+          STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb_s_inner(j,i)+1, nzt
+!
+!--                Subgrid horizontal heat fluxes u"pt", v"pt"
+                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
+                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
+                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
+                                                 * ddx * rmask(j,i,sr)
+                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
+                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
+                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
+                                                 * ddy * rmask(j,i,sr)
+!
+!--                Resolved horizontal heat fluxes u*pt*, v*pt*
+                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
+                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
+                                       * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
+                                                 pt(k,j,i)   - hom(k,1,4,sr) )
+                   pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
+                                 pt(k,j,i)   - hom(k,1,4,sr) )
+                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
+                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
+                                       * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
+                                                 pt(k,j,i)   - hom(k,1,4,sr) )
+                ENDDO
+             ENDDO
+          ENDDO
+!
+!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
+          sums_l(nzb,58,tn) = 0.0
+          sums_l(nzb,59,tn) = 0.0
+          sums_l(nzb,60,tn) = 0.0
+          sums_l(nzb,61,tn) = 0.0
+          sums_l(nzb,62,tn) = 0.0
+          sums_l(nzb,63,tn) = 0.0
+
+       ENDIF
+
+!
+!--    Calculate the user-defined profiles
+       CALL user_statistics( 'profiles', sr, tn )
+       !$OMP END PARALLEL
+
+!
+!--    Summation of thread sums
+       IF ( threads_per_task > 1 )  THEN
+          STOP '+++ openACC porting for threads_per_task > 1 in flow_statistics is still missing'
+          DO  i = 1, threads_per_task-1
+             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
+             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
+             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
+                                      sums_l(:,45:pr_palm,i)
+             IF ( max_pr_user > 0 )  THEN
+                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
+                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
+                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
+             ENDIF
+          ENDDO
+       ENDIF
+
+       !$acc update host( hom, sums, sums_l )
+
+#if defined( __parallel )
+
+!
+!--    Compute total sum from local sums
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
+                           MPI_SUM, comm2d, ierr )
+#else
+       sums = sums_l(:,:,0)
+#endif
+
+!
+!--    Final values are obtained by division by the total number of grid points
+!--    used for summation. After that store profiles.
+!--    Profiles:
+       DO  k = nzb, nzt+1
+          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
+          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
+          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
+          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
+          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
+          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
+          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
+          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
+          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
+          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
+          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
+          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
+          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
+          sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
+       ENDDO
+
+!--    Upstream-parts
+       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
+!--    u* and so on
+!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
+!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
+!--    above the topography, they are being divided by ngp_2dh(sr)
+       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
+                                    ngp_2dh(sr)
+       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
+                                    ngp_2dh(sr)
+!--    eges, e*
+       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
+                                    ngp_3d(sr)
+!--    Old and new divergence
+       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
+                                    ngp_3d_inner(sr)
+
+!--    User-defined profiles
+       IF ( max_pr_user > 0 )  THEN
+          DO  k = nzb, nzt+1
+             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
+                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
+                                    ngp_2dh_s_inner(k,sr)
+          ENDDO
+       ENDIF
+
+!
+!--    Collect horizontal average in hom.
+!--    Compute deduced averages (e.g. total heat flux)
+       hom(:,1,3,sr)  = sums(:,3)      ! w
+       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
+       hom(:,1,9,sr)  = sums(:,9)      ! km
+       hom(:,1,10,sr) = sums(:,10)     ! kh
+       hom(:,1,11,sr) = sums(:,11)     ! l
+       hom(:,1,12,sr) = sums(:,12)     ! w"u"
+       hom(:,1,13,sr) = sums(:,13)     ! w*u*
+       hom(:,1,14,sr) = sums(:,14)     ! w"v"
+       hom(:,1,15,sr) = sums(:,15)     ! w*v*
+       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
+       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
+       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
+       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
+       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
+       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
+       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
+                                       ! profile 24 is initial profile (sa)
+                                       ! profiles 25-29 left empty for initial
+                                       ! profiles
+       hom(:,1,30,sr) = sums(:,30)     ! u*2
+       hom(:,1,31,sr) = sums(:,31)     ! v*2
+       hom(:,1,32,sr) = sums(:,32)     ! w*2
+       hom(:,1,33,sr) = sums(:,33)     ! pt*2
+       hom(:,1,34,sr) = sums(:,34)     ! e*
+       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
+       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
+       hom(:,1,37,sr) = sums(:,37)     ! w*e*
+       hom(:,1,38,sr) = sums(:,38)     ! w*3
+       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5   ! Sw
+       hom(:,1,40,sr) = sums(:,40)     ! p
+       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
+       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*
+       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
+       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
+       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
+       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
+       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
+       hom(:,1,52,sr) = sums(:,52)     ! w*qv*
+       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
+       hom(:,1,54,sr) = sums(:,54)     ! ql
+       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
+       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
+       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
+       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
+       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
+       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
+       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
+       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
+       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
+       hom(:,1,64,sr) = sums(:,64)     ! rho
+       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
+       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
+       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
+       hom(:,1,68,sr) = sums(:,68)     ! w*p*
+       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
+       hom(:,1,70,sr) = sums(:,70)     ! q*2
+       hom(:,1,71,sr) = sums(:,71)     ! prho
+       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
+       hom(:,1,73,sr) = sums(:,73)     ! nr
+       hom(:,1,74,sr) = sums(:,74)     ! qr
+       hom(:,1,75,sr) = sums(:,75)     ! qc
+       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
+                                       ! 77 is initial density profile
+
+       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
+                                       ! upstream-parts u_x, u_y, u_z, v_x,
+                                       ! v_y, usw. (in last but one profile)
+       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
+                                       ! u*, w'u', w'v', t* (in last profile)
+
+       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
+          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
+                               sums(:,pr_palm+1:pr_palm+max_pr_user)
+       ENDIF
+
+!
+!--    Determine the boundary layer height using two different schemes.
+!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
+!--    first relative minimum (maximum) of the total heat flux.
+!--    The corresponding height is assumed as the boundary layer height, if it
+!--    is less than 1.5 times the height where the heat flux becomes negative
+!--    (positive) for the first time.
+       z_i(1) = 0.0
+       first = .TRUE.
+
+       IF ( ocean )  THEN
+          DO  k = nzt, nzb+1, -1
+             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
+                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
+                first = .FALSE.
+                height = zw(k)
+             ENDIF
+             IF ( hom(k,1,18,sr) < 0.0  .AND. &
+                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
+                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
+                IF ( zw(k) < 1.5 * height )  THEN
+                   z_i(1) = zw(k)
+                ELSE
+                   z_i(1) = height
+                ENDIF
+                EXIT
+             ENDIF
+          ENDDO
+       ELSE
+          DO  k = nzb, nzt-1
+             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
+                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
+                first = .FALSE.
+                height = zw(k)
+             ENDIF
+             IF ( hom(k,1,18,sr) < 0.0  .AND. &
+                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
+                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
+                IF ( zw(k) < 1.5 * height )  THEN
+                   z_i(1) = zw(k)
+                ELSE
+                   z_i(1) = height
+                ENDIF
+                EXIT
+             ENDIF
+          ENDDO
+       ENDIF
+
+!
+!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
+!--    by Uhlenbrock(2006). The boundary layer height is the height with the
+!--    maximal local temperature gradient: starting from the second (the last
+!--    but one) vertical gridpoint, the local gradient must be at least
+!--    0.2K/100m and greater than the next four gradients.
+!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
+!--             ocean case!
+       z_i(2) = 0.0
+       DO  k = nzb+1, nzt+1
+          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
+       ENDDO
+       dptdz_threshold = 0.2 / 100.0
+
+       IF ( ocean )  THEN
+          DO  k = nzt+1, nzb+5, -1
+             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
+                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
+                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
+                z_i(2) = zw(k-1)
+                EXIT
+             ENDIF
+          ENDDO
+       ELSE
+          DO  k = nzb+1, nzt-3
+             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
+                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
+                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
+                z_i(2) = zw(k-1)
+                EXIT
+             ENDIF
+          ENDDO
+       ENDIF
+
+       hom(nzb+6,1,pr_palm,sr) = z_i(1)
+       hom(nzb+7,1,pr_palm,sr) = z_i(2)
+
+!
+!--    Computation of both the characteristic vertical velocity and
+!--    the characteristic convective boundary layer temperature.
+!--    The horizontal average at nzb+1 is input for the average temperature.
+       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
+           .AND.  z_i(1) /= 0.0 )  THEN
+          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
+                                       hom(nzb,1,18,sr) *      &
+                                       ABS( z_i(1) ) )**0.333333333
+!--       so far this only works if Prandtl layer is used
+          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
+       ELSE
+          hom(nzb+8,1,pr_palm,sr)  = 0.0
+          hom(nzb+11,1,pr_palm,sr) = 0.0
+       ENDIF
+
+!
+!--    Collect the time series quantities
+       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
+       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
+       ts_value(3,sr) = dt_3d
+       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
+       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
+       ts_value(6,sr) = u_max
+       ts_value(7,sr) = v_max
+       ts_value(8,sr) = w_max
+       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
+       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
+       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
+       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
+       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
+       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
+       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
+       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
+       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
+       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
+       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
+       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
+       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
+
+       IF ( ts_value(5,sr) /= 0.0 )  THEN
+          ts_value(22,sr) = ts_value(4,sr)**2 / &
+                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
+       ELSE
+          ts_value(22,sr) = 10000.0
+       ENDIF
+
+       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
+!
+!--    Calculate additional statistics provided by the user interface
+       CALL user_statistics( 'time_series', sr, 0 )
+
+    ENDDO    ! loop of the subregions
+
+    !$acc end data
+
+!
+!-- If required, sum up horizontal averages for subsequent time averaging
+    IF ( do_sum )  THEN
+       IF ( average_count_pr == 0 )  hom_sum = 0.0
+       hom_sum = hom_sum + hom(:,1,:,:)
+       average_count_pr = average_count_pr + 1
+       do_sum = .FALSE.
+    ENDIF
+
+!
+!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
+!-- This flag is reset after each time step in time_integration.
+    flow_statistics_called = .TRUE.
+
+    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
+
+
+ END SUBROUTINE flow_statistics
+#endif
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 1220)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 1221)
@@ -23,5 +23,6 @@
 ! Current revisions:
 ! ------------------
-!
+! +rflags_s_inner in copyin statement, use copyin for most arrays instead of
+! copy
 !
 ! Former revisions:
@@ -1499,5 +1500,8 @@
        CALL disturb_field( nzb_v_inner, tend, v )
        n_sor = nsor_ini
-       !$acc data copy( d, ddzu, ddzw, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, p, tri, tric, u, v, w, weight_pres, weight_substep, tend )
+       !$acc data copyin( d, ddzu, ddzw, nzb_s_inner, nzb_u_inner )            &
+       !$acc      copyin( nzb_v_inner, nzb_w_inner, p, rflags_s_inner, tend )  &
+       !$acc      copyin( weight_pres, weight_substep )                        &
+       !$acc      copy( tri, tric, u, v, w )
        CALL pres
        !$acc end data
Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 1220)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 1221)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! -----------------
-! 
+! wall_flags_00 introduced to hold bits 32-63,
+! additional 3D-flag arrays for replacing the 2D-index array nzb_s_inner in
+! loops optimized for openACC (pres + flow_statistics)
 !
 ! Former revisions:
@@ -387,4 +389,6 @@
               nzb_diff_v(nysg:nyng,nxlg:nxrg),                              &
               nzb_2d(nysg:nyng,nxlg:nxrg),                                  &
+              rflags_s_inner(nzb:nzt+2,nysg:nyng,nxlg:nxrg),                &
+              rflags_invers(nysg:nyng,nxlg:nxrg,nzb:nzt+2),                 &
               wall_e_x(nysg:nyng,nxlg:nxrg),                                &
               wall_e_y(nysg:nyng,nxlg:nxrg),                                &
@@ -403,4 +407,7 @@
     nzb_v_inner = nzb;  nzb_v_outer = nzb
     nzb_w_inner = nzb;  nzb_w_outer = nzb
+
+    rflags_s_inner = 1.0
+    rflags_invers  = 1.0
 
 !
@@ -854,4 +861,14 @@
           
        ENDIF
+!
+!--    Set flag arrays to be used for masking of grid points
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             DO  k = nzb, nzt+1
+                IF ( k <= nzb_s_inner(j,i) )  rflags_s_inner(k,j,i) = 0.0
+                IF ( k <= nzb_s_inner(j,i) )  rflags_invers(j,i,k)  = 0.0
+             ENDDO
+          ENDDO
+       ENDDO
 #endif
     ENDIF
@@ -1107,6 +1124,8 @@
 !
 !-- Allocate flags needed for masking walls.
-    ALLOCATE( wall_flags_0(nzb:nzt,nys:nyn,nxl:nxr) )
-    wall_flags_0 = 0
+    ALLOCATE( wall_flags_0(nzb:nzt,nys:nyn,nxl:nxr), &
+              wall_flags_00(nzb:nzt,nys:nyn,nxl:nxr) )
+    wall_flags_0  = 0
+    wall_flags_00 = 0
 
     IF ( scalar_advec == 'ws-scheme' )  THEN
@@ -1311,5 +1330,5 @@
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 )
                 ELSE
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 32 )
+                   wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 0 )
                 ENDIF
 !
@@ -1325,11 +1344,11 @@
 !--                because flux_t(nzb_w_inner(j,i)) is used for the tendency
 !--                at k == nzb_w_inner(j,i)+1.
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 33 )
+                   wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 1 )
                    flag_set = .TRUE.
                 ELSEIF ( k == nzb_w_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 34 )
+                   wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 2 )
                    flag_set = .TRUE.
                 ELSEIF ( k > nzb_w_inner(j,i) .AND. .NOT. flag_set )  THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 35 )
+                   wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 3 )
                 ENDIF
 
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 1220)
+++ palm/trunk/SOURCE/modules.f90	(revision 1221)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-!
+! wall_flags_0 changed to 32bit int, +wall_flags_00,
+! +rflags_s_inner, rflags_invers
 !
 ! Former revisions:
@@ -1088,5 +1089,5 @@
 ! Description:
 ! ------------
-! Definition of array bounds and number of gridpoints
+! Definition of array bounds, number of gridpoints, and wall flag arrays
 !------------------------------------------------------------------------------!
 
@@ -1112,6 +1113,5 @@
     INTEGER, DIMENSION(:,:,:), POINTER ::  flags
 
-    INTEGER( KIND = SELECTED_INT_KIND(11) ), DIMENSION(:,:,:), ALLOCATABLE ::  &
-                wall_flags_0  ! need to have 34 Bit
+    INTEGER, DIMENSION(:,:,:), ALLOCATABLE ::  wall_flags_0, wall_flags_00
 
     INTEGER, DIMENSION(:,:,:), ALLOCATABLE,  TARGET ::                         &
@@ -1120,4 +1120,5 @@
                 wall_flags_9, wall_flags_10
 
+    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  rflags_s_inner, rflags_invers
 
     SAVE
Index: palm/trunk/SOURCE/palm.f90
===================================================================
--- palm/trunk/SOURCE/palm.f90	(revision 1220)
+++ palm/trunk/SOURCE/palm.f90	(revision 1221)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! +wall_flags_00, rflags_invers, rflags_s_inner in copyin statement
 !
 ! Former revisions:
@@ -252,5 +252,6 @@
     !$acc       copyin( nzb_diff_v, nzb_s_inner, nzb_s_outer, nzb_u_inner )    &
     !$acc       copyin( nzb_u_outer, nzb_v_inner, nzb_v_outer, nzb_w_inner )   &
-    !$acc       copyin( nzb_w_outer, wall_heatflux, wall_e_x, wall_e_y, wall_u, wall_v, wall_w_x, wall_w_y, wall_flags_0 )  &
+    !$acc       copyin( nzb_w_outer, rflags_invers, rflags_s_inner, rmask, wall_heatflux, wall_e_x, wall_e_y, wall_u, wall_v, wall_w_x, wall_w_y, wall_flags_0, wall_flags_00 )  &
+    !$acc       copyin( ngp_2dh, ngp_2dh_s_inner )  &
     !$acc       copyin( weight_pres, weight_substep )
 !
Index: palm/trunk/SOURCE/pres.f90
===================================================================
--- palm/trunk/SOURCE/pres.f90	(revision 1220)
+++ palm/trunk/SOURCE/pres.f90	(revision 1221)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-!
+! openACC porting of reduction operations, loops for calculating d are
+! using the rflags_s_inner multiply flag instead of the nzb_s_inner loop index
 !
 ! Former revisions:
@@ -360,16 +361,13 @@
     !$OMP PARALLEL PRIVATE (i,j,k)
     !$OMP DO SCHEDULE( STATIC )
-    !$acc kernels present( d, ddzw, nzb_s_inner, u, v, w )
-    !$acc loop
+    !$acc kernels present( d, ddzw, rflags_s_inner, u, v, w )
+    !$acc loop collapse( 3 )
     DO  i = nxl, nxr
        DO  j = nys, nyn
-          !$acc loop vector(32)
           DO  k = 1, nzt
-             IF ( k > nzb_s_inner(j,i) )  THEN
-                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
-                           ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
-                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
-                           * d_weight_pres
-             ENDIF
+             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx +             &
+                        ( v(k,j+1,i) - v(k,j,i) ) * ddy +               &
+                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d  &
+                        * d_weight_pres * rflags_s_inner(k,j,i)
           ENDDO
        ENDDO
@@ -382,4 +380,5 @@
     !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
     !$OMP DO SCHEDULE( STATIC )
+    !$acc parallel loop collapse(3) present( d ) reduction(+:threadsum)
     DO  i = nxl, nxr
        DO  j = nys, nyn
@@ -389,4 +388,5 @@
        ENDDO
     ENDDO
+    !$acc end parallel
     localsum = localsum + threadsum * dt_3d * &
                           weight_pres(intermediate_timestep_count)
@@ -594,5 +594,4 @@
 !-- Correction of the provisional velocities with the current perturbation 
 !-- pressure just computed
-    !$acc update host( u, v, w )
     IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
        volume_flow_l(1) = 0.0
@@ -748,15 +747,13 @@
     ENDDO
 #else
-    !$acc kernels present( d, ddzw, nzb_s_inner, u, v, w )
-    !$acc loop
+    !$acc kernels present( d, ddzw, rflags_s_inner, u, v, w )
+    !$acc loop collapse( 3 )
     DO  i = nxl, nxr
        DO  j = nys, nyn
-          !$acc loop vector( 32 )
           DO  k = 1, nzt
-             IF ( k > nzb_s_inner(j,i) )  THEN
-                d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
-                           ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
-                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
-             ENDIF
+                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx +   &
+                             ( v(k,j+1,i) - v(k,j,i) ) * ddy +   &
+                             ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
+                           ) * rflags_s_inner(k,j,i)
           ENDDO
        ENDDO
@@ -765,11 +762,13 @@
 !
 !-- Compute possible PE-sum of divergences for flow_statistics
+    !$acc parallel loop collapse(3) present( d ) reduction(+:threadsum)
     DO  i = nxl, nxr
        DO  j = nys, nyn
-          DO  k = nzb_s_inner(j,i)+1, nzt
+          DO  k = nzb+1, nzt
              threadsum = threadsum + ABS( d(k,j,i) )
           ENDDO
        ENDDO
     ENDDO
+    !$acc end parallel
 #endif
 
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 1220)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 1221)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! host update of arrays before timestep is called
 !
 ! Former revisions:
@@ -211,7 +211,12 @@
 
        CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
+
 !
 !--    Determine size of next time step
-       IF ( simulated_time /= 0.0 )  CALL timestep
+       IF ( simulated_time /= 0.0 )  THEN
+          !$acc update host( kh, km, u, v, w )
+          CALL timestep
+       ENDIF
+
 !
 !--    Execute the user-defined actions
Index: palm/trunk/SOURCE/tridia_solver.f90
===================================================================
--- palm/trunk/SOURCE/tridia_solver.f90	(revision 1220)
+++ palm/trunk/SOURCE/tridia_solver.f90	(revision 1221)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-!
+! dummy argument tri in 1d-routines replaced by tri_for_1d because of name
+! conflict with arry tri in module arrays_3d
 !
 ! Former revisions:
@@ -382,5 +383,5 @@
 
 
-    SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri )
+    SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri_for_1d )
 
 !------------------------------------------------------------------------------!
@@ -408,5 +409,5 @@
 
        REAL, DIMENSION(0:nx,1:nz)     ::  ar
-       REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
+       REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d
 
 
@@ -424,6 +425,6 @@
        DO  k = 0, nz-1
           DO  i = 0,nx
-             tri(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
-             tri(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
+             tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
+             tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
           ENDDO
        ENDDO
@@ -433,5 +434,5 @@
        IF ( j <= nnyh )  THEN
 #if defined( __intel11 )
-          CALL maketri_1dd( j, tri )
+          CALL maketri_1dd( j, tri_for_1d )
 #else
           CALL maketri_1dd( j )
@@ -439,5 +440,5 @@
        ELSE
 #if defined( __intel11 )
-          CALL maketri_1dd( ny+1-j, tri )
+          CALL maketri_1dd( ny+1-j, tri_for_1d )
 #else
           CALL maketri_1dd( ny+1-j )
@@ -445,14 +446,14 @@
        ENDIF
 #if defined( __intel11 )
-       CALL split_1dd( tri )
+       CALL split_1dd( tri_for_1d )
 #else
        CALL split_1dd
 #endif
-       CALL substi_1dd( ar, tri )
+       CALL substi_1dd( ar, tri_for_1d )
 
     CONTAINS
 
 #if defined( __intel11 )
-       SUBROUTINE maketri_1dd( j, tri )
+       SUBROUTINE maketri_1dd( j, tri_for_1d )
 #else
        SUBROUTINE maketri_1dd( j )
@@ -473,5 +474,5 @@
 
 #if defined( __intel11 )
-          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
+          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d
 #endif
 
@@ -502,15 +503,15 @@
                 a = -1.0 * ddzu_pres(k+2) * ddzw(k+1)
                 c = -1.0 * ddzu_pres(k+1) * ddzw(k+1)
-                tri(1,i,k) = a + c - l(i)
+                tri_for_1d(1,i,k) = a + c - l(i)
              ENDDO
           ENDDO
           IF ( ibc_p_b == 1 )  THEN
              DO  i = 0, nx
-                tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
+                tri_for_1d(1,i,0) = tri_for_1d(1,i,0) + tri_for_1d(2,i,0)
              ENDDO
           ENDIF
           IF ( ibc_p_t == 1 )  THEN
              DO  i = 0, nx
-                tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
+                tri_for_1d(1,i,nz-1) = tri_for_1d(1,i,nz-1) + tri_for_1d(3,i,nz-1)
              ENDDO
           ENDIF
@@ -520,5 +521,5 @@
 
 #if defined( __intel11 )
-       SUBROUTINE split_1dd( tri )
+       SUBROUTINE split_1dd( tri_for_1d )
 #else
        SUBROUTINE split_1dd
@@ -534,5 +535,5 @@
 
 #if defined( __intel11 )
-          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
+          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d
 #endif
 
@@ -541,10 +542,10 @@
 !--       Splitting
           DO  i = 0, nx
-             tri(4,i,0) = tri(1,i,0)
+             tri_for_1d(4,i,0) = tri_for_1d(1,i,0)
           ENDDO
           DO  k = 1, nz-1
              DO  i = 0, nx
-                tri(5,i,k) = tri(2,i,k) / tri(4,i,k-1)
-                tri(4,i,k) = tri(1,i,k) - tri(3,i,k-1) * tri(5,i,k)
+                tri_for_1d(5,i,k) = tri_for_1d(2,i,k) / tri_for_1d(4,i,k-1)
+                tri_for_1d(4,i,k) = tri_for_1d(1,i,k) - tri_for_1d(3,i,k-1) * tri_for_1d(5,i,k)
              ENDDO
           ENDDO
@@ -553,5 +554,5 @@
 
 
-       SUBROUTINE substi_1dd( ar, tri )
+       SUBROUTINE substi_1dd( ar, tri_for_1d )
 
 !------------------------------------------------------------------------------!
@@ -565,5 +566,5 @@
           REAL, DIMENSION(0:nx,nz)       ::  ar
           REAL, DIMENSION(0:nx,0:nz-1)   ::  ar1
-          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
+          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d
 
 !
@@ -574,5 +575,5 @@
           DO  k = 1, nz-1
              DO  i = 0, nx
-                ar1(i,k) = ar(i,k+1) - tri(5,i,k) * ar1(i,k-1)
+                ar1(i,k) = ar(i,k+1) - tri_for_1d(5,i,k) * ar1(i,k-1)
              ENDDO
           ENDDO
@@ -584,10 +585,10 @@
 !--       the model domain.
           DO  i = 0, nx
-             ar(i,nz) = ar1(i,nz-1) / ( tri(4,i,nz-1) + 1.0E-20 )
+             ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20 )
           ENDDO
           DO  k = nz-2, 0, -1
              DO  i = 0, nx
-                ar(i,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,k+2) ) &
-                            / tri(4,i,k)
+                ar(i,k+1) = ( ar1(i,k) - tri_for_1d(3,i,k) * ar(i,k+2) ) &
+                            / tri_for_1d(4,i,k)
              ENDDO
           ENDDO
