Index: /palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- /palm/trunk/SOURCE/boundary_conds.f90	(revision 1112)
+++ /palm/trunk/SOURCE/boundary_conds.f90	(revision 1113)
@@ -1,3 +1,3 @@
- SUBROUTINE boundary_conds( range )
+ SUBROUTINE boundary_conds
 
 !--------------------------------------------------------------------------------!
@@ -20,8 +20,12 @@
 ! Current revisions:
 ! -----------------
+! GPU-porting
+! dummy argument "range" removed
+! Bugfix: wrong index in loops of radiation boundary condition
 !  
 ! Former revisions:
 ! -----------------
 ! $Id$
+!
 !
 ! 1053 2012-11-13 17:11:03Z hoffmann
@@ -104,6 +108,4 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  range
-
     INTEGER ::  i, j, k
 
@@ -111,100 +113,142 @@
 
 
-    IF ( range == 'main')  THEN
-!
-!--    Bottom boundary 
-       IF ( ibc_uv_b == 1 )  THEN
-          u_p(nzb,:,:) = u_p(nzb+1,:,:)
-          v_p(nzb,:,:) = v_p(nzb+1,:,:)
-       ENDIF
+!
+!-- Bottom boundary 
+    IF ( ibc_uv_b == 1 )  THEN
+       !$acc kernels present( u_p, v_p )
+       u_p(nzb,:,:) = u_p(nzb+1,:,:)
+       v_p(nzb,:,:) = v_p(nzb+1,:,:)
+       !$acc end kernels
+    ENDIF
+
+    !$acc kernels present( nzb_w_inner, w_p )
+    DO  i = nxlg, nxrg
+       DO  j = nysg, nyng
+          w_p(nzb_w_inner(j,i),j,i) = 0.0
+       ENDDO
+    ENDDO
+    !$acc end kernels
+
+!
+!-- Top boundary
+    IF ( ibc_uv_t == 0 )  THEN
+       !$acc kernels present( u_init, u_p, v_init, v_p )
+        u_p(nzt+1,:,:) = u_init(nzt+1)
+        v_p(nzt+1,:,:) = v_init(nzt+1)
+       !$acc end kernels
+    ELSE
+       !$acc kernels present( u_p, v_p )
+        u_p(nzt+1,:,:) = u_p(nzt,:,:)
+        v_p(nzt+1,:,:) = v_p(nzt,:,:)
+       !$acc end kernels
+    ENDIF
+    !$acc kernels present( w_p )
+    w_p(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
+    !$acc end kernels
+
+!
+!-- Temperature at bottom boundary.
+!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
+!-- the sea surface temperature of the coupled ocean model.
+    IF ( ibc_pt_b == 0 )  THEN
+       !$acc kernels present( nzb_s_inner, pt, pt_p )
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
-             w_p(nzb_w_inner(j,i),j,i) = 0.0
+             pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
           ENDDO
        ENDDO
-
-!
-!--    Top boundary
-       IF ( ibc_uv_t == 0 )  THEN
-           u_p(nzt+1,:,:) = u_init(nzt+1)
-           v_p(nzt+1,:,:) = v_init(nzt+1)
-       ELSE
-           u_p(nzt+1,:,:) = u_p(nzt,:,:)
-           v_p(nzt+1,:,:) = v_p(nzt,:,:)
+       !$acc end kernels
+    ELSEIF ( ibc_pt_b == 1 )  THEN
+       !$acc kernels present( nzb_s_inner, pt_p )
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
+          ENDDO
+       ENDDO
+      !$acc end kernels
+    ENDIF
+
+!
+!-- Temperature at top boundary
+    IF ( ibc_pt_t == 0 )  THEN
+       !$acc kernels present( pt, pt_p )
+        pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
+       !$acc end kernels
+    ELSEIF ( ibc_pt_t == 1 )  THEN
+       !$acc kernels present( pt_p )
+        pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
+       !$acc end kernels
+    ELSEIF ( ibc_pt_t == 2 )  THEN
+       !$acc kernels present( dzu, pt_p )
+        pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
+       !$acc end kernels
+    ENDIF
+
+!
+!-- Boundary conditions for TKE
+!-- Generally Neumann conditions with de/dz=0 are assumed
+    IF ( .NOT. constant_diffusion )  THEN
+       !$acc kernels present( e_p, nzb_s_inner )
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
+          ENDDO
+       ENDDO
+       e_p(nzt+1,:,:) = e_p(nzt,:,:)
+       !$acc end kernels
+    ENDIF
+
+!
+!-- Boundary conditions for salinity
+    IF ( ocean )  THEN
+!
+!--    Bottom boundary: Neumann condition because salinity flux is always
+!--    given
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i)
+          ENDDO
+       ENDDO
+
+!
+!--    Top boundary: Dirichlet or Neumann
+       IF ( ibc_sa_t == 0 )  THEN
+           sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
+       ELSEIF ( ibc_sa_t == 1 )  THEN
+           sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
        ENDIF
