Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 1109)
+++ palm/trunk/SOURCE/Makefile	(revision 1111)
@@ -20,5 +20,6 @@
 # Current revisions:
 # ------------------
-# 
+# dependencies removed from init_pegrid
+# bugfix: dependency added for cuda_fft_interfaces
 #
 # Former revisions:
@@ -267,5 +268,5 @@
 cpu_log.o: modules.o
 cpu_statistics.o: modules.o
-cuda_fft_interfaces.o: cuda_fft_interfaces.f90
+cuda_fft_interfaces.o: cuda_fft_interfaces.f90 modules.o
 data_log.o: modules.o
 data_output_dvrp.o: modules.o
@@ -303,5 +304,5 @@
 init_masks.o: modules.o
 init_ocean.o: modules.o eqn_state_seawater.o
-init_pegrid.o: modules.o fft_xy.o poisfft.o poisfft_hybrid.o
+init_pegrid.o: modules.o
 init_pt_anomaly.o: modules.o
 init_rankine.o: modules.o
Index: palm/trunk/SOURCE/Makefile_check
===================================================================
--- palm/trunk/SOURCE/Makefile_check	(revision 1109)
+++ palm/trunk/SOURCE/Makefile_check	(revision 1111)
@@ -20,5 +20,5 @@
 # Current revisions:
 # ------------------
-# 
+# dependencies removed from init_pegrid
 #
 # Former revisions:
@@ -134,5 +134,5 @@
 init_grid.o: modules.o
 init_masks.o: modules.o
-init_pegrid.o: modules.o fft_xy.o poisfft.o poisfft_hybrid.o
+init_pegrid.o: modules.o
 local_stop.o: modules.o
 message.o: modules.o
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 1109)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 1111)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! ibc_p_b = 2 removed
 !
 ! Former revisions:
@@ -1601,6 +1601,4 @@
     ELSEIF ( bc_p_b == 'neumann' )  THEN
        ibc_p_b = 1
-    ELSEIF ( bc_p_b == 'neumann+inhomo' )  THEN
-       ibc_p_b = 2
     ELSE
        message_string = 'unknown boundary condition: bc_p_b = "' // &
@@ -1608,9 +1606,5 @@
        CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( ibc_p_b == 2  .AND.  .NOT. prandtl_layer )  THEN
-       message_string = 'boundary condition: bc_p_b = "' // TRIM( bc_p_b ) // &
-                        '" not allowed with prandtl_layer = .FALSE.'
-       CALL message( 'check_parameters', 'PA0060', 1, 2, 0, 6, 0 )
-    ENDIF
+
     IF ( bc_p_t == 'dirichlet' )  THEN
        ibc_p_t = 0
Index: palm/trunk/SOURCE/cpu_statistics.f90
===================================================================
--- palm/trunk/SOURCE/cpu_statistics.f90	(revision 1109)
+++ palm/trunk/SOURCE/cpu_statistics.f90	(revision 1111)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! output of grid point numbers and average CPU time per grid point and timestep
 !
 ! Former revisions:
@@ -69,11 +69,13 @@
 !------------------------------------------------------------------------------!
 
+    USE control_parameters
     USE cpulog
+    USE indices,  ONLY: nx, ny, nz
     USE pegrid
-    USE control_parameters
 
     IMPLICIT NONE
 
     INTEGER    ::  i, ii(1), iii, sender
+    REAL       ::  average_cputime
     REAL, SAVE ::  norm = 1.0
     REAL, DIMENSION(:),   ALLOCATABLE ::  pe_max, pe_min, pe_rms, sum
@@ -157,16 +159,29 @@
 
 !
+!--    Get total time in order to calculate CPU-time per gridpoint and timestep
+       IF ( nr_timesteps_this_run /= 0 )  THEN
+          average_cputime = log_point(1)%sum / REAL( (nx+1) * (ny+1) * nz ) / &
+                            REAL( nr_timesteps_this_run ) * 1E6  ! in micro-sec
+       ELSE
+          average_cputime = -1.0
+       ENDIF
+
+!
 !--    Write cpu-times sorted by size
        CALL check_open( 18 )
 #if defined( __parallel )
-       WRITE ( 18, 100 )  TRIM( run_description_header ),        &
-                          numprocs * threads_per_task, pdims(1), pdims(2), &
-                          threads_per_task
+       WRITE ( 18, 100 )  TRIM( run_description_header ),                          &
+                          numprocs * threads_per_task, pdims(1), pdims(2),         &
+                          threads_per_task, nx+1, ny+1, nz, nr_timesteps_this_run, &
+                          average_cputime
+                          
        IF ( num_acc_per_node /= 0 )  WRITE ( 18, 108 )  num_acc_per_node
        WRITE ( 18, 110 )
 #else
-       WRITE ( 18, 100 )  TRIM( run_description_header ),        &
-                          numprocs * threads_per_task, 1, 1, &
-                          threads_per_task
+       WRITE ( 18, 100 )  TRIM( run_description_header ),                          &
+                          numprocs * threads_per_task, 1, 1,                       &
+                          threads_per_task, nx+1, ny+1, nz, nr_timesteps_this_run, &
+                          average_cputime
+
        IF ( num_acc_per_node /= 0 )  WRITE ( 18, 109 )  num_acc_per_node
        WRITE ( 18, 110 )
@@ -308,5 +323,8 @@
 
 100 FORMAT (A/11('-')//'CPU measures for ',I5,' PEs (',I5,'(x) * ',I5,'(y', &
-            &') tasks *',I5,' threads):')
+            &') tasks *',I5,' threads):'//                                  &
+            'gridpoints (x/y/z): ',20X,I5,' * ',I5,' * ',I5/                &
+            'nr of timesteps: ',22X,I6/                                     &
+            'cpu time per grid point and timestep: ',5X,F8.5,' * 10**-6 s')
 
 101 FORMAT (/'special measures:'/ &
@@ -320,6 +338,6 @@
 106 FORMAT (/'Exchange of ghostpoints via MPI_ISEND/MPI_IRECV')
 107 FORMAT (//)
-108 FORMAT ('Accelerator boards per node: ',I2)
-109 FORMAT ('Accelerator boards: ',I2)
+108 FORMAT ('Accelerator boards per node: ',14X,I2)
+109 FORMAT ('Accelerator boards: ',23X,I2)
 110 FORMAT ('----------------------------------------------------------',   &
             &'------------'//&
Index: palm/trunk/SOURCE/cuda_fft_interfaces.f90
===================================================================
--- palm/trunk/SOURCE/cuda_fft_interfaces.f90	(revision 1109)
+++ palm/trunk/SOURCE/cuda_fft_interfaces.f90	(revision 1111)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! idata and odata changed from 1d- to 3d-arrays
 !
 ! Former revisions:
@@ -87,6 +87,6 @@
 
           INTEGER(C_INT), value ::  plan
-          COMPLEX(dpk), device  ::  idata(*)
-          REAL(dpk), device     ::  odata(*)
+          COMPLEX(dpk), device  ::  idata(:,:,:)
+          REAL(dpk), device     ::  odata(:,:,:)
 
        END SUBROUTINE CUFFTEXECZ2D
@@ -103,6 +103,6 @@
 
           INTEGER(C_INT), value ::  plan
-          REAL(dpk), device     ::  idata(*)
-          COMPLEX(dpk), device  ::  odata(*)
+          REAL(dpk), device     ::  idata(:,:,:)
+          COMPLEX(dpk), device  ::  odata(:,:,:)
 
        END SUBROUTINE CUFFTEXECD2Z
Index: palm/trunk/SOURCE/fft_xy.f90
===================================================================
--- palm/trunk/SOURCE/fft_xy.f90	(revision 1109)
+++ palm/trunk/SOURCE/fft_xy.f90	(revision 1111)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! -----------------
-! 
+! further openACC statements added, CUDA branch completely runs on GPU
+! bugfix: CUDA fft plans adjusted for domain decomposition (before they always
+! used total domain)
 !
 ! Former revisions:
@@ -213,8 +215,8 @@
           total_points_x_transpo = (nx+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1)
           total_points_y_transpo = (ny+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1)
-          CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (ny+1)*nz )
-          CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (ny+1)*nz )
-          CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nx+1)*nz )
-          CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nx+1)*nz )
+          CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
+          CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
+          CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
+          CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
 #else
           message_string = 'no system-specific fft-call available'
@@ -261,5 +263,5 @@
 
        CHARACTER (LEN=*) ::  direction
-       INTEGER ::  i, ishape(1), j, k, m
+       INTEGER ::  i, ishape(1), j, k
 
        LOGICAL ::  forward_fft
@@ -273,7 +275,6 @@
        REAL, DIMENSION(6*(nx+1)) ::  work2
 #elif defined( __cuda_fft )
-       REAL(dpk), DEVICE, DIMENSION(:), ALLOCATABLE    ::  cuda_a_device
-       COMPLEX(dpk), DEVICE, DIMENSION(:), ALLOCATABLE ::  cuda_b_device
-       COMPLEX(dpk), DIMENSION(:), ALLOCATABLE         ::  cuda_host
+       !$acc declare create( ar_tmp )
+       COMPLEX(dpk), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) ::  ar_tmp
 #endif
        REAL, DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::  ar
