       dev_maxMinArray.getPtr(),
       dev_a1.getPtr(), dev_a2.getPtr(), 
       dev_a3.getPtr(), dev_a4.getPtr(),
-      NY * NX1, numParticles, particlesToBuffer
+      NY * NX1, numRemainingParticles, particlesToBuffer
       );
 
    checkCuda(cudaStreamSynchronize(stream1));

      );
   checkForCudaError("calcA failed");

   if(sortParticleArray)
   {
      assert(numParticles < 1000000);
      //std::cout << "Sorting!" << std::endl;
      threadsInBlock = 256;
      blockSize = dim3(threadsInBlock);
      numBlocks = dim3(static_cast<unsigned int>(calcNumBlocks(threadsInBlock, numParticles)));
      sharedMemoryBytes = d_partLoc.getElementSize() * (threadsInBlock + 1);
      DevMem<unsigned int, DevMemReuse> d_firstOob(1);
      d_firstOob.copyArrayToDev(&numParticles, 1);
      //std::cout << "numParticles=" << numParticles << " OOB_PARTICLE=" << OOB_PARTICLE << std::endl;
      //std::cout << "blockSize.x=" << blockSize.x << " numBlocks=" << numBlocks.x << " sharedMem=" << sharedMemoryBytes << std::endl;
      //std::cout << "dev_particleLocations.size()=" << dev_particleLocations.size() << std::endl;
      findFirstOob<<<numBlocks, blockSize, sharedMemoryBytes, stream2>>>(
         d_partLoc.getPtr(),
         numParticles,
         OOB_PARTICLE,
         d_firstOob.getPtr()
         );
      checkForCudaError("findFirstOob failed");
      checkCuda(cudaStreamSynchronize(stream2));
      d_firstOob.copyArrayToHost(&numParticles, 1);
      assert(numParticles < 1000000);
   }



      dev_a3Sum.getPtr(),
      dev_a4Sum.getPtr(),
      cold,
      numRemainingParticles,
      NIJ
      );
   topRowBlockBoundaries.freeMem();