-       w_p(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
-
-!
-!--    Temperature at bottom boundary.
-!--    In case of coupled runs (ibc_pt_b = 2) the temperature is given by
-!--    the sea surface temperature of the coupled ocean model.
-       IF ( ibc_pt_b == 0 )  THEN
+
+    ENDIF
+
+!
+!-- Boundary conditions for total water content or scalar,
+!-- bottom and top boundary (see also temperature)
+    IF ( humidity  .OR.  passive_scalar )  THEN
+!
+!--    Surface conditions for constant_humidity_flux
+       IF ( ibc_q_b == 0 ) THEN
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
-             ENDDO
-          ENDDO
-       ELSEIF ( ibc_pt_b == 1 )  THEN
+                q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
+             ENDDO
+          ENDDO
+       ELSE
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
+                q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
              ENDDO
           ENDDO
        ENDIF
-
-!
-!--    Temperature at top boundary
-       IF ( ibc_pt_t == 0 )  THEN
-           pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
-       ELSEIF ( ibc_pt_t == 1 )  THEN
-           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
-       ELSEIF ( ibc_pt_t == 2 )  THEN
-           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
-       ENDIF
-
-!
-!--    Boundary conditions for TKE
-!--    Generally Neumann conditions with de/dz=0 are assumed
-       IF ( .NOT. constant_diffusion )  THEN
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-                e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
-             ENDDO
-          ENDDO
-          e_p(nzt+1,:,:) = e_p(nzt,:,:)
-       ENDIF
-
-!
-!--    Boundary conditions for salinity
-       IF ( ocean )  THEN
-!
-!--       Bottom boundary: Neumann condition because salinity flux is always
-!--       given
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-                sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i)
-             ENDDO
-          ENDDO
-
-!
-!--       Top boundary: Dirichlet or Neumann
-          IF ( ibc_sa_t == 0 )  THEN
-              sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
-          ELSEIF ( ibc_sa_t == 1 )  THEN
-              sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
-          ENDIF
-
-       ENDIF
-
-!
-!--    Boundary conditions for total water content or scalar,
-!--    bottom and top boundary (see also temperature)
-       IF ( humidity  .OR.  passive_scalar )  THEN
-!
+!
+!--    Top boundary
+       q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
+
+       IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+!             
 !--       Surface conditions for constant_humidity_flux
-          IF ( ibc_q_b == 0 ) THEN
+          IF ( ibc_qr_b == 0 ) THEN
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
+                   qr_p(nzb_s_inner(j,i),j,i) = qr(nzb_s_inner(j,i),j,i)
                 ENDDO
              ENDDO
@@ -212,5 +256,5 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
+                   qr_p(nzb_s_inner(j,i),j,i) = qr_p(nzb_s_inner(j,i)+1,j,i)
                 ENDDO
              ENDDO
@@ -218,45 +262,23 @@
 !
 !--       Top boundary
-          q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
-
-          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+          qr_p(nzt+1,:,:) = qr_p(nzt,:,:) + bc_qr_t_val * dzu(nzt+1)
 !             