@@ -502,58 +503,53 @@
 #elif defined( __cuda_fft )
 
-          ALLOCATE( cuda_a_device(0:total_points_x_transpo-1) )
-          ALLOCATE( cuda_b_device(0:((nx+1)/2+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) - 1) )
-          ALLOCATE( cuda_host(0:((nx+1)/2+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) - 1) )
-
-          m = 0
-
           IF ( forward_fft )  THEN
 
-             cuda_a_device = ar(0:total_points_x_transpo-1,nys_x,nzb_x)
-
-             CALL CUFFTEXECD2Z( plan_xf, cuda_a_device, cuda_b_device )
-             cuda_host = cuda_b_device
-
+             !$acc data present( ar )
+             CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
+
+             !$acc kernels
+             !$acc loop
              DO  k = nzb_x, nzt_x
                 DO  j = nys_x, nyn_x
 
+                   !$acc loop vector( 32 )
                    DO  i = 0, (nx+1)/2
-                      ar(i,j,k)      = REAL( cuda_host(m+i) )  * dnx
-                   ENDDO
-
+                      ar(i,j,k)      = REAL( ar_tmp(i,j,k) )  * dnx
+                   ENDDO
+
+                   !$acc loop vector( 32 )
                    DO  i = 1, (nx+1)/2 - 1
-                      ar(nx+1-i,j,k) = AIMAG( cuda_host(m+i) ) * dnx
-                   ENDDO
-
-                   m = m + (nx+1)/2 + 1
-
-                ENDDO
-             ENDDO
-
-          ELSE
-
+                      ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
+                   ENDDO
+
+                ENDDO
+             ENDDO
+             !$acc end kernels
+             !$acc end data
+
+          ELSE
+
+             !$acc data present( ar )
+             !$acc kernels
+             !$acc loop
              DO  k = nzb_x, nzt_x
                 DO  j = nys_x, nyn_x
 
-                   cuda_host(m) = CMPLX( ar(0,j,k), 0.0 )
-
+                   ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
+
+                   !$acc loop vector( 32 )
                    DO  i = 1, (nx+1)/2 - 1
-                      cuda_host(m+i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
-                   ENDDO
-                   cuda_host(m+(nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
-
-                   m = m + (nx+1)/2 + 1
-
-                ENDDO
-             ENDDO
-
-             cuda_b_device = cuda_host
-             CALL CUFFTEXECZ2D( plan_xi, cuda_b_device, cuda_a_device )
-
-             ar(0:total_points_x_transpo-1,nys_x,nzb_x) = cuda_a_device
-
-          ENDIF
-
-          DEALLOCATE( cuda_a_device, cuda_b_device, cuda_host )
+                      ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
+                   ENDDO
+                   ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
+
+                ENDDO
+             ENDDO
+             !$acc end kernels
+
+             CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
+             !$acc end data
+
+          ENDIF
 
 #else
@@ -775,5 +771,5 @@
 
        CHARACTER (LEN=*) ::  direction
-       INTEGER ::  i, j, jshape(1), k, m
+       INTEGER ::  i, j, jshape(1), k
 
        LOGICAL ::  forward_fft
@@ -787,7 +783,6 @@
        REAL, DIMENSION(6*(ny+1)) ::  work2
 #elif defined( __cuda_fft )
-       REAL(dpk), DEVICE, DIMENSION(:), ALLOCATABLE    ::  cuda_a_device
-       COMPLEX(dpk), DEVICE, DIMENSION(:), ALLOCATABLE ::  cuda_b_device
-       COMPLEX(dpk), DIMENSION(:), ALLOCATABLE         ::  cuda_host
+       !$acc declare create( ar_tmp )
+       COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) ::  ar_tmp
 #endif
        REAL, DIMENSION(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) ::  ar
@@ -1013,58 +1008,53 @@
 #elif defined( __cuda_fft )
 
-          ALLOCATE( cuda_a_device(0:total_points_y_transpo-1) )
-          ALLOCATE( cuda_b_device(0:((ny+1)/2+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) - 1) )
-          ALLOCATE( cuda_host(0:((ny+1)/2+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) - 1) )
-
-          m = 0
-
           IF ( forward_fft )  THEN
 
-             cuda_a_device = ar(0:total_points_y_transpo-1,nxl_y,nzb_y)
-
-             CALL CUFFTEXECD2Z( plan_yf, cuda_a_device, cuda_b_device )
-             cuda_host = cuda_b_device
-
+             !$acc data present( ar )
+             CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
+
+             !$acc kernels
+             !$acc loop
              DO  k = nzb_y, nzt_y
                 DO  i = nxl_y, nxr_y
 
+                   !$acc loop vector( 32 )
                    DO  j = 0, (ny+1)/2
-                      ar(j,i,k)      = REAL( cuda_host(m+j) )  * dny
-                   ENDDO
-
+                      ar(j,i,k)      = REAL( ar_tmp(j,i,k) )  * dny
+                   ENDDO
+
+                   !$acc loop vector( 32 )
                    DO  j = 1, (ny+1)/2 - 1
-                      ar(ny+1-j,i,k) = AIMAG( cuda_host(m+j) ) * dny
-                   ENDDO
-
-                   m = m + (ny+1)/2 + 1
-
-                ENDDO
-             ENDDO
-
-          ELSE
-
+                      ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
+                   ENDDO
+
+                ENDDO
+             ENDDO
+             !$acc end kernels
+             !$acc end data
+
+          ELSE
+
+             !$acc data present( ar )
+             !$acc kernels
+             !$acc loop
              DO  k = nzb_y, nzt_y
                 DO  i = nxl_y, nxr_y
 
-                   cuda_host(m) = CMPLX( ar(0,i,k), 0.0 )
-
+                   ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
+
+                   !$acc loop vector( 32 )
                    DO  j = 1, (ny+1)/2 - 1
-                      cuda_host(m+j) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
-                   ENDDO
-                   cuda_host(m+(ny+1)/2) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
-
-                   m = m + (ny+1)/2 + 1
-
-                ENDDO
-             ENDDO
-
-             cuda_b_device = cuda_host
-             CALL CUFFTEXECZ2D( plan_yi, cuda_b_device, cuda_a_device )
-
-             ar(0:total_points_y_transpo-1,nxl_y,nzb_y) = cuda_a_device
-
-          ENDIF
-
-          DEALLOCATE( cuda_a_device, cuda_b_device, cuda_host )
+                      ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
+                   ENDDO
+                   ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
+
+                ENDDO
+             ENDDO
+             !$acc end kernels
+
+             CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
+             !$acc end data
+
+          ENDIF
 
 #else
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 1109)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 1111)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! openACC directive added
 !
 ! Former revisions:
@@ -182,4 +183,6 @@
 
     ENDIF
+
+    !$acc update host( km, kh, e, pt, qs, qsws, rif, shf, ts, u, v, w )
 
 !
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 1109)
+++ palm/trunk/SOURCE/header.f90	(revision 1111)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! output of accelerator board information
+! ibc_p_b = 2 removed
 !
 ! Former revisions:
@@ -291,4 +292,5 @@
                           threads_per_task, pdims(1), pdims(2), TRIM( char1 )
     ENDIF
+    IF ( num_acc_per_node /= 0 )  WRITE ( io, 117 )  num_acc_per_node    
     IF ( ( host(1:3) == 'ibm'  .OR.  host(1:3) == 'nec'  .OR.    &
            host(1:2) == 'lc'   .OR.  host(1:3) == 'dec' )  .AND. &
@@ -305,4 +307,6 @@
        WRITE ( io, 108 )  maximum_parallel_io_streams
     ENDIF
+#else
+    IF ( num_acc_per_node /= 0 )  WRITE ( io, 120 )  num_acc_per_node
 #endif
     WRITE ( io, 99 )
@@ -593,6 +597,4 @@
     ELSEIF ( ibc_p_b == 1 )  THEN
        runten = 'p(0)     = p(1)   |'
-    ELSE
-       runten = 'p(0)     = p(1) +R|'
     ENDIF
     IF ( ibc_p_t == 0 )  THEN
@@ -1613,4 +1615,5 @@
             37X,'independent precursor runs'/             &
             37X,42('-'))
+117 FORMAT (' Accelerator boards / node:  ',I2)
 #endif
 110 FORMAT (/' Numerical Schemes:'/ &
@@ -1627,4 +1630,5 @@
             '     translation velocity = ',A/ &
             '     distance advected ',A,':  ',F8.3,' km(x)  ',F8.3,' km(y)')
+120 FORMAT (' Accelerator boards: ',8X,I2)
 122 FORMAT (' --> Time differencing scheme: ',A)
 123 FORMAT (' --> Rayleigh-Damping active, starts ',A,' z = ',F8.2,' m'/ &
@@ -1680,5 +1684,5 @@
              ' CPU-time used:       ',F9.3,' s     per timestep:               ', &
                '  ',F9.3,' s'/                                                    &
-             '                                   per second of simulated tim',    &
+             '                                      per second of simulated tim', &
                'e: ',F9.3,' s')
 207 FORMAT ( ' Coupling start time: ',F9.3,' s')
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 1109)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 1111)
@@ -23,5 +23,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! openACC directives added for pres
+! array diss allocated only if required
 !
 ! Former revisions:
