 
 thrust::device_vector<unsigned int> *keys;       //!< Temporary vector of sort keys
 
+unsigned int *d_n_copy_ghosts;     //! Counter for ghost list construction
+
 void gpu_allocate_tmp_storage()
     {
     keys = new thrust::device_vector<unsigned int>;
+
+    cudaMalloc(&d_n_copy_ghosts, sizeof(unsigned int));
     }
 
 void gpu_deallocate_tmp_storage()
     {
     delete keys;
+
+    cudaFree(d_n_copy_ghosts);
     }
 
 //! GPU Kernel to find incomplete bonds

        return plan;
        }
     };

//! Select ghost particles for sending in one direction
struct select_particle_ghost
    {
    const unsigned int dir; //!< Current direction
    const unsigned int N;   //!< Number of local particles

    //! Constructor
    /*! \param _dir Direction of the neighboring domain
     */
    select_particle_ghost(unsigned int _dir, unsigned int _N)
        : dir(_dir), N(_N)
        {
        }

    //! Select particles for sending
    /*! \param plan Particle exchange plan
        \returns true if particle is selected for sending
     */
    __host__ __device__ bool operator()(const thrust::tuple<unsigned int, unsigned char>& t)
        {
        unsigned int idx = thrust::get<0>(t);
        unsigned char plan = thrust::get<1>(t);

        bool result = (plan & (1 << dir));

        // do not send back ghost particles just received from the opposite direction
        if ((dir%2) && (plan & (1 << (dir-1))) && idx >= N)
            result = false;

        return result;
        }
     };


//! Structure to pack a particle data element into
struct __align__(128) pdata_element_gpu

                                  unsigned int dir,
                                  unsigned char *d_plan,
                                  unsigned int *d_global_tag,
                                  unsigned int *d_copy_ghosts,
                                  unsigned int &n_copy_ghosts)
    {
    n_copy_ghosts = 0;
    cudaMemcpy(d_n_copy_ghosts, &n_copy_ghosts, sizeof(unsigned int), cudaMemcpyHostToDevice);

    unsigned int block_size = 512;
    gpu_make_exchange_ghost_list_kernel<<<n_total/block_size + 1, block_size>>>(n_total,
                                                                                N,
                                                                                dir,
                                                                                d_plan,
                                                                                d_global_tag,
                                                                                d_copy_ghosts,
                                                                                d_n_copy_ghosts);
    cudaMemcpy(&n_copy_ghosts, d_n_copy_ghosts, sizeof(unsigned int), cudaMemcpyDeviceToHost);
    }

//! Fill send buffers of particles we are sending as ghost particles with partial particle data