-!--          Surface conditions for constant_humidity_flux
-             IF ( ibc_qr_b == 0 ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      qr_p(nzb_s_inner(j,i),j,i) = qr(nzb_s_inner(j,i),j,i)
-                   ENDDO
+!--       Surface conditions for constant_humidity_flux
+          IF ( ibc_nr_b == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   nr_p(nzb_s_inner(j,i),j,i) = nr(nzb_s_inner(j,i),j,i)
                 ENDDO
-             ELSE
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      qr_p(nzb_s_inner(j,i),j,i) = qr_p(nzb_s_inner(j,i)+1,j,i)
-                   ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   nr_p(nzb_s_inner(j,i),j,i) = nr_p(nzb_s_inner(j,i)+1,j,i)
                 ENDDO
-             ENDIF
-!
-!--          Top boundary
-             qr_p(nzt+1,:,:) = qr_p(nzt,:,:) + bc_qr_t_val * dzu(nzt+1)
-!             
-!--          Surface conditions for constant_humidity_flux
-             IF ( ibc_nr_b == 0 ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      nr_p(nzb_s_inner(j,i),j,i) = nr(nzb_s_inner(j,i),j,i)
-                   ENDDO
-                ENDDO
-             ELSE
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      nr_p(nzb_s_inner(j,i),j,i) = nr_p(nzb_s_inner(j,i)+1,j,i)
-                   ENDDO
-                ENDDO
-             ENDIF
-!
-!--          Top boundary
-             nr_p(nzt+1,:,:) = nr_p(nzt,:,:) + bc_nr_t_val * dzu(nzt+1)
-          ENDIF
-
+             ENDDO
+          ENDIF
+!
+!--       Top boundary
+          nr_p(nzt+1,:,:) = nr_p(nzt,:,:) + bc_nr_t_val * dzu(nzt+1)
        ENDIF
 
@@ -707,5 +729,5 @@
 !--       Calculate the new velocities
           DO  k = nzb+1, nzt+1
-             DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) *            &
                                        ( u(k,j,0) - u(k,j,1) ) * ddx
@@ -846,5 +868,5 @@
 !--       Calculate the new velocities
           DO  k = nzb+1, nzt+1
-             DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) *      &
                                        ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
Index: /palm/trunk/SOURCE/exchange_horiz.f90
===================================================================
--- /palm/trunk/SOURCE/exchange_horiz.f90	(revision 1112)
+++ /palm/trunk/SOURCE/exchange_horiz.f90	(revision 1113)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! GPU-porting for single-core (1PE) mode
 !
 ! Former revisions:
@@ -82,5 +82,5 @@
 
 
-    INTEGER ::  nbgp_local
+    INTEGER ::  i, j, k, nbgp_local
     REAL, DIMENSION(nzb:nzt+1,nys-nbgp_local:nyn+nbgp_local, &
                     nxl-nbgp_local:nxr+nbgp_local) ::  ar
@@ -199,13 +199,47 @@
 
 !
-!-- Lateral boundary conditions in the non-parallel case
+!-- Lateral boundary conditions in the non-parallel case.
+!-- Case dependent, because in GPU mode still not all arrays are on device. This
+!-- workaround has to be removed later. Also, since PGI compiler 12.5 has problems
+!-- with array syntax, explicit loops are used.
     IF ( bc_lr == 'cyclic' )  THEN
-        ar(:,:,nxl-nbgp_local:nxl-1) = ar(:,:,nxr-nbgp_local+1:nxr)
-        ar(:,:,nxr+1:nxr+nbgp_local) = ar(:,:,nxl:nxl+nbgp_local-1)
+       IF ( on_device )  THEN
+          !$acc kernels present( ar )
+          !$acc loop independent
+          DO  i = 0, nbgp_local-1
+             DO  j = nys-nbgp_local, nyn+nbgp_local
+                !$acc loop vector( 32 )
+                DO  k = nzb, nzt+1
+                   ar(k,j,nxl-nbgp_local+i) = ar(k,j,nxr-nbgp_local+1+i)
+                   ar(k,j,nxr+1+i)          = ar(k,j,nxl+i)
+                ENDDO
+             ENDDO
+          ENDDO
+          !$acc end kernels
+       ELSE
+          ar(:,:,nxl-nbgp_local:nxl-1) = ar(:,:,nxr-nbgp_local+1:nxr)
+          ar(:,:,nxr+1:nxr+nbgp_local) = ar(:,:,nxl:nxl+nbgp_local-1)
+       ENDIF
     ENDIF
 
     IF ( bc_ns == 'cyclic' )  THEN
-        ar(:,nys-nbgp_local:nys-1,:) = ar(:,nyn-nbgp_local+1:nyn,:)
-        ar(:,nyn+1:nyn+nbgp_local,:) = ar(:,nys:nys+nbgp_local-1,:)
+       IF ( on_device )  THEN
+          !$acc kernels present( ar )
+          !$acc loop
+          DO  i = nxl-nbgp_local, nxr+nbgp_local
+             !$acc loop independent
+             DO  j = 0, nbgp_local-1
+                !$acc loop vector( 32 )
+                DO  k = nzb, nzt+1
+                   ar(k,nys-nbgp_local+j,i) = ar(k,nyn-nbgp_local+1+j,i)
+                     ar(k,nyn+1+j,i)          = ar(k,nys+j,i)
+                ENDDO
+             ENDDO
+          ENDDO
+          !$acc end kernels
+       ELSE
+          ar(:,nys-nbgp_local:nys-1,:) = ar(:,nyn-nbgp_local+1:nyn,:)
+          ar(:,nyn+1:nyn+nbgp_local,:) = ar(:,nys:nys+nbgp_local-1,:)
+       ENDIF
     ENDIF
 
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 1112)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 1113)
@@ -23,5 +23,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! openACC directive modified
 !
 ! Former revisions:
@@ -1435,5 +1435,5 @@
        CALL disturb_field( nzb_v_inner, tend, v )
        n_sor = nsor_ini
-       !$acc data copy( d, ddzw, nzb_s_inner, tric, u, v, w, tend )
+       !$acc data copy( d, ddzu, ddzw, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, p, tric, u, v, w, weight_pres, weight_substep, tend )
        CALL pres
        !$acc end data
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 1112)
+++ /palm/trunk/SOURCE/modules.f90	(revision 1113)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! +on_device
 !
 ! Former revisions:
@@ -752,4 +752,5 @@
                 masking_method = .FALSE., mg_switch_to_pe0 = .FALSE., &
                 netcdf_output = .FALSE., neutral = .FALSE., ocean = .FALSE., &
+                on_device = .FALSE., &
                 outflow_l = .FALSE., outflow_n = .FALSE., outflow_r = .FALSE., &
                 outflow_s = .FALSE., passive_scalar = .FALSE., &
Index: /palm/trunk/SOURCE/palm.f90
===================================================================
--- /palm/trunk/SOURCE/palm.f90	(revision 1112)
+++ /palm/trunk/SOURCE/palm.f90	(revision 1113)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! openACC statements modified
 !
 ! Former revisions:
@@ -237,11 +237,12 @@
 !-- Declare and initialize variables in the accelerator memory with their
 !-- host values
-    !$acc  data copyin( d, diss, e, e_p, kh, km, pt, pt_p, q, ql, tend, te_m, tpt_m, tu_m, tv_m, tw_m, u, u_p, v, vpt, v_p, w, w_p )          &
-    !$acc       copyin( tric, ddzu, ddzw, dd2zu, l_grid, l_wall, ptdf_x, ptdf_y, pt_init, rdf, rdf_sc, ug, vg, zu, zw )   &
+    !$acc  data copyin( d, diss, e, e_p, kh, km, p, pt, pt_p, q, ql, tend, te_m, tpt_m, tu_m, tv_m, tw_m, u, u_p, v, vpt, v_p, w, w_p )          &
+    !$acc       copyin( tric, dzu, ddzu, ddzw, dd2zu, l_grid, l_wall, ptdf_x, ptdf_y, pt_init, rdf, rdf_sc, ug, u_init, vg, v_init, zu, zw )   &
     !$acc       copyin( hom, qs, qsws, qswst, rif, rif_wall, shf, ts, tswst, us, usws, uswst, vsws, vswst, z0, z0h )      &
     !$acc       copyin( fxm, fxp, fym, fyp, fwxm, fwxp, fwym, fwyp, nzb_diff_s_inner, nzb_diff_s_outer, nzb_diff_u )       &
     !$acc       copyin( nzb_diff_v, nzb_s_inner, nzb_s_outer, nzb_u_inner )    &
     !$acc       copyin( nzb_u_outer, nzb_v_inner, nzb_v_outer, nzb_w_inner )   &