@@ -234,4 +235,5 @@
     USE random_function_mod
     USE statistics
+    USE transpose_indices
 
     IMPLICIT NONE
@@ -334,4 +336,10 @@
     ENDIF
 
+!
+!-- Array for storing constant coeffficients of the tridiagonal solver
+    IF ( psolver == 'poisfft' )  THEN
+       ALLOCATE( tric(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) )
+    ENDIF
+
     IF ( humidity  .OR.  passive_scalar ) THEN
 !
@@ -474,7 +482,4 @@
     IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
        ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-    ELSE
-       ALLOCATE ( diss(2,2,2) )  ! required because diss is used as a
-                                 ! formal parameter
     ENDIF
 
@@ -1427,5 +1432,7 @@
        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 )
        CALL pres
+       !$acc end data
        n_sor = nsor
     ENDIF
Index: palm/trunk/SOURCE/init_pegrid.f90
===================================================================
--- palm/trunk/SOURCE/init_pegrid.f90	(revision 1109)
+++ palm/trunk/SOURCE/init_pegrid.f90	(revision 1111)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! initialization of poisfft moved to poisfft
 !
 ! Former revisions:
@@ -150,10 +150,7 @@
 
     USE control_parameters
-    USE fft_xy
     USE grid_variables
     USE indices
     USE pegrid
-    USE poisfft_mod
-    USE poisfft_hybrid_mod
     USE statistics
     USE transpose_indices
@@ -1135,10 +1132,4 @@
     ENDIF
 
-    IF ( psolver == 'poisfft_hybrid' )  THEN
-       CALL poisfft_hybrid_ini
-    ELSEIF ( psolver == 'poisfft' )  THEN
-       CALL poisfft_init
-    ENDIF
-
 !
 !-- Allocate wall flag arrays used in the multigrid solver
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 1109)
+++ palm/trunk/SOURCE/modules.f90	(revision 1111)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! +tric, nr_timesteps_this_run
 !
 ! Former revisions:
@@ -407,12 +407,12 @@
 
     REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
-          c_u, c_v, c_w, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q, diss_s_qr,  &
-          diss_s_sa, diss_s_u, diss_s_v, diss_s_w, dzu_mg, dzw_mg, flux_s_e,   &
-          flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr, flux_s_sa, flux_s_u,      &
-          flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg, mean_inflow_profiles, nrs,  &
-          nrsws, nrswst, pt_slope_ref, qs, qsws, qswst, qswst_remote, qrs,     &
-          qrsws, qrswst, rif, saswsb, saswst, shf, total_2d_a, total_2d_o, ts, &
-          tswst, us, usws, uswst, vsws, vswst, z0, z0h
-          
+          c_u, c_v, c_w, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q,             &
+          diss_s_qr, diss_s_sa, diss_s_u, diss_s_v, diss_s_w, dzu_mg, dzw_mg,  &
+          flux_s_e, flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr, flux_s_sa,      &
+          flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg,                   &
+          mean_inflow_profiles, nrs, nrsws, nrswst, pt_slope_ref, qs, qsws,    &
+          qswst, qswst_remote, qrs, qrsws, qrswst, rif, saswsb, saswst, shf,   &
+          total_2d_a, total_2d_o, ts, tswst, us, usws, uswst, vsws, vswst, z0, &
+          z0h
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
@@ -422,6 +422,6 @@
           flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km, lad_s,   &
           lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, tend_pt,    &
-          tend_nr, tend_q, tend_qr, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n,  &
-          v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
+          tend_nr, tend_q, tend_qr, tric, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l,   &
+          v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
             
 
@@ -684,6 +684,6 @@
                 maximum_parallel_io_streams = -1, max_pr_user = 0, &
                 mgcycles = 0, mg_cycles = -1, mg_switch_to_pe0_level = 0, mid, &
-                netcdf_data_format = 2, ngsrb = 2, nsor = 20, &
-                nsor_ini = 100, n_sor, normalizing_region = 0, &
+                netcdf_data_format = 2, ngsrb = 2, nr_timesteps_this_run = 0, &
+                nsor = 20, nsor_ini = 100, n_sor, normalizing_region = 0, &
                 nz_do3d = -9999, pch_index = 0, prt_time_count = 0, &
                 recycling_plane, runnr = 0, &
Index: palm/trunk/SOURCE/palm.f90
===================================================================
--- palm/trunk/SOURCE/palm.f90	(revision 1109)
+++ palm/trunk/SOURCE/palm.f90	(revision 1111)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! openACC statements updated
 !
 ! Former revisions:
@@ -234,6 +234,6 @@
 !-- Declare and initialize variables in the accelerator memory with their
 !-- host values
-    !$acc  data copyin( diss, e, e_p, kh, km, pt, pt_p, q, ql, tend, te_m, tpt_m, tu_m, tv_m, tw_m, u, u_p, v, vpt, v_p, w, w_p )          &
-    !$acc       copyin( ddzu, ddzw, dd2zu, l_grid, l_wall, ptdf_x, ptdf_y, pt_init, rdf, rdf_sc, ug, vg, zu, zw )   &
+    !$acc  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       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 )       &
Index: palm/trunk/SOURCE/poisfft.f90
===================================================================
--- palm/trunk/SOURCE/poisfft.f90	(revision 1109)
+++ palm/trunk/SOURCE/poisfft.f90	(revision 1111)
@@ -20,5 +20,15 @@
 ! Current revisions:
 ! -----------------
-! 
+! further openACC porting of non-parallel (MPI) branch:
+! tridiagonal routines split into extermal subroutines (instead using CONTAINS),
+! no distinction between parallel/non-parallel in poisfft and tridia any more,
+! tridia routines moved to end of file because of probable bug in PGI compiler
+! (otherwise "invalid device function" is indicated during runtime),
+! optimization of tridia routines: constant elements and coefficients of tri are
+! stored in seperate arrays ddzuw and tric, last dimension of tri reduced from 5
+! to 2,
+! poisfft_init is now called internally from poisfft, maketri is called from
+! poisfft_init,
+! ibc_p_b = 2 removed
 !
 ! Former revisions:
@@ -157,4 +167,8 @@
     IMPLICIT NONE
 
+    LOGICAL, SAVE ::  poisfft_initialized = .FALSE.
+
+    REAL, DIMENSION(:,:), ALLOCATABLE ::  ddzuw
+
     PRIVATE
 
@@ -181,7 +195,31 @@
     SUBROUTINE poisfft_init
 
+       USE arrays_3d,  ONLY:  ddzu_pres, ddzw
+
+       IMPLICIT NONE
+
+       INTEGER ::  k
+
+
        CALL fft_init
 
+       ALLOCATE( ddzuw(0:nz-1,3) )
+
+       DO  k = 0, nz-1
+          ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1)
+          ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1)
+          ddzuw(k,3) = -1.0 * &
+                       ( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) )
+       ENDDO
+!
+!--    Calculate constant coefficients of the tridiagonal matrix
+#if ! defined ( __check )
+       CALL maketri
+#endif
+
+       poisfft_initialized = .TRUE.
+
     END SUBROUTINE poisfft_init
+
 
 #if ! defined ( __check )
@@ -199,8 +237,9 @@
        CALL cpu_log( log_point_s(3), 'poisfft', 'start' )
 
+       IF ( .NOT. poisfft_initialized )  CALL poisfft_init
+
 !
 !--    Two-dimensional Fourier Transformation in x- and y-direction.
-#if defined( __parallel )
-       IF ( pdims(2) == 1 )  THEN
+       IF ( pdims(2) == 1  .AND.  pdims(1) > 1 )  THEN
 
 !
@@ -217,5 +256,5 @@
           CALL tr_xy_ffty( ar, work, ar )
 
-       ELSEIF ( pdims(1) == 1 )  THEN
+       ELSEIF ( pdims(1) == 1  .AND.  pdims(2) > 1 )  THEN
 
 !
@@ -235,5 +274,5 @@
 
 !
-!--       2d-domain-decomposition
+!--       2d-domain-decomposition or no decomposition (1 PE run)
 !--       Transposition z --> x
           CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
@@ -296,63 +335,4 @@
        ENDIF
 
-#else
-
-!
-!--    Two-dimensional Fourier Transformation along x- and y-direction.
-       CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
-       !$acc data copyin( ar, work )
-       CALL transpose_zx( ar, work, ar )
-       !$acc update host( ar )
-       CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
-       
-
-       CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
-       CALL fft_x( ar, 'forward' )
-       CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
-
-       CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
-       CALL transpose_xy( ar, work, ar )
-       CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
-
-       CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
-       CALL fft_y( ar, 'forward' )
-       CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
-
-!
-!--    Solve the tridiagonal equation system along z
-       CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
-       CALL transpose_yz( ar, work, ar )
-       CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
-
-       CALL cpu_log( log_point_s(6), 'tridia', 'start' )
-       CALL tridia( ar )
-       CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
-
-       CALL cpu_log( log_point_s(8), 'transpo invers', 'start' )
-       CALL transpose_zy( ar, work, ar )
-       CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
-
-!
-!--    Inverse Fourier Transformation.
-       CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
-       CALL fft_y( ar, 'backward' )
-       CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
-
-       CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
-       CALL transpose_yx( ar, work, ar )
-       CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
-
-       CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
-       CALL fft_x( ar, 'backward' )
-       CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
-
-       CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
-       CALL transpose_xz( ar, work, ar )
-       CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
-
-       !$acc end data
-
-#endif
-
        CALL cpu_log( log_point_s(3), 'poisfft', 'stop' )
 
@@ -361,235 +341,4 @@
 
 
-    SUBROUTINE tridia( ar )
-
-!------------------------------------------------------------------------------!
-! solves the linear system of equations:
-!
-! -(4 pi^2(i^2/(dx^2*nnx^2)+j^2/(dy^2*nny^2))+
-!   1/(dzu(k)*dzw(k))+1/(dzu(k-1)*dzw(k)))*p(i,j,k)+
-! 1/(dzu(k)*dzw(k))*p(i,j,k+1)+1/(dzu(k-1)*dzw(k))*p(i,j,k-1)=d(i,j,k)
-!
-! by using the Thomas algorithm
-!------------------------------------------------------------------------------!
-
-       USE arrays_3d
-
-       IMPLICIT NONE
-
-       INTEGER ::  i, j, k, nnyh
-
-       REAL, DIMENSION(nxl_z:nxr_z,0:nz-1)   ::  ar1
-       REAL, DIMENSION(5,nxl_z:nxr_z,0:nz-1) ::  tri
-
-       REAL    ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
-
-
-       nnyh = (ny+1) / 2
-
-!
-!--    Define constant elements of the tridiagonal matrix.
-!$OMP  PARALLEL PRIVATE ( k, i )
-!$OMP  DO
-       DO  k = 0, nz-1
-          DO  i = nxl_z, nxr_z
-             tri(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
-             tri(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
-          ENDDO
-       ENDDO
-!$OMP  END PARALLEL
-
-#if defined( __parallel )
-!
-!--    Repeat for all y-levels.
-!$OMP  PARALLEL FIRSTPRIVATE( tri ) PRIVATE ( ar1, j )
-!$OMP  DO
-       DO  j = nys_z, nyn_z
-          IF ( j <= nnyh )  THEN
-             CALL maketri( j )
-          ELSE
-             CALL maketri( ny+1-j )
-          ENDIF
-          CALL split
-          CALL substi( j )
-       ENDDO
-!$OMP  END PARALLEL
-#else
-!
-!--    First y-level.
-       CALL maketri( nys_z )
-       CALL split
-       CALL substi( 0 )
-
-!
-!--    Further y-levels.
-       DO  j = 1, nnyh - 1
-          CALL maketri( j )
-          CALL split
-          CALL substi( j )
-          CALL substi( ny+1-j )
-       ENDDO
-       CALL maketri( nnyh )
-       CALL split
-       CALL substi( nnyh+nys )
-#endif
-
-    CONTAINS
-
-       SUBROUTINE maketri( j )
-
-!------------------------------------------------------------------------------!
-! Computes the i- and j-dependent component of the matrix
-!------------------------------------------------------------------------------!
-
-          USE arrays_3d
-          USE constants
-          USE control_parameters
-          USE grid_variables
-
-          IMPLICIT NONE
-
-          INTEGER ::  i, j, k, nnxh
-          REAL    ::  a, c
-          REAL    ::  ll(nxl_z:nxr_z)
-
-
-          nnxh = ( nx + 1 ) / 2
-
-!
-!--       Provide the tridiagonal matrix for solution of the Poisson equation in
-!--       Fourier space. The coefficients are computed following the method of
-!--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
-!--       Siano's original version by discretizing the Poisson equation,
-!--       before it is Fourier-transformed
-#if defined( __parallel )
-          DO  i = nxl_z, nxr_z
-             IF ( i >= 0 .AND. i <= nnxh )  THEN
-                ll(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
-                                          REAL( nx+1 ) ) ) / ( dx * dx ) + &
-                        2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
-                                          REAL( ny+1 ) ) ) / ( dy * dy )
-             ELSE
-                ll(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
-                                          REAL( nx+1 ) ) ) / ( dx * dx ) + &
-                        2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
-                                          REAL( ny+1 ) ) ) / ( dy * dy )
-             ENDIF
-             DO  k = 0,nz-1
-                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 - ll(i)
-             ENDDO
-          ENDDO
-#else
-          DO  i = 0, nnxh
-             ll(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / REAL( nx+1 ) ) ) / &
-                           ( dx * dx ) + &
-                     2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / REAL( ny+1 ) ) ) / &
-                           ( dy * dy )
-             DO  k = 0, nz-1
-                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 - ll(i)
-                IF ( i >= 1 .and. i < nnxh )  THEN
-                   tri(1,nx+1-i,k) = tri(1,i,k)
-                ENDIF
-             ENDDO
-          ENDDO
-#endif
-          IF ( ibc_p_b == 1  .OR.  ibc_p_b == 2 )  THEN
-             DO  i = nxl_z, nxr_z
-                tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
-             ENDDO
-          ENDIF
-          IF ( ibc_p_t == 1 )  THEN
-             DO  i = nxl_z, nxr_z
-                tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
-             ENDDO
-          ENDIF
-
-       END SUBROUTINE maketri
-
-
-       SUBROUTINE substi( j )
-
-!------------------------------------------------------------------------------!
-! Substitution (Forward and Backward) (Thomas algorithm)
-!------------------------------------------------------------------------------!
-
-          USE control_parameters
-
-          IMPLICIT NONE
-
-          INTEGER ::  i, j, k
-
-!
-!--       Forward substitution.
-          DO  i = nxl_z, nxr_z
-             ar1(i,0) = ar(i,j,1)
-          ENDDO
-          DO  k = 1, nz - 1
-             DO  i = nxl_z, nxr_z
-                ar1(i,k) = ar(i,j,k+1) - tri(5,i,k) * ar1(i,k-1)
-             ENDDO
-          ENDDO
-
-!
-!--       Backward substitution
-!--       Note, the 1.0E-20 in the denominator is due to avoid divisions
-!--       by zero appearing if the pressure bc is set to neumann at the top of
-!--       the model domain.
-          DO  i = nxl_z, nxr_z
-             ar(i,j,nz) = ar1(i,nz-1) / ( tri(4,i,nz-1) + 1.0E-20 )
-          ENDDO
-          DO  k = nz-2, 0, -1
-             DO  i = nxl_z, nxr_z
-                ar(i,j,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,j,k+2) ) &
-                              / tri(4,i,k)
-             ENDDO
-          ENDDO
-
-!
-!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
-!--       The respective values of ar should be zero at all k-levels if
-!--       acceleration of horizontally averaged vertical velocity is zero.
-          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
-             IF ( j == 0  .AND.  nxl_z == 0 )  THEN
-                DO  k = 1, nz
-                   ar(nxl_z,j,k) = 0.0
-                ENDDO
-             ENDIF
-          ENDIF
-
-       END SUBROUTINE substi
-
-
-       SUBROUTINE split
-
-!------------------------------------------------------------------------------!
-! Splitting of the tridiagonal matrix (Thomas algorithm)
-!------------------------------------------------------------------------------!
-
-          IMPLICIT NONE
-
-          INTEGER ::  i, k
-
-!
-!--       Splitting.
-          DO  i = nxl_z, nxr_z
-             tri(4,i,0) = tri(1,i,0)
-          ENDDO
-          DO  k = 1, nz-1
-             DO  i = nxl_z, nxr_z
-                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)
-             ENDDO
-          ENDDO
-
-       END SUBROUTINE split
-
-    END SUBROUTINE tridia
-
-
-#if defined( __parallel )
     SUBROUTINE ffty_tr_yx( f_in, work, f_out )
 
@@ -706,4 +455,5 @@
 !
 !--    Transpose array
+#if defined( __parallel )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
@@ -712,4 +462,5 @@
                           comm1dx, ierr )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
+#endif
 
     END SUBROUTINE ffty_tr_yx
@@ -745,4 +496,5 @@
 !
 !--    Transpose array
+#if defined( __parallel )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
@@ -751,4 +503,5 @@
                           comm1dx, ierr )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
+#endif
 
 !
@@ -1061,4 +814,5 @@
 !
 !--    Transpose array
+#if defined( __parallel )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
@@ -1067,4 +821,5 @@
                           comm1dy, ierr )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
+#endif
 
     END SUBROUTINE fftx_tr_xy
@@ -1096,4 +851,5 @@
 !
 !--    Transpose array
+#if defined( __parallel )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
@@ -1102,4 +858,5 @@
                           comm1dy, ierr )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
+#endif
 
 !
@@ -1426,5 +1183,5 @@
              ENDDO
           ENDDO
-          IF ( ibc_p_b == 1  .OR.  ibc_p_b == 2 )  THEN
+          IF ( ibc_p_b == 1 )  THEN
              DO  i = 0, nx
                 tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
@@ -1530,5 +1287,249 @@
     END SUBROUTINE tridia_1dd
 