-    !$acc       copyin( nzb_w_outer, wall_heatflux, wall_e_x, wall_e_y, wall_u, wall_v, wall_w_x, wall_w_y, wall_flags_0 )
+    !$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( weight_pres, weight_substep )
 !
 !-- Integration of the model equations using timestep-scheme
Index: /palm/trunk/SOURCE/pres.f90
===================================================================
--- /palm/trunk/SOURCE/pres.f90	(revision 1112)
+++ /palm/trunk/SOURCE/pres.f90	(revision 1113)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! GPU-porting of several loops, some loops rearranged
 !
 ! Former revisions:
@@ -403,5 +403,5 @@
 !--       Solver for 2d-decomposition
           CALL poisfft( d, tend )
-          !$acc update host( d )
+
        ELSEIF ( psolver == 'poisfft_hybrid' )  THEN 
 !
@@ -410,4 +410,5 @@
 !--       are some optimization problems in poisfft
           CALL poisfft_hybrid( d )
+
        ENDIF
 
@@ -416,6 +417,9 @@
 !--    z-direction
        !$OMP PARALLEL DO
+       !$acc kernels present( d, tend )
+       !$acc loop
        DO  i = nxl, nxr
           DO  j = nys, nyn
+             !$acc loop vector( 32 )
              DO  k = nzb+1, nzt
                 tend(k,j,i) = d(k,j,i)
@@ -423,4 +427,5 @@
           ENDDO
        ENDDO
+       !$acc end kernels
 
 !
@@ -432,4 +437,5 @@
 !--       Neumann (dp/dz = 0)
           !$OMP PARALLEL DO
+          !$acc kernels present( nzb_s_inner, tend )
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
@@ -437,4 +443,5 @@
              ENDDO
           ENDDO
+          !$acc end kernels
 
        ELSE
@@ -442,4 +449,5 @@
 !--       Dirichlet
           !$OMP PARALLEL DO
+          !$acc kernels present( tend )
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
@@ -447,4 +455,5 @@
              ENDDO
           ENDDO
+          !$acc end kernels
 
        ENDIF
@@ -456,4 +465,5 @@
 !--       Neumann
           !$OMP PARALLEL DO
+          !$acc kernels present( tend )
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
@@ -461,4 +471,5 @@
              ENDDO
           ENDDO
+          !$acc end kernels
 
        ELSE
@@ -466,4 +477,5 @@
 !--       Dirichlet
           !$OMP PARALLEL DO
+          !$acc kernels present( tend )
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
@@ -471,4 +483,5 @@
              ENDDO
           ENDDO
+          !$acc end kernels
 
        ENDIF
@@ -476,5 +489,15 @@
 !
 !--    Exchange boundaries for p
+       IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
+          on_device = .TRUE.         ! to be removed after complete porting
+       ELSE                          ! of ghost point exchange
+          !$acc update host( tend )
+       ENDIF
        CALL exchange_horiz( tend, nbgp )
+       IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
+          on_device = .FALSE.        ! to be removed after complete porting
+       ELSE                          ! of ghost point exchange
+          !$acc update device( tend )
+       ENDIF
       
     ELSEIF ( psolver == 'sor' )  THEN
@@ -531,6 +554,9 @@
        !$OMP PARALLEL PRIVATE (i,j,k)
        !$OMP DO
+       !$acc kernels present( p, tend, weight_substep )
+       !$acc loop
        DO  i = nxl-1, nxr+1
           DO  j = nys-1, nyn+1
+             !$acc loop vector( 32 )
              DO  k = nzb, nzt+1
                 p(k,j,i) = tend(k,j,i) * &
@@ -539,4 +565,5 @@
           ENDDO
        ENDDO
+       !$acc end kernels
        !$OMP END PARALLEL
 
@@ -544,6 +571,9 @@
        !$OMP PARALLEL PRIVATE (i,j,k)
        !$OMP DO
+       !$acc kernels present( p, tend, weight_substep )
+       !$acc loop
        DO  i = nxl-1, nxr+1
           DO  j = nys-1, nyn+1
+             !$acc loop vector( 32 )
              DO  k = nzb, nzt+1
                 p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
@@ -552,4 +582,5 @@
           ENDDO
        ENDDO
+       !$acc end kernels
        !$OMP END PARALLEL
 