-#endif
-#endif
+
+    SUBROUTINE tridia( ar )
+
+!------------------------------------------------------------------------------!
+! solves the linear system of equations:
+!
+! -(4 pi^2(i^2/(dx^2*nnx^2)+j^2/(dy^2*nny^2))+
+!   1/(dzu(k)*dzw(k))+1/(dzu(k-1)*dzw(k)))*p(i,j,k)+
+! 1/(dzu(k)*dzw(k))*p(i,j,k+1)+1/(dzu(k-1)*dzw(k))*p(i,j,k-1)=d(i,j,k)
+!
+! by using the Thomas algorithm
+!------------------------------------------------------------------------------!
+
+       USE arrays_3d
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+
+       !$acc declare create( tri )
+       REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) ::  tri
+
+       REAL    ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
+
+
+       CALL split( tri )
+       CALL substi( ar, tri )
+
+    END SUBROUTINE tridia
+
+
+    SUBROUTINE maketri
+
+!------------------------------------------------------------------------------!
+! Computes the i- and j-dependent component of the matrix
+!------------------------------------------------------------------------------!
+
+          USE arrays_3d,  ONLY: tric
+          USE constants
+          USE control_parameters
+          USE grid_variables
+
+          IMPLICIT NONE
+
+          INTEGER ::  i, j, k, nnxh, nnyh
+
+          !$acc declare create( ll )
+          REAL    ::  ll(nxl_z:nxr_z,nys_z:nyn_z)
+
+
+          nnxh = ( nx + 1 ) / 2
+          nnyh = ( ny + 1 ) / 2
+
+!
+!--       Provide the constant coefficients of the tridiagonal matrix for solution
+!--       of the Poisson equation in Fourier space.
+!--       The coefficients are computed following the method of
+!--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
+!--       Siano's original version by discretizing the Poisson equation,
+!--       before it is Fourier-transformed.
+
+          !$acc kernels present( tric )
+          !$acc loop vector( 32 )
+          DO  j = nys_z, nyn_z
+             DO  i = nxl_z, nxr_z
+                IF ( j >= 0  .AND.  j <= nnyh )  THEN
+                   IF ( i >= 0  .AND.  i <= nnxh )  THEN
+                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
+                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
+                                2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
+                                                  REAL( ny+1 ) ) ) / ( dy * dy )
+                   ELSE
+                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
+                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
+                                2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
+                                                  REAL( ny+1 ) ) ) / ( dy * dy )
+                   ENDIF
+                ELSE
+                   IF ( i >= 0  .AND.  i <= nnxh )  THEN
+                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
+                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
+                                2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
+                                                  REAL( ny+1 ) ) ) / ( dy * dy )
+                   ELSE
+                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
+                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
+                                2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
+                                                  REAL( ny+1 ) ) ) / ( dy * dy )
+                   ENDIF
+                ENDIF
+             ENDDO
+          ENDDO
+
+          !$acc loop
+          DO  k = 0, nz-1
+             DO  j = nys_z, nyn_z
+                !$acc loop vector( 32 )
+                DO  i = nxl_z, nxr_z
+                   tric(i,j,k) = ddzuw(k,3) - ll(i,j)
+                ENDDO
+             ENDDO
+          ENDDO
+          !$acc end kernels
+
+          IF ( ibc_p_b == 1 )  THEN
+             !$acc kernels present( tric )
+             !$acc loop
+             DO  j = nys_z, nyn_z
+                DO  i = nxl_z, nxr_z
+                   tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
+                ENDDO
+             ENDDO
+             !$acc end kernels
+          ENDIF
+          IF ( ibc_p_t == 1 )  THEN
+             !$acc kernels present( tric )
+             !$acc loop
+             DO  j = nys_z, nyn_z
+                DO  i = nxl_z, nxr_z
+                   tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2)
+                ENDDO
+             ENDDO
+             !$acc end kernels
+          ENDIF
+
+    END SUBROUTINE maketri
+
+
+    SUBROUTINE substi( ar, tri )
+
+!------------------------------------------------------------------------------!
+! Substitution (Forward and Backward) (Thomas algorithm)
+!------------------------------------------------------------------------------!
+
+          USE control_parameters
+
+          IMPLICIT NONE
+
+          INTEGER ::  i, j, k
+
+          REAL    ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
+          REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) ::  tri
+
+          !$acc declare create( ar1 )
+          REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1)   ::  ar1
+
+!
+!--       Forward substitution
+          DO  k = 0, nz - 1
+             !$acc kernels present( ar, tri )
+             !$acc loop
+             DO  j = nys_z, nyn_z
+                DO  i = nxl_z, nxr_z
+
+                   IF ( k == 0 )  THEN
+                      ar1(i,j,k) = ar(i,j,k+1)
+                   ELSE
+                      ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
+                   ENDIF
+
+                ENDDO
+             ENDDO
+             !$acc end kernels
+          ENDDO
+
+!
+!--       Backward substitution
+!--       Note, the 1.0E-20 in the denominator is due to avoid divisions
+!--       by zero appearing if the pressure bc is set to neumann at the top of
+!--       the model domain.
+          DO  k = nz-1, 0, -1
+             !$acc kernels present( ar, tri )
+             !$acc loop
+             DO  j = nys_z, nyn_z
+                DO  i = nxl_z, nxr_z
+
+                   IF ( k == nz-1 )  THEN
+                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20 )
+                   ELSE
+                      ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
+                              / tri(i,j,k,1)
+                   ENDIF
+                ENDDO
+             ENDDO
+             !$acc end kernels
+          ENDDO
+
+!
+!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
+!--       The respective values of ar should be zero at all k-levels if
+!--       acceleration of horizontally averaged vertical velocity is zero.
+          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
+             IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
+                !$acc kernels loop present( ar )
+                DO  k = 1, nz
+                   ar(nxl_z,nys_z,k) = 0.0
+                ENDDO
+             ENDIF
+          ENDIF
+
+    END SUBROUTINE substi
+
+
+    SUBROUTINE split( tri )
+
+!------------------------------------------------------------------------------!
+! Splitting of the tridiagonal matrix (Thomas algorithm)
+!------------------------------------------------------------------------------!
+
+          USE arrays_3d,  ONLY: tric
+
+          IMPLICIT NONE
+
+          INTEGER ::  i, j, k
+
+          REAL, DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,2) ::  tri
+
+!
+!--       Splitting
+          !$acc kernels present( tri, tric )
+          !$acc loop
+          DO  j = nys_z, nyn_z
+             !$acc loop vector( 32 )
+             DO  i = nxl_z, nxr_z
+                tri(i,j,0,1) = tric(i,j,0)
+             ENDDO
+          ENDDO
+          !$acc end kernels
+
+          DO  k = 1, nz-1
+             !$acc kernels present( tri, tric )
+             !$acc loop
+             DO  j = nys_z, nyn_z
+                !$acc loop vector( 32 )
+                DO  i = nxl_z, nxr_z
+                   tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
+                   tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
+                ENDDO
+             ENDDO
+             !$acc end kernels
+          ENDDO
+
+    END SUBROUTINE split
+
+#endif
+
  END MODULE poisfft_mod
Index: palm/trunk/SOURCE/poisfft_hybrid.f90
===================================================================
--- palm/trunk/SOURCE/poisfft_hybrid.f90	(revision 1109)
+++ palm/trunk/SOURCE/poisfft_hybrid.f90	(revision 1111)
@@ -20,5 +20,6 @@
 ! Current revisions: 
 ! ----------------- 
-! 
+! poisfft_hybrid_ini is now called internally from poisfft_hybrid,
+! ibc_p_b = 2 removed
 !
 ! Former revisions: 
@@ -115,4 +116,6 @@
                      tasks_per_logical_node = -1    ! default no cluster
 
+    LOGICAL, SAVE ::  poisfft_initialized = .FALSE.
+
 
     PRIVATE
@@ -283,4 +286,6 @@
        ENDIF
 
+       poisfft_initialized = .TRUE.
+
     END SUBROUTINE poisfft_hybrid_ini
 
@@ -294,4 +299,6 @@
 
        REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar
+
+       IF ( .NOT. poisfft_initialized )  CALL poisfft_hybrid_ini
 
        IF ( host(1:3) == 'nec' )  THEN
@@ -961,5 +968,5 @@
              ENDDO
           ENDDO
-          IF ( ibc_p_b == 1  .OR.  ibc_p_b == 2 )  THEN
+          IF ( ibc_p_b == 1 )  THEN
              DO  i = 0,nx
                 tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
Index: palm/trunk/SOURCE/pres.f90
===================================================================
--- palm/trunk/SOURCE/pres.f90	(revision 1109)
+++ palm/trunk/SOURCE/pres.f90	(revision 1111)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! openACC statements added,
+! ibc_p_b = 2 removed
 !
 ! Former revisions:
@@ -304,6 +305,9 @@
     ELSE
        !$OMP PARALLEL DO SCHEDULE( STATIC )
+       !$acc kernels present( d )
+       !$acc loop
        DO  i = nxl, nxr
           DO  j = nys, nyn
+             !$acc loop vector(32)
              DO  k = nzb+1, nzt
                 d(k,j,i) = 0.0