@@ -571,43 +602,71 @@
     !$OMP PARALLEL PRIVATE (i,j,k)
     !$OMP DO
+    !$acc kernels present( ddzu, nzb_u_inner, nzb_v_inner, nzb_w_inner, tend, u, v, w, weight_pres )
+    !$acc loop
     DO  i = nxl, nxr   
        DO  j = nys, nyn
-          DO  k = nzb_w_inner(j,i)+1, nzt
-             w(k,j,i) = w(k,j,i) - dt_3d *                                 &
-                        ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) *      &
-                        weight_pres(intermediate_timestep_count)
-          ENDDO
-          DO  k = nzb_u_inner(j,i)+1, nzt
-             u(k,j,i) = u(k,j,i) - dt_3d *                                 &
-                        ( tend(k,j,i) - tend(k,j,i-1) ) * ddx *            &
-                        weight_pres(intermediate_timestep_count)
-          ENDDO
-          DO  k = nzb_v_inner(j,i)+1, nzt
-             v(k,j,i) = v(k,j,i) - dt_3d *                                 &
-                        ( tend(k,j,i) - tend(k,j-1,i) ) * ddy *            &
-                        weight_pres(intermediate_timestep_count)
+          !$acc loop vector( 32 )
+          DO  k = 1, nzt
+             IF ( k > nzb_w_inner(j,i) )  THEN
+                w(k,j,i) = w(k,j,i) - dt_3d *                                 &
+                           ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) *      &
+                           weight_pres(intermediate_timestep_count)
+             ENDIF
+          ENDDO
+          !$acc loop vector( 32 )
+          DO  k = 1, nzt
+             IF ( k > nzb_u_inner(j,i) )  THEN
+                u(k,j,i) = u(k,j,i) - dt_3d *                                 &
+                           ( tend(k,j,i) - tend(k,j,i-1) ) * ddx *            &
+                           weight_pres(intermediate_timestep_count)
+             ENDIF
+          ENDDO
+          !$acc loop vector( 32 )
+          DO  k = 1, nzt
+             IF ( k > nzb_v_inner(j,i) )  THEN
+                v(k,j,i) = v(k,j,i) - dt_3d *                                 &
+                           ( tend(k,j,i) - tend(k,j-1,i) ) * ddy *            &
+                           weight_pres(intermediate_timestep_count)
+             ENDIF
           ENDDO                                                         
-!
-!--       Sum up the volume flow through the right and north boundary
-          IF ( conserve_volume_flow  .AND.  bc_lr_cyc  .AND.  bc_ns_cyc  .AND. &
-               i == nx )  THEN
-             !$OMP CRITICAL
-             DO  k = nzb_2d(j,i) + 1, nzt
-                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
-             ENDDO
-             !$OMP END CRITICAL
-          ENDIF
-          IF ( conserve_volume_flow  .AND.  bc_ns_cyc  .AND.  bc_lr_cyc  .AND. &
-               j == ny )  THEN
-             !$OMP CRITICAL
-             DO  k = nzb_2d(j,i) + 1, nzt
-                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
-             ENDDO
-             !$OMP END CRITICAL
-          ENDIF
 
        ENDDO
     ENDDO
+    !$acc end kernels
     !$OMP END PARALLEL
+
+!
+!-- Sum up the volume flow through the right and north boundary
+    IF ( conserve_volume_flow  .AND.  bc_lr_cyc  .AND.  bc_ns_cyc  .AND.  &
+         nxr == nx )  THEN
+
+       !$OMP PARALLEL PRIVATE (j,k)
+       !$OMP DO
+       DO  j = nys, nyn
+          !$OMP CRITICAL
+          DO  k = nzb_2d(j,nx) + 1, nzt
+             volume_flow_l(1) = volume_flow_l(1) + u(k,j,nx) * dzw(k)
+          ENDDO
+          !$OMP END CRITICAL
+       ENDDO
+       !$OMP END PARALLEL
+
+    ENDIF
+
+    IF ( conserve_volume_flow  .AND.  bc_ns_cyc  .AND.  bc_lr_cyc  .AND.  &
+         nyn == ny )  THEN
+
+       !$OMP PARALLEL PRIVATE (i,k)
+       !$OMP DO
+       DO  i = nxl, nxr
+          !$OMP CRITICAL
+          DO  k = nzb_2d(ny,i) + 1, nzt
+             volume_flow_l(2) = volume_flow_l(2) + v(k,ny,i) * dzw(k)
+           ENDDO
+          !$OMP END CRITICAL
+       ENDDO
+       !$OMP END PARALLEL
+
+    ENDIF
     
 !