@@ -311,4 +315,5 @@
           ENDDO
        ENDDO
+       !$acc end kernels
     ENDIF
 
@@ -328,20 +333,4 @@
           ENDDO
 !
-!--       Additional pressure boundary condition at the bottom boundary for 
-!--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
-!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
-!--       This condition must not be applied at the start of a run, because then
-!--       flow_statistics has not yet been called and thus sums = 0.
-          IF ( ibc_p_b == 2  .AND.  sums(nzb+1,4) /= 0.0 )  THEN
-             k = nzb_s_inner(j,i)
-             d(k+1,j,i) = d(k+1,j,i) + (                                     &
-                                         ( usws(j,i+1) - usws(j,i) ) * ddx   &
-                                       + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
-                                       - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
-                                         sums(k+1,4)                         &
-                                       ) * ddzw(k+1) * ddt_3d * d_weight_pres 
-          ENDIF
-
-!
 !--       Compute possible PE-sum of divergences for flow_statistics
           DO  k = nzb_s_inner(j,i)+1, nzt
@@ -357,52 +346,24 @@
     !$OMP END PARALLEL
 #else
-    IF ( ibc_p_b == 2 .AND. sums(nzb+1,4) /= 0.0 )  THEN
-       !$OMP PARALLEL PRIVATE (i,j,k)
-       !$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) ) * ddt_3d      &
-                                                                * d_weight_pres 
-             ENDDO
-          ENDDO
-!
-!--       Additional pressure boundary condition at the bottom boundary for
-!--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
-!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
-!--       This condition must not be applied at the start of a run, because then
-!--       flow_statistics has not yet been called and thus sums = 0.
-          DO  j = nys, nyn
-              k = nzb_s_inner(j,i)
-              d(k+1,j,i) = d(k+1,j,i) + (                        &
-                             ( usws(j,i+1) - usws(j,i) ) * ddx   &
-                           + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
-                           - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
-                             sums(k+1,4)                         &
-                                        ) * ddzw(k+1) * ddt_3d   &
-                                          * d_weight_pres 
-          ENDDO
-       ENDDO
-       !$OMP END PARALLEL
-
-    ELSE
-
-       !$OMP PARALLEL PRIVATE (i,j,k)
-       !$OMP DO SCHEDULE( STATIC )
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt
+
+    !$OMP PARALLEL PRIVATE (i,j,k)
+    !$OMP DO SCHEDULE( STATIC )
+    !$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) ) * ddt_3d      &
-                                                                * d_weight_pres 
-             ENDDO
-          ENDDO
-       ENDDO
-       !$OMP END PARALLEL
-
-    ENDIF
+                           ( 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
+          ENDDO
+       ENDDO
+    ENDDO
+    !$acc end kernels
+    !$OMP END PARALLEL
 
 !
@@ -439,4 +400,5 @@
 !--       Solver for 2d-decomposition
           CALL poisfft( d, tend )
+          !$acc update host( d )
        ELSEIF ( psolver == 'poisfft_hybrid' )  THEN 
 !
@@ -466,17 +428,4 @@
 !
 !--       Neumann (dp/dz = 0)
-          !$OMP PARALLEL DO
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
-             ENDDO
-          ENDDO
-
-       ELSEIF ( ibc_p_b == 2 )  THEN
-!
-!--       Neumann condition for inhomogeneous surfaces,
-!--       here currently still in the form of a zero gradient. Actually 
-!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0 would have to be used for
-!--       the computation (cf. above: computation of divergences).
           !$OMP PARALLEL DO
           DO  i = nxlg, nxrg
@@ -611,4 +560,5 @@
 !-- 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,4 +698,6 @@
     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/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 1109)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 1111)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! update directives for prognostic quantities removed
 !
 ! Former revisions:
@@ -1574,5 +1574,4 @@
 
     CALL cpu_log( log_point(5), 'u-equation', 'stop' )
-    !$acc update host( u_p )
 
 !
@@ -1637,5 +1636,4 @@
 
     CALL cpu_log( log_point(6), 'v-equation', 'stop' )
-    !$acc update host( v_p )
 
 !
@@ -1700,5 +1698,4 @@
 
     CALL cpu_log( log_point(7), 'w-equation', 'stop' )
-    !$acc update host( w_p )
 
 
@@ -1796,5 +1793,4 @@
 
        CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
-       !$acc update host( pt_p )
 
     ENDIF
@@ -2046,5 +2042,4 @@
 
        CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
-       !$acc update host( e_p )
 
     ENDIF
Index: palm/trunk/SOURCE/swap_timelevel.f90
===================================================================
--- palm/trunk/SOURCE/swap_timelevel.f90	(revision 1109)
+++ palm/trunk/SOURCE/swap_timelevel.f90	(revision 1111)
@@ -20,4 +20,5 @@
 ! Current revisions:
 ! -----------------
+! openACC directives added
 ! 
 ! Former revisions:
@@ -82,10 +83,14 @@
     CALL cpu_log( log_point(28), 'swap_timelevel (nop)', 'start' )
 
+    !$acc kernels present( pt, pt_p, u, u_p, v, v_p, w, w_p )
     u  = u_p
     v  = v_p
     w  = w_p
     pt = pt_p
+    !$acc end kernels
     IF ( .NOT. constant_diffusion )  THEN
+       !$acc kernels present( e, e_p )
        e = e_p
+       !$acc end kernels
     ENDIF
     IF ( ocean )  THEN
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 1109)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 1111)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! +internal timestep counter for cpu statistics added,
+! openACC directives updated
 !
 ! Former revisions:
@@ -238,4 +239,5 @@
 !--       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 )
           CALL exchange_horiz( u_p, nbgp )
           CALL exchange_horiz( v_p, nbgp )
@@ -272,4 +274,5 @@
 !
 !--       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
 
@@ -298,4 +301,5 @@
              time_disturb = time_disturb + dt_3d
              IF ( time_disturb >= dt_disturb )  THEN
+                !$acc update host( u, v )
                 IF ( hom(nzb+5,1,pr_palm,0) < disturbance_energy_limit )  THEN
                    CALL disturb_field( nzb_u_inner, tend, u )
@@ -310,4 +314,5 @@
                    dist_range = 0
                 ENDIF
+                !$acc update device( u, v )
                 time_disturb = time_disturb - dt_disturb
              ENDIF
@@ -321,9 +326,4 @@
              CALL pres
           ENDIF
-!
-!--       Update device memory for calculating diffusion quantities and for next
-!--       timestep
-          !$acc update device( e, pt, u, v, w )
-          !$acc update device( q )  if ( allocated( q ) )
 
 !
@@ -351,7 +351,4 @@
                 CALL prandtl_fluxes
                 CALL cpu_log( log_point(19), 'prandtl_fluxes', 'stop' )
-!
-!++             Statistics still require updates on host
-                !$acc update host( qs, qsws, rif, shf, ts )
              ENDIF
 
@@ -369,7 +366,4 @@
              ENDIF
              CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
-!
-!++          Statistics still require update of diffusivities on host
-             !$acc update host( kh, km )
 
           ENDIF
@@ -379,4 +373,5 @@
 !
 !--    Increase simulation time and output times
+       nr_timesteps_this_run      = nr_timesteps_this_run + 1
        current_timestep_number    = current_timestep_number + 1
        simulated_time             = simulated_time   + dt_3d
Index: palm/trunk/SOURCE/transpose.f90
===================================================================
--- palm/trunk/SOURCE/transpose.f90	(revision 1109)
+++ palm/trunk/SOURCE/transpose.f90	(revision 1111)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! openACC directives added,
+! resorting data from/to work changed, work got 4 dimensions instead of 1
 !
 ! Former revisions:
@@ -81,10 +82,13 @@
     IMPLICIT NONE
 
-    INTEGER ::  i, j, k, l, m, ys
+    INTEGER ::  i, j, k, l, ys
     
-    REAL ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x),   &
-             f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx),  &
-             f_out(0:ny,nxl_y:nxr_y,nzb_y:nzt_y),  &
-             work(nnx*nny*nnz)
+    REAL ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x), f_out(0:ny,nxl_y:nxr_y,nzb_y:nzt_y)
+
+    REAL, DIMENSION(nyn_x-nys_x+1,nzb_y:nzt_y,nxl_y:nxr_y,0:pdims(2)-1) ::  work
+
+    !$acc declare create( f_inv )
+    REAL ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx)
+
 
 !
@@ -93,6 +97,9 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+    !$acc kernels present( f_in )
+    !$acc loop
     DO  i = 0, nx
        DO  k = nzb_x, nzt_x
+          !$acc loop vector( 32 )
           DO  j = nys_x, nyn_x
              f_inv(j,k,i) = f_in(i,j,k)
@@ -100,4 +107,5 @@
        ENDDO
     ENDDO
+    !$acc end kernels
 !$OMP  END PARALLEL
 