@@ -645,7 +704,17 @@
 !
 !-- Exchange of boundaries for the velocities
+    IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
+       on_device = .TRUE.         ! to be removed after complete porting
+    ELSE                          ! of ghost point exchange
+       !$acc update host( u, v, w )
+    ENDIF
     CALL exchange_horiz( u, nbgp )
     CALL exchange_horiz( v, nbgp )
     CALL exchange_horiz( w, nbgp )
+    IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
+       on_device = .FALSE.        ! to be removed after complete porting
+    ELSE                          ! of ghost point exchange
+       !$acc update device( u, v, w )
+    ENDIF
 
 !
@@ -679,10 +748,26 @@
     ENDDO
 #else
+    !$acc kernels present( d, ddzw, nzb_s_inner, u, v, w )
+    !$acc loop
     DO  i = nxl, nxr
        DO  j = nys, nyn
+          !$acc loop vector( 32 )
+          DO  k = 1, nzt
+             IF ( k > nzb_s_inner(j,i) )  THEN
+                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
+          ENDDO
+       ENDDO
+    ENDDO
+    !$acc end kernels
+!
+!-- Compute possible PE-sum of divergences for flow_statistics
+    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
+    !$OMP DO SCHEDULE( STATIC )
+    DO  i = nxl, nxr
+       DO  j = nys, nyn
           DO  k = nzb_s_inner(j,i)+1, nzt
-             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)
              threadsum = threadsum + ABS( d(k,j,i) )
           ENDDO
@@ -701,8 +786,5 @@
     CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
 
-    !$acc update device( u, v, w )
-
     CALL cpu_log( log_point(8), 'pres', 'stop' )
-    
 
 
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 1112)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 1113)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! ------------------
-! 
+! GPU-porting of boundary conditions,
+! openACC directives updated
+! formal parameter removed from routine boundary_conds
 !
 ! Former revisions:
@@ -242,5 +244,11 @@
 !--       Exchange of ghost points (lateral boundary conditions)
           CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
-          !$acc update host( e_p, pt_p, u_p, v_p, w_p )
+
+          IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
+             on_device = .TRUE.         ! to be removed after complete porting
+          ELSE                          ! of ghost point exchange
+             !$acc update host( e_p, pt_p, u_p, v_p, w_p )
+          ENDIF
+
           CALL exchange_horiz( u_p, nbgp )
           CALL exchange_horiz( v_p, nbgp )
@@ -268,4 +276,10 @@
           IF ( wang_kernel  .OR.  turbulence )  CALL exchange_horiz( diss, nbgp )
 
+          IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
+             on_device = .FALSE.        ! to be removed after complete porting
+          ELSE                          ! of ghost point exchange
+             !$acc update device( e_p, pt_p, u_p, v_p, w_p )
+          ENDIF
+
           CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
 
@@ -273,9 +287,8 @@
 !--       Boundary conditions for the prognostic quantities (except of the
 !--       velocities at the outflow in case of a non-cyclic lateral wall)
-          CALL boundary_conds( 'main' )
+          CALL boundary_conds
 
 !
 !--       Swap the time levels in preparation for the next time step.
-          !$acc update device( e_p, pt_p, u_p, v_p, w_p )
           CALL swap_timelevel
 
@@ -305,4 +318,5 @@
              IF ( time_disturb >= dt_disturb )  THEN
                 !$acc update host( u, v )
+                IF ( numprocs == 1 )  on_device = .FALSE.  ! workaround, remove later
                 IF ( hom(nzb+5,1,pr_palm,0) < disturbance_energy_limit )  THEN
                    CALL disturb_field( nzb_u_inner, tend, u )
@@ -317,4 +331,5 @@
                    dist_range = 0
                 ENDIF
+                IF ( numprocs == 1 )  on_device = .TRUE.  ! workaround, remove later
                 !$acc update device( u, v )
                 time_disturb = time_disturb - dt_disturb