@@ -109,25 +117,28 @@
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLTOALL( f_inv(nys_x,nzb_x,0), sendrecvcount_xy, MPI_REAL, &
-                          work(1),              sendrecvcount_xy, MPI_REAL, &
+       !$acc update host( f_inv )
+       CALL MPI_ALLTOALL( f_inv(nys_x,nzb_x,0),  sendrecvcount_xy, MPI_REAL, &
+                          work(1,nzb_y,nxl_y,0), sendrecvcount_xy, MPI_REAL, &
                           comm1dy, ierr )
+       !$acc update device( work )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
 
 !
 !--    Reorder transposed array
-!$OMP  PARALLEL PRIVATE ( i, j, k, l, m, ys )
+!$OMP  PARALLEL PRIVATE ( i, j, k, l, ys )
 !$OMP  DO
        DO  l = 0, pdims(2) - 1
-          m  = l * ( nxr_y - nxl_y + 1 ) * ( nzt_y - nzb_y + 1 ) * &
-                   ( nyn_x - nys_x + 1 )
           ys = 0 + l * ( nyn_x - nys_x + 1 )
+          !$acc kernels present( f_out, work )
+          !$acc loop
           DO  i = nxl_y, nxr_y
              DO  k = nzb_y, nzt_y
+                !$acc loop vector( 32 )
                 DO  j = ys, ys + nyn_x - nys_x
-                   m = m + 1
-                   f_out(j,i,k) = work(m)
+                   f_out(j,i,k) = work(j-ys+1,k,i,l)
                 ENDDO
              ENDDO
           ENDDO
+          !$acc end kernels
        ENDDO
 !$OMP  END PARALLEL
@@ -140,6 +151,9 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+       !$acc kernels present( f_out )
+       !$acc loop
        DO  k = nzb_y, nzt_y
           DO  i = nxl_y, nxr_y
+             !$acc loop vector( 32 )
              DO  j = 0, ny
                 f_out(j,i,k) = f_inv(j,k,i)
@@ -147,4 +161,5 @@
           ENDDO
        ENDDO
+       !$acc end kernels
 !$OMP  END PARALLEL
 
@@ -172,10 +187,12 @@
     IMPLICIT NONE
 
-    INTEGER ::  i, j, k, l, m, xs
+    INTEGER ::  i, j, k, l, xs
     
-    REAL ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x),  &
-             f_inv(nys:nyn,nxl:nxr,1:nz),         &
-             f_out(1:nz,nys:nyn,nxl:nxr),         &
-             work(nnx*nny*nnz)
+    REAL ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x), f_out(1:nz,nys:nyn,nxl:nxr)
+
+    REAL, DIMENSION(nys_x:nyn_x,nnx,nzb_x:nzt_x,0:pdims(1)-1) ::  work
+
+    !$acc declare create( f_inv )
+    REAL ::  f_inv(nys:nyn,nxl:nxr,1:nz)
 
 
@@ -188,17 +205,19 @@
 !
 !--    Reorder input array for transposition
-!$OMP  PARALLEL PRIVATE ( i, j, k, l, m, xs )
+!$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
 !$OMP  DO
        DO  l = 0, pdims(1) - 1
-          m  = l * ( nzt_x - nzb_x + 1 ) * nnx * ( nyn_x - nys_x + 1 )
           xs = 0 + l * nnx
+          !$acc kernels present( f_in, work )
+          !$acc loop
           DO  k = nzb_x, nzt_x
              DO  i = xs, xs + nnx - 1
+                !$acc loop vector( 32 )
                 DO  j = nys_x, nyn_x
-                   m = m + 1
-                   work(m) = f_in(i,j,k)
+                   work(j,i-xs+1,k,l) = f_in(i,j,k)
                 ENDDO
              ENDDO
           ENDDO
+          !$acc end kernels
        ENDDO
 !$OMP  END PARALLEL
@@ -208,7 +227,9 @@
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLTOALL( work(1),          sendrecvcount_zx, MPI_REAL, &
-                          f_inv(nys,nxl,1), sendrecvcount_zx, MPI_REAL, &
+       !$acc update host( work )
+       CALL MPI_ALLTOALL( work(nys_x,1,nzb_x,0), sendrecvcount_zx, MPI_REAL, &
+                          f_inv(nys,nxl,1),      sendrecvcount_zx, MPI_REAL, &
                           comm1dx, ierr )
+       !$acc update device( f_inv )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
 
@@ -217,6 +238,9 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+       !$acc kernels present( f_out )
+       !$acc loop
        DO  k = 1, nz
           DO  i = nxl, nxr
+             !$acc loop vector( 32 )
              DO  j = nys, nyn
                 f_out(k,j,i) = f_inv(j,i,k)
@@ -224,4 +248,5 @@
           ENDDO
        ENDDO
+       !$acc end kernels
 !$OMP  END PARALLEL
 #endif
@@ -233,6 +258,9 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+       !$acc kernels present( f_in )
+       !$acc loop
        DO  i = nxl, nxr
           DO  j = nys, nyn
+             !$acc loop vector( 32 )
              DO  k = 1, nz
                 f_inv(j,i,k) = f_in(i,j,k)
@@ -240,10 +268,14 @@
           ENDDO
        ENDDO
-!$OMP  END PARALLEL
-
-!$OMP  PARALLEL PRIVATE ( i, j, k )
-!$OMP  DO
+       !$acc end kernels
+!$OMP  END PARALLEL
+
+!$OMP  PARALLEL PRIVATE ( i, j, k )
+!$OMP  DO
+       !$acc kernels present( f_out )
+       !$acc loop
        DO  k = 1, nz
           DO  i = nxl, nxr
+             !$acc loop vector( 32 )
              DO  j = nys, nyn
                 f_out(k,j,i) = f_inv(j,i,k)
@@ -251,4 +283,5 @@
           ENDDO
        ENDDO
+       !$acc end kernels
 !$OMP  END PARALLEL
 
@@ -276,10 +309,13 @@
     IMPLICIT NONE
 
-    INTEGER ::  i, j, k, l, m, ys
+    INTEGER ::  i, j, k, l, ys
     
-    REAL ::  f_in(0:ny,nxl_y:nxr_y,nzb_y:nzt_y),  &
-             f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx), &
-             f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x), &
-             work(nnx*nny*nnz)
+    REAL ::  f_in(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x)
+
+    REAL, DIMENSION(nyn_x-nys_x+1,nzb_y:nzt_y,nxl_y:nxr_y,0:pdims(2)-1) ::  work
+
+    !$acc declare create( f_inv )
+    REAL ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx)
+
 
     IF ( numprocs /= 1 )  THEN
@@ -288,18 +324,19 @@
 !
 !--    Reorder input array for transposition
-!$OMP  PARALLEL PRIVATE ( i, j, k, l, m, ys )
+!$OMP  PARALLEL PRIVATE ( i, j, k, l, ys )
 !$OMP  DO
        DO  l = 0, pdims(2) - 1
-          m  = l * ( nxr_y - nxl_y + 1 ) * ( nzt_y - nzb_y + 1 ) * &
-                   ( nyn_x - nys_x + 1 )
           ys = 0 + l * ( nyn_x - nys_x + 1 )
+          !$acc kernels present( f_in, work )
+          !$acc loop
           DO  i = nxl_y, nxr_y
              DO  k = nzb_y, nzt_y
+                !$acc loop vector( 32 )
                 DO  j = ys, ys + nyn_x - nys_x
-                   m = m + 1
-                   work(m) = f_in(j,i,k)
+                   work(j-ys+1,k,i,l) = f_in(j,i,k)
                 ENDDO
              ENDDO
           ENDDO
+          !$acc end kernels
        ENDDO
 !$OMP  END PARALLEL
@@ -309,7 +346,9 @@
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLTOALL( work(1),              sendrecvcount_xy, MPI_REAL, &
-                          f_inv(nys_x,nzb_x,0), sendrecvcount_xy, MPI_REAL, &
+       !$acc update host( work )
+       CALL MPI_ALLTOALL( work(1,nzb_y,nxl_y,0), sendrecvcount_xy, MPI_REAL, &
+                          f_inv(nys_x,nzb_x,0),  sendrecvcount_xy, MPI_REAL, &
                           comm1dy, ierr )
+       !$acc update device( f_inv )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
 #endif
@@ -321,6 +360,9 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+       !$acc kernels present( f_in )
+       !$acc loop
        DO  i = nxl_y, nxr_y
           DO  k = nzb_y, nzt_y
+             !$acc loop vector( 32 )
              DO  j = 0, ny
                 f_inv(j,k,i) = f_in(j,i,k)
@@ -328,4 +370,5 @@
           ENDDO
        ENDDO
+       !$acc end kernels
 !$OMP  END PARALLEL
 
@@ -336,6 +379,9 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+    !$acc kernels present( f_out )
+    !$acc loop
     DO  i = 0, nx
        DO  k = nzb_x, nzt_x
+          !$acc loop vector( 32 )
           DO  j = nys_x, nyn_x
              f_out(i,j,k) = f_inv(j,k,i)
@@ -343,4 +389,5 @@
        ENDDO
     ENDDO
+    !$acc end kernels
 !$OMP  END PARALLEL
 
@@ -434,10 +481,13 @@
     IMPLICIT NONE
 
-    INTEGER ::  i, j, k, l, m, zs
+    INTEGER ::  i, j, k, l, zs
     
-    REAL ::  f_in(0:ny,nxl_y:nxr_y,nzb_y:nzt_y),  &
-             f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny), &
-             f_out(nxl_z:nxr_z,nys_z:nyn_z,1:nz), &
-             work(nnx*nny*nnz)
+    REAL ::  f_in(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), f_out(nxl_z:nxr_z,nys_z:nyn_z,1:nz)
+
+    REAL, DIMENSION(nxl_z:nxr_z,nzt_y-nzb_y+1,nys_z:nyn_z,0:pdims(1)-1) ::  work
+
+    !$acc declare create( f_inv )
+    REAL ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny)
+
 
 !
@@ -446,6 +496,9 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+    !$acc kernels present( f_in )
+    !$acc loop
     DO  j = 0, ny
        DO  k = nzb_y, nzt_y
+          !$acc loop vector( 32 )
           DO  i = nxl_y, nxr_y
              f_inv(i,k,j) = f_in(j,i,k)
@@ -453,4 +506,5 @@
        ENDDO
     ENDDO
+    !$acc end kernels
 !$OMP  END PARALLEL
 
@@ -464,6 +518,9 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+       !$acc kernels present( f_out )
+       !$acc loop
        DO  j = 0, ny
           DO  k = nzb_y, nzt_y
+             !$acc loop vector( 32 )
              DO  i = nxl_y, nxr_y
                 f_out(i,j,k) = f_inv(i,k,j)
@@ -471,4 +528,5 @@
           ENDDO
        ENDDO
+       !$acc end kernels
 !$OMP  END PARALLEL
 
@@ -480,25 +538,28 @@
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLTOALL( f_inv(nxl_y,nzb_y,0), sendrecvcount_yz, MPI_REAL, &
-                          work(1),              sendrecvcount_yz, MPI_REAL, &
+       !$acc update host( f_inv )
+       CALL MPI_ALLTOALL( f_inv(nxl_y,nzb_y,0),  sendrecvcount_yz, MPI_REAL, &
+                          work(nxl_z,1,nys_z,0), sendrecvcount_yz, MPI_REAL, &
                           comm1dx, ierr )
+       !$acc update device( work )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
 
 !
 !--    Reorder transposed array
-!$OMP  PARALLEL PRIVATE ( i, j, k, l, m, zs )
+!$OMP  PARALLEL PRIVATE ( i, j, k, l, zs )
 !$OMP  DO
        DO  l = 0, pdims(1) - 1
-          m  = l * ( nyn_z - nys_z + 1 ) * ( nzt_y - nzb_y + 1 ) * &
-                   ( nxr_z - nxl_z + 1 )
           zs = 1 + l * ( nzt_y - nzb_y + 1 )
+          !$acc kernels present( f_out, work )
+          !$acc loop
           DO  j = nys_z, nyn_z
              DO  k = zs, zs + nzt_y - nzb_y
+                !$acc loop vector( 32 )
                 DO  i = nxl_z, nxr_z
-                   m = m + 1
-                   f_out(i,j,k) = work(m)
+                   f_out(i,j,k) = work(i,k-zs+1,j,l)
                 ENDDO
              ENDDO
           ENDDO
+          !$acc end kernels
        ENDDO
 !$OMP  END PARALLEL
@@ -528,10 +589,11 @@
     IMPLICIT NONE
 
-    INTEGER ::  i, j, k, l, m, xs
+    INTEGER ::  i, j, k, l, xs
     
-    REAL ::  f_in(1:nz,nys:nyn,nxl:nxr), f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x), &
-             work(nnx*nny*nnz)
-
-    !$acc declare create ( f_inv )
+    REAL ::  f_in(1:nz,nys:nyn,nxl:nxr), f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x)
+
+    REAL, DIMENSION(nys_x:nyn_x,nnx,nzb_x:nzt_x,0:pdims(1)-1) ::  work
+
+    !$acc declare create( f_inv )
     REAL ::  f_inv(nys:nyn,nxl:nxr,1:nz)
 
@@ -552,4 +614,5 @@
        ENDDO
     ENDDO
+    !$acc end kernels
 !$OMP  END PARALLEL
 
@@ -573,4 +636,5 @@
           ENDDO
        ENDDO
+       !$acc end kernels
 !$OMP  END PARALLEL
 
@@ -582,24 +646,28 @@
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLTOALL( f_inv(nys,nxl,1), sendrecvcount_zx, MPI_REAL, &
-                          work(1),          sendrecvcount_zx, MPI_REAL, &
+       !$acc update host( f_inv )
+       CALL MPI_ALLTOALL( f_inv(nys,nxl,1),      sendrecvcount_zx, MPI_REAL, &
+                          work(nys_x,1,nzb_x,0), sendrecvcount_zx, MPI_REAL, &
                           comm1dx, ierr )
+       !$acc update device( work )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
 
 !
 !--    Reorder transposed array
-!$OMP  PARALLEL PRIVATE ( i, j, k, l, m, xs )
+!$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
 !$OMP  DO
        DO  l = 0, pdims(1) - 1
-          m  = l * ( nzt_x - nzb_x + 1 ) * nnx * ( nyn_x - nys_x + 1 )
           xs = 0 + l * nnx
+          !$acc kernels present( f_out, work )
+          !$acc loop
           DO  k = nzb_x, nzt_x
              DO  i = xs, xs + nnx - 1
+                !$acc loop vector( 32 )
                 DO  j = nys_x, nyn_x
-                   m = m + 1
-                   f_out(i,j,k) = work(m)
+                   f_out(i,j,k) = work(j,i-xs+1,k,l)
                 ENDDO
              ENDDO
           ENDDO
+          !$acc end kernels
        ENDDO
 !$OMP  END PARALLEL
@@ -629,10 +697,13 @@
     IMPLICIT NONE
 
-    INTEGER ::  i, j, k, l, m, zs
+    INTEGER ::  i, j, k, l, zs
     
-    REAL ::  f_in(nxl_z:nxr_z,nys_z:nyn_z,1:nz),  &
-             f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny), &
-             f_out(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), &
-             work(nnx*nny*nnz)
+    REAL ::  f_in(nxl_z:nxr_z,nys_z:nyn_z,1:nz), f_out(0:ny,nxl_y:nxr_y,nzb_y:nzt_y)
+
+    REAL, DIMENSION(nxl_z:nxr_z,nzt_y-nzb_y+1,nys_z:nyn_z,0:pdims(1)-1) ::  work
+
+    !$acc declare create( f_inv )
+    REAL ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny)
+
 
 !
@@ -644,18 +715,19 @@
 !
 !--    Reorder input array for transposition
-!$OMP  PARALLEL PRIVATE ( i, j, k, l, m, zs )
+!$OMP  PARALLEL PRIVATE ( i, j, k, l, zs )
 !$OMP  DO
        DO  l = 0, pdims(1) - 1
-          m  = l * ( nyn_z - nys_z + 1 ) * ( nzt_y - nzb_y + 1 ) * &
-                   ( nxr_z - nxl_z + 1 )
           zs = 1 + l * ( nzt_y - nzb_y + 1 )
+          !$acc kernels present( f_in, work )
+          !$acc loop
           DO  j = nys_z, nyn_z
              DO  k = zs, zs + nzt_y - nzb_y
+                !$acc loop vector( 32 )
                 DO  i = nxl_z, nxr_z
-                   m = m + 1
-                   work(m) = f_in(i,j,k)
+                   work(i,k-zs+1,j,l) = f_in(i,j,k)
                 ENDDO
              ENDDO
           ENDDO
+          !$acc end kernels
        ENDDO
 !$OMP  END PARALLEL
@@ -665,7 +737,9 @@
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-       CALL MPI_ALLTOALL( work(1),              sendrecvcount_yz, MPI_REAL, &
-                          f_inv(nxl_y,nzb_y,0), sendrecvcount_yz, MPI_REAL, &
+       !$acc update host( work )
+       CALL MPI_ALLTOALL( work(nxl_z,1,nys_z,0), sendrecvcount_yz, MPI_REAL, &
+                          f_inv(nxl_y,nzb_y,0),  sendrecvcount_yz, MPI_REAL, &
                           comm1dx, ierr )
+       !$acc update device( f_inv )
        CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
 #endif
@@ -676,6 +750,9 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+       !$acc kernels present( f_in )
+       !$acc loop
        DO  k = nzb_y, nzt_y
           DO  j = 0, ny
+             !$acc loop vector( 32 )
              DO  i = nxl_y, nxr_y
                 f_inv(i,k,j) = f_in(i,j,k)
@@ -683,4 +760,5 @@
           ENDDO
        ENDDO
+       !$acc end kernels
 !$OMP  END PARALLEL
 
@@ -691,6 +769,9 @@
 !$OMP  PARALLEL PRIVATE ( i, j, k )
 !$OMP  DO
+    !$acc kernels present( f_out )
+    !$acc loop
     DO  k = nzb_y, nzt_y
        DO  i = nxl_y, nxr_y
+          !$acc loop vector( 32 )
           DO  j = 0, ny
              f_out(j,i,k) = f_inv(i,k,j)
@@ -698,4 +779,5 @@
        ENDDO
     ENDDO
+    !$acc end kernels
 !$OMP  END PARALLEL
 
