diff --git a/CMakeLists.txt b/CMakeLists.txt
index 000798d..96bc953 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -171,11 +171,12 @@ endif()
 # Note that these options may not play nice with nvcc wrapper
 if (TEST_INTEL_OPTIMIZATION)
   add_compile_options(-fp-model fast=2 -qopt_report5 -vec-threshold0 -qopt_report_phase=vec)
-endif()	       
+endif()
 
 add_subdirectory(src)
 add_subdirectory(example/calculate_pi)
 add_subdirectory(example/face_fields)
 add_subdirectory(example/advection)
+add_subdirectory(example/particles)
 
 include(cmake/CheckCopyright.cmake)
diff --git a/example/particles/CMakeLists.txt b/example/particles/CMakeLists.txt
new file mode 100644
index 0000000..7a87429
--- /dev/null
+++ b/example/particles/CMakeLists.txt
@@ -0,0 +1,19 @@
+#=========================================================================================
+# (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+#
+# This program was produced under U.S. Government contract 89233218CNA000001 for Los
+# Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+# for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+# in the program are reserved by Triad National Security, LLC, and the U.S. Department
+# of Energy/National Nuclear Security Administration. The Government is granted for
+# itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+# license in this material to reproduce, prepare derivative works, distribute copies to
+# the public, perform publicly and display publicly, and to permit others to do so.
+#=========================================================================================
+
+add_executable(
+    particles-example
+        main.cpp
+        particles.cpp
+)
+target_link_libraries(particles-example PRIVATE parthenon)
diff --git a/example/particles/main.cpp b/example/particles/main.cpp
new file mode 100644
index 0000000..fb11582
--- /dev/null
+++ b/example/particles/main.cpp
@@ -0,0 +1,60 @@
+//========================================================================================
+// (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+#include "parthenon_manager.hpp"
+
+#include "particles.hpp"
+
+int main(int argc, char *argv[]) {
+  using parthenon::ParthenonManager;
+  using parthenon::ParthenonStatus;
+  ParthenonManager pman;
+
+  // call ParthenonInit to initialize MPI and Kokkos, parse the input deck, and set up
+  auto manager_status = pman.ParthenonInit(argc, argv);
+  if (manager_status == ParthenonStatus::complete) {
+    pman.ParthenonFinalize();
+    return 0;
+  }
+  if (manager_status == ParthenonStatus::error) {
+    pman.ParthenonFinalize();
+    return 1;
+  }
+  printf("File: %s Line: %i\n", __FILE__, __LINE__);
+  // Now that ParthenonInit has been called and setup succeeded, the code can now
+  // make use of MPI and Kokkos
+
+  // Initialize the driver
+  particles_example::ParticleDriver driver(pman.pinput.get(), pman.pmesh.get(),
+                                            pman.pouts.get());
+  printf("File: %s Line: %i\n", __FILE__, __LINE__);
+
+  // start a timer
+  pman.PreDriver();
+  printf("File: %s Line: %i\n", __FILE__, __LINE__);
+
+  // This line actually runs the simulation
+  auto driver_status = driver.Execute();
+  printf("File: %s Line: %i\n", __FILE__, __LINE__);
+
+  // Make final outputs, print diagnostics
+  pman.PostDriver(driver_status);
+  printf("File: %s Line: %i\n", __FILE__, __LINE__);
+
+  // call MPI_Finalize and Kokkos::finalize if necessary
+  pman.ParthenonFinalize();
+
+  // MPI and Kokkos can no longer be used
+
+  return (0);
+}
diff --git a/example/particles/parthinput.particles b/example/particles/parthinput.particles
new file mode 100644
index 0000000..fca535b
--- /dev/null
+++ b/example/particles/parthinput.particles
@@ -0,0 +1,60 @@
+# ========================================================================================
+#  Athena++ astrophysical MHD code
+#  Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
+#  Licensed under the 3-clause BSD License, see LICENSE file for details
+# ========================================================================================
+#  (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+#
+#  This program was produced under U.S. Government contract 89233218CNA000001 for Los
+#  Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+#  for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+#  in the program are reserved by Triad National Security, LLC, and the U.S. Department
+#  of Energy/National Nuclear Security Administration. The Government is granted for
+#  itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+#  license in this material to reproduce, prepare derivative works, distribute copies to
+#  the public, perform publicly and display publicly, and to permit others to do so.
+# ========================================================================================
+
+<comment>
+problem = Particle streaming
+
+<job>
+problem_id = particles
+
+<mesh>
+refinement = none
+#numlevel = 3
+
+nx1 = 64
+x1min = -0.5
+x1max = 0.5
+ix1_bc = periodic
+ox1_bc = periodic
+
+nx2 = 64
+x2min = -0.5
+x2max = 0.5
+ix2_bc = periodic
+ox2_bc = periodic
+
+nx3 = 1
+x3min = -0.5
+x3max = 0.5
+
+<meshblock>
+nx1 = 16
+nx2 = 16
+#nx3 = 1
+
+<output1>
+file_type = hdf5
+dt = 0.05
+variable = prim
+
+<time>
+tlim = 1.0
+integrator = rk2
+
+<Particles>
+num_particles = 100
+particle_speed = 1.0
diff --git a/example/particles/particles.cpp b/example/particles/particles.cpp
new file mode 100644
index 0000000..1880a7b
--- /dev/null
+++ b/example/particles/particles.cpp
@@ -0,0 +1,450 @@
+//========================================================================================
+// (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+#include "particles.hpp"
+
+#include <algorithm>
+#include <limits>
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "bvals/boundary_conditions.hpp"
+#include "bvals/bvals.hpp"
+#include "driver/multistage.hpp"
+#include "interface/params.hpp"
+#include "interface/state_descriptor.hpp"
+#include "mesh/mesh.hpp"
+#include "parthenon_manager.hpp"
+#include "reconstruct/reconstruction.hpp"
+#include "refinement/refinement.hpp"
+
+using parthenon::BlockStageNamesIntegratorTask;
+using parthenon::BlockStageNamesIntegratorTaskFunc;
+using parthenon::BlockTask;
+using parthenon::CellVariable;
+using parthenon::Integrator;
+using parthenon::Metadata;
+using parthenon::Params;
+using parthenon::ParArrayND;
+using parthenon::ParthenonManager;
+
+// *************************************************//
+// redefine some weakly linked parthenon functions *//
+// *************************************************//
+
+namespace parthenon {
+
+Packages_t ParthenonManager::ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
+  Packages_t packages;
+  packages["Particles"] = particles_example::Particles::Initialize(pin.get());
+  return packages;
+}
+
+void MeshBlock::ProblemGenerator(ParameterInput *pin) {
+  SwarmContainer &sc = real_containers.GetSwarm();
+  printf("%s %i\n", __FILE__, __LINE__);
+  Swarm &s = sc.Get("particles");
+  //CellVariable<Real> &q = rc.Get("advected");
+
+  s.AddParticle();
+
+
+
+  /*for (int k = 0; k < ncells3; k++) {
+    for (int j = 0; j < ncells2; j++) {
+      for (int i = 0; i < ncells1; i++) {
+        Real rsq = std::pow(pcoord->x1v(i), 2) + std::pow(pcoord->x2v(j), 2);
+        q(k, j, i) = (rsq < 0.15 * 0.15 ? 1.0 : 0.0);
+      }
+    }
+  }*/
+}
+
+void ParthenonManager::SetFillDerivedFunctions() {
+  FillDerivedVariables::SetFillDerivedFunctions(particles_example::Particles::PreFill,
+                                                particles_example::Particles::PostFill);
+}
+
+} // namespace parthenon
+
+// *************************************************//
+// define the "physics" package Particles, which   *//
+// includes defining various functions that control*//
+// how parthenon functions and any tasks needed to *//
+// implement the "physics"                         *//
+// *************************************************//
+
+namespace particles_example {
+namespace Particles {
+
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
+  auto pkg = std::make_shared<StateDescriptor>("Particles");
+
+  int num_particles = pin->GetOrAddInteger("Particles", "num_particles", 100);
+  pkg->AddParam<>("num_particles", num_particles);
+  Real particle_speed = pin->GetOrAddReal("Particles", "particle_speed", 1.0);
+  pkg->AddParam<>("particle_speed", particle_speed);
+
+  std::string swarm_name = "particles";
+  Metadata m;
+  pkg->AddSwarm(swarm_name, m);
+
+  //pkg->FillDerived = SquareIt;
+  //pkg->CheckRefinement = CheckRefinement;
+  //pkg->EstimateTimestep = EstimateTimestep;
+
+  return pkg;
+}
+
+AmrTag CheckRefinement(Container<Real> &rc) {
+  /*MeshBlock *pmb = rc.pmy_block;
+  // refine on advected, for example.  could also be a derived quantity
+  CellVariable<Real> &v = rc.Get("advected");
+  Real vmin = 1.0;
+  Real vmax = 0.0;
+  for (int k = 0; k < pmb->ncells3; k++) {
+    for (int j = 0; j < pmb->ncells2; j++) {
+      for (int i = 0; i < pmb->ncells1; i++) {
+        vmin = (v(k, j, i) < vmin ? v(k, j, i) : vmin);
+        vmax = (v(k, j, i) > vmax ? v(k, j, i) : vmax);
+      }
+    }
+  }
+  auto pkg = pmb->packages["Advection"];
+  const auto &refine_tol = pkg->Param<Real>("refine_tol");
+  const auto &derefine_tol = pkg->Param<Real>("derefine_tol");
+
+  if (vmax > refine_tol && vmin < derefine_tol) return AmrTag::refine;
+  if (vmax < derefine_tol) return AmrTag::derefine;*/
+  return AmrTag::same;
+}
+
+// demonstrate usage of a "pre" fill derived routine
+void PreFill(Container<Real> &rc) {
+  /*MeshBlock *pmb = rc.pmy_block;
+  int is = 0;
+  int js = 0;
+  int ks = 0;
+  int ie = pmb->ncells1 - 1;
+  int je = pmb->ncells2 - 1;
+  int ke = pmb->ncells3 - 1;
+  CellVariable<Real> &qin = rc.Get("advected");
+  CellVariable<Real> &qout = rc.Get("one_minus_advected");
+  for (int i = is; i <= ie; i++) {
+    for (int j = js; j <= je; j++) {
+      for (int k = ks; k <= ke; k++) {
+        qout(k, j, i) = 1.0 - qin(k, j, i);
+      }
+    }
+  }*/
+}
+
+// this is the package registered function to fill derived
+void SquareIt(Container<Real> &rc) {
+  /*MeshBlock *pmb = rc.pmy_block;
+  int is = 0;
+  int js = 0;
+  int ks = 0;
+  int ie = pmb->ncells1 - 1;
+  int je = pmb->ncells2 - 1;
+  int ke = pmb->ncells3 - 1;
+  CellVariable<Real> &qin = rc.Get("one_minus_advected");
+  CellVariable<Real> &qout = rc.Get("one_minus_advected_sq");
+  for (int i = is; i <= ie; i++) {
+    for (int j = js; j <= je; j++) {
+      for (int k = ks; k <= ke; k++) {
+        qout(k, j, i) = qin(k, j, i) * qin(k, j, i);
+      }
+    }
+  }*/
+}
+
+// demonstrate usage of a "post" fill derived routine
+void PostFill(Container<Real> &rc) {
+  /*MeshBlock *pmb = rc.pmy_block;
+  int is = 0;
+  int js = 0;
+  int ks = 0;
+  int ie = pmb->ncells1 - 1;
+  int je = pmb->ncells2 - 1;
+  int ke = pmb->ncells3 - 1;
+  CellVariable<Real> &qin = rc.Get("one_minus_advected_sq");
+  // get component 12
+  CellVariable<Real> &q0 = rc.Get("one_minus_sqrt_one_minus_advected_sq", 12);
+  // and component 37
+  CellVariable<Real> &q1 = rc.Get("one_minus_sqrt_one_minus_advected_sq", 37);
+  for (int i = is; i <= ie; i++) {
+    for (int j = js; j <= je; j++) {
+      for (int k = ks; k <= ke; k++) {
+        // this will make component 12 = advected
+        q0(k, j, i) = 1.0 - sqrt(qin(k, j, i));
+        // and this will make component 37 = 1 - advected
+        q1(k, j, i) = 1.0 - q0(k, j, i);
+      }
+    }
+  }*/
+}
+
+// provide the routine that estimates a stable timestep for this package
+Real EstimateTimestep(Container<Real> &rc) {
+  /*MeshBlock *pmb = rc.pmy_block;
+  auto pkg = pmb->packages["Advection"];
+  const auto &cfl = pkg->Param<Real>("cfl");
+  const auto &vx = pkg->Param<Real>("vx");
+  const auto &vy = pkg->Param<Real>("vy");
+
+  int is = pmb->is;
+  int js = pmb->js;
+  int ks = pmb->ks;
+  int ie = pmb->ie;
+  int je = pmb->je;
+  int ke = pmb->ke;
+
+  Real min_dt = std::numeric_limits<Real>::max();
+  ParArrayND<Real> dx0("dx0", pmb->ncells1);
+  ParArrayND<Real> dx1("dx1", pmb->ncells1);
+
+  // this is obviously overkill for this constant velocity problem
+  for (int k = ks; k <= ke; k++) {
+    for (int j = js; j <= je; j++) {
+      pmb->pcoord->CenterWidth1(k, j, is, ie, dx0);
+      pmb->pcoord->CenterWidth2(k, j, is, ie, dx1);
+      for (int i = is; i <= ie; i++) {
+        min_dt = std::min(min_dt, dx0(i) / std::abs(vx));
+        min_dt = std::min(min_dt, dx1(i) / std::abs(vy));
+      }
+    }
+  }
+
+  return cfl * min_dt;*/
+
+  MeshBlock *pmb = rc.pmy_block;
+  auto pkg = pmb->packages["Particles"];
+  const auto &tlim = pkg->Param<Real>("tlim");
+
+  return 0.1*tlim;
+}
+
+// Compute fluxes at faces given the constant velocity field and
+// some field "advected" that we are pushing around.
+// This routine implements all the "physics" in this example
+TaskStatus CalculateFluxes(Container<Real> &rc) {
+  /*MeshBlock *pmb = rc.pmy_block;
+  int is = pmb->is;
+  int js = pmb->js;
+  int ks = pmb->ks;
+  int ie = pmb->ie;
+  int je = pmb->je;
+  int ke = pmb->ke;
+
+  CellVariable<Real> &q = rc.Get("advected");
+  auto pkg = pmb->packages["Advection"];
+  const auto &vx = pkg->Param<Real>("vx");
+  const auto &vy = pkg->Param<Real>("vy");
+
+  int maxdim = std::max(std::max(pmb->ncells1, pmb->ncells2), pmb->ncells3);
+  ParArrayND<Real> ql("ql", maxdim);
+  ParArrayND<Real> qr("qr", maxdim);
+  ParArrayND<Real> qltemp("qltemp", maxdim);
+
+  // get x-fluxes
+  for (int k = ks; k <= ke; k++) {
+    for (int j = js; j <= je; j++) {
+      // get reconstructed state on faces
+      pmb->precon->DonorCellX1(k, j, is - 1, ie + 1, q.data, ql, qr);
+      if (vx > 0.0) {
+        for (int i = is; i <= ie + 1; i++) {
+          q.flux[0](k, j, i) = ql(i) * vx;
+        }
+      } else {
+        for (int i = is; i <= ie + 1; i++) {
+          q.flux[0](k, j, i) = qr(i) * vx;
+        }
+      }
+    }
+  }
+  // get y-fluxes
+  if (pmb->pmy_mesh->ndim >= 2) {
+    for (int k = ks; k <= ke; k++) {
+      pmb->precon->DonorCellX2(k, js - 1, is, ie, q.data, ql, qr);
+      for (int j = js; j <= je + 1; j++) {
+        pmb->precon->DonorCellX2(k, j, is, ie, q.data, qltemp, qr);
+        if (vy > 0.0) {
+          for (int i = is; i <= ie; i++) {
+            q.flux[1](k, j, i) = ql(i) * vy;
+          }
+        } else {
+          for (int i = is; i <= ie; i++) {
+            q.flux[1](k, j, i) = qr(i) * vy;
+          }
+        }
+        auto temp = ql;
+        ql = qltemp;
+        qltemp = temp;
+      }
+    }
+  }
+
+  // TODO(jcd): implement z-fluxes
+  */
+  return TaskStatus::complete;
+}
+
+} // namespace Particles
+
+// *************************************************//
+// define the application driver. in this case,    *//
+// that just means defining the MakeTaskList       *//
+// function.                                       *//
+// *************************************************//
+// first some helper tasks
+TaskStatus UpdateContainer(MeshBlock *pmb, int stage,
+                           std::vector<std::string> &stage_name, Integrator *integrator) {
+  // const Real beta = stage_wghts[stage-1].beta;
+  const Real beta = integrator->beta[stage - 1];
+  Container<Real> &base = pmb->real_containers.Get();
+  Container<Real> &cin = pmb->real_containers.Get(stage_name[stage - 1]);
+  Container<Real> &cout = pmb->real_containers.Get(stage_name[stage]);
+  Container<Real> &dudt = pmb->real_containers.Get("dUdt");
+  parthenon::Update::AverageContainers(cin, base, beta);
+  parthenon::Update::UpdateContainer(cin, dudt, beta * pmb->pmy_mesh->dt, cout);
+  return TaskStatus::complete;
+}
+
+TaskStatus UpdateSwarm(MeshBlock *pmb, int stage, std::vector<std::string> &stage_name,
+                       Integrator *integrator) {
+  SwarmContainer &base = pmb->real_containers.GetSwarm();
+  return TaskStatus::complete;
+}
+
+// See the advection.hpp declaration for a description of how this function gets called.
+TaskList ParticleDriver::MakeTaskList(MeshBlock *pmb, int stage) {
+  TaskList tl;
+  // we're going to populate our list with multiple kinds of tasks
+  // these lambdas just clean up the interface to adding tasks of the relevant kinds
+  /*auto AddMyTask = [&tl, pmb, stage, this](BlockStageNamesIntegratorTaskFunc func,
+                                           TaskID dep) {
+    return tl.AddTask<BlockStageNamesIntegratorTask>(func, dep, pmb, stage, stage_name,
+                                                     integrator);
+  };
+  auto AddContainerTask = [&tl](ContainerTaskFunc func, TaskID dep, Container<Real> &rc) {
+    return tl.AddTask<ContainerTask>(func, dep, rc);
+  };
+  auto AddTwoContainerTask = [&tl](TwoContainerTaskFunc f, TaskID dep,
+                                   Container<Real> &rc1, Container<Real> &rc2) {
+    return tl.AddTask<TwoContainerTask>(f, dep, rc1, rc2);
+  };*/
+
+  TaskID none(0);
+  // first make other useful containers
+  printf("stage = %i\n", stage);
+  if (stage == 1) {
+    SwarmContainer &base = pmb->real_containers.GetSwarm();
+    pmb->real_containers.Add("particles", base);
+    //Container<Real> &base = pmb->real_containers.Get();
+    //pmb->real_containers.Add("dUdt", base);
+    //for (int i = 1; i < integrator->nstages; i++)
+    //  pmb->real_containers.Add(stage_name[i], base);
+  }
+
+  // pull out the container we'll use to get fluxes and/or compute RHSs
+  //Container<Real> &sc0 = pmb->real_containers.Get(stage_name[stage - 1]);
+  // pull out a container we'll use to store dU/dt.
+  // This is just -flux_divergence in this example
+  //Container<Real> &dudt = pmb->real_containers.Get("dUdt");
+  // pull out the container that will hold the updated state
+  // effectively, sc1 = sc0 + dudt*dt
+  //Container<Real> &sc1 = pmb->real_containers.Get(stage_name[stage]);
+
+  SwarmContainer &particles = pmb->real_containers.GetSwarm("particles");
+
+  Swarm &swarm = particles.Get("particles");
+
+  /*auto start_recv = AddContainerTask(Container<Real>::StartReceivingTask, none, sc1);
+
+  auto advect_flux = AddContainerTask(Advection::CalculateFluxes, none, sc0);
+
+  auto send_flux =
+      AddContainerTask(Container<Real>::SendFluxCorrectionTask, advect_flux, sc0);
+  auto recv_flux =
+      AddContainerTask(Container<Real>::ReceiveFluxCorrectionTask, advect_flux, sc0);
+
+  // compute the divergence of fluxes of conserved variables
+  auto flux_div =
+      AddTwoContainerTask(parthenon::Update::FluxDivergence, recv_flux, sc0, dudt)*/;
+
+  // Push all particles by v*dt
+  auto update_swarm = tl.AddTask<TwoSwarmTask>(parthenon::Update::TransportSwarm, none,
+                                               swarm, swarm);
+
+  // update ghost cells
+  /*auto send =
+      AddContainerTask(Container<Real>::SendBoundaryBuffersTask, update_container, sc1);
+  auto recv = AddContainerTask(Container<Real>::ReceiveBoundaryBuffersTask, send, sc1);
+  auto fill_from_bufs = AddContainerTask(Container<Real>::SetBoundariesTask, recv, sc1);
+  auto clear_comm_flags =
+      AddContainerTask(Container<Real>::ClearBoundaryTask, fill_from_bufs, sc1);
+
+  auto prolongBound = tl.AddTask<BlockTask>(
+      [](MeshBlock *pmb) {
+        pmb->pbval->ProlongateBoundaries(0.0, 0.0);
+        return TaskStatus::complete;
+      },
+      fill_from_bufs, pmb);
+
+  // set physical boundaries
+  auto set_bc = AddContainerTask(parthenon::ApplyBoundaryConditions, prolongBound, sc1);
+  */
+  // fill in derived fields
+  //auto fill_derived =
+  //    AddContainerTask(parthenon::FillDerivedVariables::FillDerived, set_bc, sc1);
+
+  // estimate next time step
+  if (stage == integrator->nstages) {
+    /*auto new_dt = AddContainerTask(
+        [](Container<Real> &rc) {
+          MeshBlock *pmb = rc.pmy_block;
+          pmb->SetBlockTimestep(parthenon::Update::EstimateTimestep(rc));
+          return TaskStatus::complete;
+        },
+        fill_derived, sc1);*/
+    printf("stage == integrator->nstages, I guess\n");
+
+    // Update refinement
+    /*if (pmesh->adaptive) {
+      auto tag_refine = tl.AddTask<BlockTask>(
+          [](MeshBlock *pmb) {
+            pmb->pmr->CheckRefinementCondition();
+            return TaskStatus::complete;
+          },
+          fill_derived, pmb);
+    }*/
+    // Purge stages -- this task isn't really required.  If we don't purge the containers
+    // then the next time we go to add them, they'll already exist and it will be a no-op.
+    // Not purging the stages is more performant, but if the "base" Container changes, we
+    // need to purge the stages and recreate the non-base containers or else there will be
+    // bugs, e.g. the containers for rk stages won't have the same variables as the "base"
+    // container, likely leading to strange errors and/or segfaults.
+    /*auto purge_stages = tl.AddTask<BlockTask>(
+        [](MeshBlock *pmb) {
+          pmb->real_containers.PurgeNonBase();
+          return TaskStatus::complete;
+        },
+        fill_derived, pmb);*/
+  }
+  return tl;
+}
+
+} // namespace particles_example
diff --git a/example/particles/particles.hpp b/example/particles/particles.hpp
new file mode 100644
index 0000000..214767a
--- /dev/null
+++ b/example/particles/particles.hpp
@@ -0,0 +1,99 @@
+//========================================================================================
+// (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef EXAMPLE_PARTICLES_PARTICLES_HPP_
+#define EXAMPLE_PARTICLES_PARTICLES_HPP_
+
+#include <memory>
+
+#include "driver/driver.hpp"
+#include "driver/multistage.hpp"
+#include "interface/container.hpp"
+#include "interface/state_descriptor.hpp"
+#include "mesh/mesh.hpp"
+#include "task_list/tasks.hpp"
+
+using parthenon::AmrTag;
+using parthenon::BaseTask;
+using parthenon::Container;
+using parthenon::SwarmContainer;
+using parthenon::Swarm;
+using parthenon::Mesh;
+using parthenon::MeshBlock;
+using parthenon::MultiStageBlockTaskDriver;
+using parthenon::Outputs;
+using parthenon::ParameterInput;
+using parthenon::Real;
+using parthenon::StateDescriptor;
+using parthenon::TaskID;
+using parthenon::TaskList;
+using parthenon::TaskStatus;
+
+namespace particles_example {
+
+class ParticleDriver : public MultiStageBlockTaskDriver {
+ public:
+  ParticleDriver(ParameterInput *pin, Mesh *pm, Outputs *pout)
+      : MultiStageBlockTaskDriver(pin, pm, pout) {}
+  // This next function essentially defines the driver.
+  // Call graph looks like
+  // main()
+  //   MultiStageBlockTaskDriver::Execute (driver.cpp)
+  //     MultiStageBlockTaskDriver::Step (driver.cpp)
+  //       DriverUtils::ConstructAndExecuteBlockTasks (driver.hpp)
+  //         AdvectionDriver::MakeTaskList (advection.cpp)
+  TaskList MakeTaskList(MeshBlock *pmb, int stage);
+};
+
+// demonstrate making a custom Task type
+/*using ContainerTaskFunc = std::function<TaskStatus(Container<Real> &)>;
+class ContainerTask : public BaseTask {
+ public:
+  ContainerTask(TaskID id, ContainerTaskFunc func, TaskID dep, Container<Real> rc)
+      : BaseTask(id, dep), _func(func), _cont(rc) {}
+  TaskStatus operator()() { return _func(_cont); }
+
+ private:
+  ContainerTaskFunc _func;
+  Container<Real> _cont;
+};
+*/
+using TwoSwarmTaskFunc =
+    std::function<TaskStatus(Swarm &, Swarm &)>;
+class TwoSwarmTask : public BaseTask {
+ public:
+  TwoSwarmTask(TaskID id, TwoSwarmTaskFunc func, TaskID dep, Swarm s1,
+                   Swarm s2)
+      : BaseTask(id, dep), _func(func), _swarm1(s1), _swarm2(s2) {}
+  TaskStatus operator()() { return _func(_swarm1, _swarm2); }
+
+ private:
+  TwoSwarmTaskFunc _func;
+  Swarm _swarm1;
+  Swarm _swarm2;
+};
+
+namespace Particles {
+
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
+AmrTag CheckRefinement(Container<Real> &rc);
+void PreFill(Container<Real> &rc);
+void SquareIt(Container<Real> &rc);
+void PostFill(Container<Real> &rc);
+Real EstimateTimestep(Container<Real> &rc);
+//TaskStatus CalculateFluxes(Container<Real> &rc);
+
+} // namespace Advection
+
+} // namespace advection_example
+
+#endif // EXAMPLE_ADVECTION_ADVECTION_HPP_
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index cd6af2f..6dee028 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -91,6 +91,8 @@ add_library(parthenon
 
   interface/container_collection.cpp
   interface/container.cpp
+  interface/swarm_container.cpp
+  interface/swarm.cpp
   interface/metadata.cpp
   interface/properties_interface.cpp
   interface/sparse_variable.cpp
diff --git a/src/driver/driver.cpp b/src/driver/driver.cpp
index 354f3f4..30d31ec 100644
--- a/src/driver/driver.cpp
+++ b/src/driver/driver.cpp
@@ -25,9 +25,12 @@ DriverStatus EvolutionDriver::Execute() {
   pmesh->mbcnt = 0;
   while ((pmesh->time < pmesh->tlim) &&
          (pmesh->nlim < 0 || pmesh->ncycle < pmesh->nlim)) {
+    printf("%s %i\n", __FILE__, __LINE__);
     if (Globals::my_rank == 0) pmesh->OutputCycleDiagnostics();
+    printf("%s %i\n", __FILE__, __LINE__);
 
     TaskListStatus status = Step();
+    printf("%s %i\n", __FILE__, __LINE__);
     if (status != TaskListStatus::complete) {
       std::cerr << "Step failed to complete all tasks." << std::endl;
       return DriverStatus::failed;
diff --git a/src/driver/driver.hpp b/src/driver/driver.hpp
index 51df104..399ce7e 100644
--- a/src/driver/driver.hpp
+++ b/src/driver/driver.hpp
@@ -58,28 +58,42 @@ namespace DriverUtils {
 
 template <typename T, class... Args>
 TaskListStatus ConstructAndExecuteBlockTasks(T *driver, Args... args) {
+  printf("%s %i\n", __FILE__, __LINE__);
 #ifdef OPENMP_PARALLEL
   int nthreads = driver->pmesh->GetNumMeshThreads();
 #endif
+  printf("%s %i\n", __FILE__, __LINE__);
   int nmb = driver->pmesh->GetNumMeshBlocksThisRank(Globals::my_rank);
+  printf("%s %i\n", __FILE__, __LINE__);
   std::vector<TaskList> task_lists;
+  printf("%s %i\n", __FILE__, __LINE__);
   MeshBlock *pmb = driver->pmesh->pblock;
+  printf("%s %i\n", __FILE__, __LINE__);
   while (pmb != nullptr) {
+  printf("%s %i\n", __FILE__, __LINE__);
     task_lists.push_back(driver->MakeTaskList(pmb, std::forward<Args>(args)...));
+  printf("%s %i\n", __FILE__, __LINE__);
     pmb = pmb->next;
   }
+  printf("%s %i\n", __FILE__, __LINE__);
   int complete_cnt = 0;
+  printf("%s %i\n", __FILE__, __LINE__);
   while (complete_cnt != nmb) {
+  printf("%s %i\n", __FILE__, __LINE__);
     // TODO(pgrete): need to let Kokkos::PartitionManager handle this
     for (auto i = 0; i < nmb; ++i) {
+  printf("%s %i\n", __FILE__, __LINE__);
       if (!task_lists[i].IsComplete()) {
+  printf("%s %i\n", __FILE__, __LINE__);
         auto status = task_lists[i].DoAvailable();
         if (status == TaskListStatus::complete) {
+  printf("%s %i\n", __FILE__, __LINE__);
           complete_cnt++;
         }
       }
     }
   }
+  printf("%s %i\n", __FILE__, __LINE__);
   return TaskListStatus::complete;
 }
 
diff --git a/src/driver/multistage.cpp b/src/driver/multistage.cpp
index 2c2efb4..4715a93 100644
--- a/src/driver/multistage.cpp
+++ b/src/driver/multistage.cpp
@@ -56,6 +56,7 @@ TaskListStatus MultiStageBlockTaskDriver::Step() {
   using DriverUtils::ConstructAndExecuteBlockTasks;
   TaskListStatus status;
   for (int stage = 1; stage <= integrator->nstages; stage++) {
+    printf("stage: %i/%i\n", stage, integrator->nstages);
     status = ConstructAndExecuteBlockTasks<>(this, stage);
     if (status != TaskListStatus::complete) break;
   }
diff --git a/src/interface/container_collection.cpp b/src/interface/container_collection.cpp
index 1e05b9e..7bbec1d 100644
--- a/src/interface/container_collection.cpp
+++ b/src/interface/container_collection.cpp
@@ -17,6 +17,29 @@
 
 namespace parthenon {
 
+template <typename T>
+void ContainerCollection<T>::Add(const std::string &name, SwarmContainer &src) {
+  // error check for duplicate names
+  auto it = swarmContainers_.find(name);
+  if (it != swarmContainers_.end()) {
+    // check to make sure they are the same
+    if (!(src == *(it->second))) {
+      throw std::runtime_error("Error attempting to add a SwarmContainer to a Collection");
+    }
+    return;
+  }
+
+  auto sc = std::make_shared<SwarmContainer>();
+  sc->pmy_block = src.pmy_block;
+  for (auto v : src.GetSwarmVector()) {
+    // just copy the (shared) pointer
+    // TODO BRR change this?
+    sc->Add(v);
+  }
+
+  swarmContainers_[name] = sc;
+}
+
 template <typename T>
 void ContainerCollection<T>::Add(const std::string &name, Container<T> &src) {
   // error check for duplicate names
diff --git a/src/interface/container_collection.hpp b/src/interface/container_collection.hpp
index 4420547..90a1acf 100644
--- a/src/interface/container_collection.hpp
+++ b/src/interface/container_collection.hpp
@@ -18,6 +18,7 @@
 #include <string>
 
 #include "interface/container.hpp"
+#include "interface/swarm_container.hpp"
 #include "interface/metadata.hpp"
 
 namespace parthenon {
@@ -26,11 +27,15 @@ template <typename T>
 class ContainerCollection {
  public:
   ContainerCollection() {
-    containers_["base"] = std::make_shared<Container<T>>(); // always add "base" container
+    // Always add "base" containers
+    containers_["base"] = std::make_shared<Container<T>>();
+    swarmContainers_["base"] = std::make_shared<SwarmContainer>();
   }
 
   void Add(const std::string &label, Container<T> &src);
 
+  void Add(const std::string &label, SwarmContainer &src);
+
   Container<T> &Get() { return *containers_["base"]; }
   Container<T> &Get(const std::string &label) {
     auto it = containers_.find(label);
@@ -40,6 +45,15 @@ class ContainerCollection {
     return *(it->second);
   }
 
+  SwarmContainer &GetSwarm() { return *swarmContainers_["base"]; }
+  SwarmContainer &GetSwarm(const std::string &label) {
+    auto it = swarmContainers_.find(label);
+    if (it == swarmContainers_.end()) {
+      throw std::runtime_error("SwarmContainer " + label + " does not exist in collection.");
+    }
+    return *(it->second);
+  }
+
   void PurgeNonBase() {
     auto c = containers_.begin();
     while (c != containers_.end()) {
@@ -49,6 +63,14 @@ class ContainerCollection {
         ++c;
       }
     }
+    auto sc = swarmContainers_.begin();
+    while (sc != swarmContainers_.end()) {
+      if (sc->first != "base") {
+        sc = swarmContainers_.erase(sc);
+      } else {
+        ++sc;
+      }
+    }
   }
 
   void Print() {
@@ -57,10 +79,16 @@ class ContainerCollection {
       c.second->Print();
       std::cout << std::endl;
     }
+    for (auto &sc : swarmContainers_) {
+      std::cout << "SwarmContainer " << sc.first << " has:" << std::endl;
+      sc.second->Print();
+      std::cout << std::endl;
+    }
   }
 
  private:
   std::map<std::string, std::shared_ptr<Container<T>>> containers_;
+  std::map<std::string, std::shared_ptr<SwarmContainer>> swarmContainers_;
 };
 
 } // namespace parthenon
diff --git a/src/interface/metadata.hpp b/src/interface/metadata.hpp
index a360974..66a120b 100644
--- a/src/interface/metadata.hpp
+++ b/src/interface/metadata.hpp
@@ -68,7 +68,13 @@
   /** Do boundary communication */                                                       \
   PARTHENON_INTERNAL_FOR_FLAG(FillGhost)                                                 \
   /** Communication arrays are a copy: hint to destructor */                             \
-  PARTHENON_INTERNAL_FOR_FLAG(SharedComms)
+  PARTHENON_INTERNAL_FOR_FLAG(SharedComms)                                               \
+  /** Integer-valued quantity */                                                         \
+  PARTHENON_INTERNAL_FOR_FLAG(Integer)                                                   \
+  /** Real-valued quantity */                                                            \
+  PARTHENON_INTERNAL_FOR_FLAG(Real)                                                      \
+  /** String-valued quantity */                                                          \
+  PARTHENON_INTERNAL_FOR_FLAG(String)                                                    \
 
 namespace parthenon {
 
@@ -198,6 +204,19 @@ class Metadata {
     return None;
   }
 
+  /// returns the type of the variable
+  MetadataFlag Type() const {
+    if (IsSet(Integer)) {
+      return Integer;
+    } else if (IsSet(Real)) {
+      return Real;
+    } else if (IsSet(String)) {
+      return String;
+    }
+    /// by default return Metadata::None
+    return None;
+  }
+
   int GetSparseId() const { return sparse_id_; }
 
   const std::vector<int> &Shape() const { return shape_; }
diff --git a/src/interface/state_descriptor.hpp b/src/interface/state_descriptor.hpp
index e534c28..dd1b2ca 100644
--- a/src/interface/state_descriptor.hpp
+++ b/src/interface/state_descriptor.hpp
@@ -19,6 +19,7 @@
 #include <vector>
 
 #include "interface/container.hpp"
+#include "interface/swarm.hpp"
 #include "interface/metadata.hpp"
 #include "interface/params.hpp"
 #include "refinement/amr_criteria.hpp"
@@ -57,6 +58,15 @@ class StateDescriptor {
   // retrieve label
   const std::string &label() { return _label; }
 
+  bool AddSwarm(const std::string &swarm_name, Metadata &m) {
+    if (_swarmMetadataMap.count(swarm_name) > 0) {
+      throw std::invalid_argument("Swarm " + swarm_name + "already exists!");
+    }
+    _swarmMetadataMap[swarm_name] = m;
+
+    return true;
+  }
+
   // field addition / retrieval routines
   // add a field with associated metadata
   bool AddField(const std::string &field_name, Metadata &m,
@@ -105,16 +115,32 @@ class StateDescriptor {
     return names;
   }
 
+  // retrieve all swarm names
+  std::vector<std::string> Swarms() {
+    std::vector<std::string> names;
+    names.reserve(_swarmMetadataMap.size());
+    for (auto &x : _swarmMetadataMap) {
+      names.push_back(x.first);
+    }
+    return names;
+  }
+
   const std::map<std::string, Metadata> &AllFields() { return _metadataMap; }
   const std::map<std::string, std::vector<Metadata>> &AllSparseFields() {
     return _sparseMetadataMap;
   }
+  const std::map<std::string, Metadata> &AllSwarms() { return _swarmMetadataMap; }
 
   // retrieve metadata for a specific field
   Metadata &FieldMetadata(const std::string &field_name) {
     return _metadataMap[field_name];
   }
 
+  // retrieve metadata for a specific swarm
+  Metadata &SwarmMetadata(const std::string &swarm_name) {
+    return _swarmMetadataMap[swarm_name];
+  }
+
   // get all metadata for this physics
   const std::map<std::string, Metadata> &AllMetadata() { return _metadataMap; }
 
@@ -128,6 +154,7 @@ class StateDescriptor {
   const std::string _label;
   std::map<std::string, Metadata> _metadataMap;
   std::map<std::string, std::vector<Metadata>> _sparseMetadataMap;
+  std::map<std::string, Metadata> _swarmMetadataMap;
 };
 
 using Packages_t = std::map<std::string, std::shared_ptr<StateDescriptor>>;
diff --git a/src/interface/swarm.cpp b/src/interface/swarm.cpp
new file mode 100644
index 0000000..ed32b44
--- /dev/null
+++ b/src/interface/swarm.cpp
@@ -0,0 +1,129 @@
+//========================================================================================
+// (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#include <cstdlib>
+#include <memory>
+#include <utility>
+#include <vector>
+#include "bvals/cc/bvals_cc.hpp"
+#include "swarm.hpp"
+#include "globals.hpp" // my_rank
+#include "mesh/mesh.hpp"
+
+namespace parthenon {
+
+void Swarm::Add(const std::vector<std::string> labelArray,
+                       const Metadata &metadata) {
+  // generate the vector and call Add
+  for (auto label : labelArray) {
+    Add(label, metadata);
+  }
+}
+
+///
+/// The internal routine for allocating a particle swarm.  This subroutine
+/// is topology aware and will allocate accordingly.
+///
+/// @param label the name of the variable
+/// @param metadata the metadata associated with the particle
+void Swarm::Add(const std::string label, const Metadata &metadata) {
+  // labels must be unique, even between different types of data
+  if (intMap_.count(label) > 0 ||
+      realMap_.count(label) > 0 ||
+      stringMap_.count(label) > 0) {
+    throw std::invalid_argument ("swarm variable " + label + " already enrolled during Add()!");
+  }
+
+  std::array<int, 6> arrDims = {nmax_pool_, 1, 1, 1, 1, 1};
+
+  if (metadata.Type() == Metadata::Integer) {
+    auto var = std::make_shared<ParticleVariable<int>>(label, nmax_pool_, metadata);
+    intVector_.push_back(var);
+    intMap_[label] = var;
+  } else if (metadata.Type() == Metadata::Real) {
+    auto var = std::make_shared<ParticleVariable<Real>>(label, nmax_pool_, metadata);
+    realVector_.push_back(var);
+    realMap_[label] = var;
+  } else if (metadata.Type() == Metadata::String) {
+    auto var = std::make_shared<ParticleVariable<std::string>>(label, nmax_pool_, metadata);
+    stringVector_.push_back(var);
+    stringMap_[label] = var;
+  } else {
+    throw std::invalid_argument ("swarm variable " + label + " does not have a valid type during Add()");
+  }
+}
+
+void Swarm::Remove(const std::string label) {
+  int idx, isize;
+  bool found = false;
+
+  // Find index of variable
+  idx = 0;
+  for (auto v : intVector_) {
+    if ( label == v->label() ) {
+      found = true;
+      break;
+    }
+    idx++;
+  }
+  if (found == true) {
+    // first delete the variable
+    intVector_[idx].reset();
+
+    // Next move the last element into idx and pop last entry
+    if ( intVector_.size() > 1) intVector_[idx] = std::move(intVector_.back());
+    intVector_.pop_back();
+
+    // Also remove variable from map
+    intMap_.erase(label);
+  }
+
+  if (found == false) {
+    idx = 0;
+    for (auto v : realVector_) {
+      if ( label == v->label() ) {
+        found = true;
+        break;
+      }
+      idx++;
+    }
+  }
+  if (found == true) {
+    realVector_[idx].reset();
+    if ( realVector_.size() > 1) realVector_[idx] = std::move(realVector_.back());
+    realVector_.pop_back();
+    realMap_.erase(label);
+  }
+
+  if (found == false) {
+    idx = 0;
+    for (auto v : stringVector_) {
+      if ( label == v->label() ) {
+        found = true;
+        break;
+      }
+      idx++;
+    }
+  }
+  if (found == true) {
+    stringVector_[idx].reset();
+    if ( stringVector_.size() > 1) stringVector_[idx] = std::move(stringVector_.back());
+    stringVector_.pop_back();
+    stringMap_.erase(label);
+  }
+
+  if (found == false) {
+    throw std::invalid_argument ("swarm not found in Remove()");
+  }
+}
+
+} // namespace parthenon
diff --git a/src/interface/swarm.hpp b/src/interface/swarm.hpp
new file mode 100644
index 0000000..3cf0894
--- /dev/null
+++ b/src/interface/swarm.hpp
@@ -0,0 +1,141 @@
+//========================================================================================
+// (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef INTERFACE_SWARM_HPP_
+#define INTERFACE_SWARM_HPP_
+
+///
+/// A swarm contains all particles of a particular species
+/// Date: August 21, 2019
+
+#include <array>
+#include <cstdint>
+#include <map>
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+#include "athena.hpp"
+#include "parthenon_arrays.hpp"
+#include "variable.hpp"
+#include "metadata.hpp"
+#include "bvals/cc/bvals_cc.hpp"
+
+namespace parthenon {
+class MeshBlock;
+
+enum class PARTICLE_TYPE {
+  INT, REAL, STRING
+};
+
+enum class PARTICLE_STATUS {
+  UNALLOCATED, ALIVE, DEAD
+};
+
+class Swarm {
+  public:
+    Swarm(const std::string label, const Metadata &metadata,
+          const int nmax_pool_in = 1000) :
+          label_(label),
+          m_(metadata),
+          nmax_pool_(nmax_pool_in),
+          mpiStatus(true) {
+      printf("%s %i\n", __FILE__, __LINE__);
+      Add("x", Metadata({Metadata::Real}));
+      Add("y", Metadata({Metadata::Real}));
+      Add("z", Metadata({Metadata::Real}));
+      Add("mask", Metadata({Metadata::Integer}));
+      auto &mask = GetInteger("mask");
+      for (int n = 0; n < nmax_pool_; n++) {
+        mask(n) = 0;
+      }
+    }
+
+  ///< Add variable to swarm
+  void Add(const std::string label, const Metadata &metadata);
+
+  ///< Add multiple variables with common metadata to swarm
+  void Add(const std::vector<std::string> labelVector, const Metadata &metadata);
+
+  ///< Remote a variable from swarm
+  void Remove(const std::string label);
+
+  ParticleVariable<Real> &GetReal(const std::string label) {
+    return *(realMap_.at(label));
+  }
+
+  ParticleVariable<int> &GetInteger(const std::string label) {
+    return *(intMap_.at(label));
+  }
+
+  ///< Assign label for swarm
+  void setLabel(const std::string label) { label_ = label; }
+
+  ///< retrieve label for swarm
+  std::string label() const { return label_; }
+
+  ///< retrieve metadata for swarm
+  const Metadata metadata() const { return m_; }
+
+  /// Assign info for swarm
+  void setInfo(const std::string info) { info_ = info; }
+
+  /// return information string
+  std::string info() { return info_; }
+
+  /// Set max pool size
+  void setPoolMax(const int nmax_pool) {
+    nmax_pool_ = nmax_pool;
+    // TODO(BRR) require that nmax_pool > nmax_pool_?
+    // TODO(BRR) resize arrays and copy data
+  }
+
+  int get_nmax_active() {
+    return nmax_active_;
+  }
+
+  bool mpiStatus;
+
+  void AddParticle() {
+    // Check that particle fits, if not double size of pool via
+    // setPoolMax(2*_nmax_pool);
+  }
+
+  // TODO(BRR) add a bunch of particles at once?
+
+  void Defrag() {
+    // TODO(BRR) Put a fast algorithm here to defrag memory
+  }
+
+ private:
+  int nmax_pool_;
+  int nmax_active_ = 0;
+  Metadata m_;
+  std::string label_;
+  std::string info_;
+  std::shared_ptr<ParArrayND<PARTICLE_STATUS>> pstatus_;
+  ParticleVariableVector<int> intVector_;
+  ParticleVariableVector<Real> realVector_;
+  ParticleVariableVector<std::string> stringVector_;
+
+  MapToParticle<int> intMap_;
+  MapToParticle<Real> realMap_;
+  MapToParticle<std::string> stringMap_;
+};
+
+using SP_Swarm = std::shared_ptr<Swarm>;
+using SwarmVector = std::vector<SP_Swarm>;
+using SwarmMap = std::map<std::string, SP_Swarm>;
+
+} // namespace parthenon
+
+#endif // INTERFACE_VARIABLE_HPP_
diff --git a/src/interface/swarm_container.cpp b/src/interface/swarm_container.cpp
new file mode 100644
index 0000000..4a21c9e
--- /dev/null
+++ b/src/interface/swarm_container.cpp
@@ -0,0 +1,96 @@
+//========================================================================================
+// (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#include <cstdlib>
+#include <memory>
+#include <utility>
+#include <vector>
+#include "bvals/cc/bvals_cc.hpp"
+#include "swarm_container.hpp"
+#include "globals.hpp" // my_rank
+#include "mesh/mesh.hpp"
+
+namespace parthenon {
+
+void SwarmContainer::Add(const std::vector<std::string> labelArray,
+                       const Metadata &metadata) {
+  // generate the vector and call Add
+  for (auto label : labelArray) {
+    Add(label, metadata);
+  }
+}
+
+///
+/// The internal routine for allocating a particle swarm.  This subroutine
+/// is topology aware and will allocate accordingly.
+///
+/// @param label the name of the variable
+/// @param metadata the metadata associated with the particle
+void SwarmContainer::Add(const std::string label,
+                       const Metadata &metadata) {
+  if (swarmMap_.find(label) != swarmMap_.end()) {
+    throw std::invalid_argument ("swarm " + label  +" already enrolled during Add()!");
+  }
+
+  auto swarm = std::make_shared<Swarm>(label, metadata);
+  swarmVector_.push_back(swarm);
+  swarmMap_[label] = swarm;
+}
+
+void SwarmContainer::Remove(const std::string label) {
+  int idx, isize;
+
+  // Find index of swarm
+  isize = swarmVector_.size();
+  idx = 0;
+  for (auto s : swarmVector_) {
+    if ( ! label.compare(s->label())) {
+      break;
+    }
+    idx++;
+  }
+  if ( idx >= isize) {
+    throw std::invalid_argument ("swarm not found in Remove()");
+  }
+
+  // first delete the variable
+  swarmVector_[idx].reset();
+
+  // Next move the last element into idx and pop last entry
+  isize--;
+  if ( isize >= 0) swarmVector_[idx] = std::move(swarmVector_.back());
+  swarmVector_.pop_back();
+
+  // Also remove swarm from map
+  swarmMap_.erase(label);
+}
+
+void SwarmContainer::SendBoundaryBuffers() {}
+
+void SwarmContainer::SetupPersistentMPI() {}
+
+bool SwarmContainer::ReceiveBoundaryBuffers() { return true; }
+
+void SwarmContainer::ReceiveAndSetBoundariesWithWait() {}
+
+void SwarmContainer::SetBoundaries() {}
+
+void SwarmContainer::StartReceiving(BoundaryCommSubset phase) {}
+
+void SwarmContainer::ClearBoundary(BoundaryCommSubset phase) {}
+
+void SwarmContainer::Print() {
+  std::cout << "Swarms are:\n";
+  for (auto s : swarmMap_) { std::cout << "  " << s.second->info() << std::endl; }
+}
+
+} // namespace parthenon
diff --git a/src/interface/swarm_container.hpp b/src/interface/swarm_container.hpp
new file mode 100644
index 0000000..1a0c8a8
--- /dev/null
+++ b/src/interface/swarm_container.hpp
@@ -0,0 +1,174 @@
+//========================================================================================
+// (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef INTERFACE_SWARM_CONTAINER_HPP_
+#define INTERFACE_SWARM_CONTAINER_HPP_
+
+#include <map>
+#include <memory>
+#include <string>
+#include <utility> // <pair>
+#include <vector>
+#include "globals.hpp"
+#include "swarm.hpp"
+
+namespace parthenon {
+///
+/// Interface to underlying infrastructure for particle data declaration and
+/// access.
+/// Date: August 22, 2019
+///
+///
+/// The SwarmContainer class is a container for the swarms of particles that
+/// make up the simulation.
+///
+/// The container class will provide the following methods:
+///
+
+class MeshBlock;
+
+class SwarmContainer {
+ public:
+  //-----------------
+  // Public Variables
+  //-----------------
+  MeshBlock* pmy_block = nullptr; // ptr to MeshBlock
+
+  //-----------------
+  //Public Methods
+  //-----------------
+  // Constructor does nothing
+  SwarmContainer() {
+  }
+
+  ///
+  /// Set the pointer to the mesh block for this container
+  void setBlock(MeshBlock *pmb) { pmy_block = pmb; }
+
+  // TODO BRR also add Add() functions for setting single int, real, string
+  // values?
+
+  ///
+  /// Allocate and add a variable<T> to the container
+  ///
+  /// This function will eventually look at the metadata flags to
+  /// identify the size of the first dimension based on the
+  /// topological location.  Dimensions will be taken from the metadata.
+  ///
+  /// @param label the name of the variable
+  /// @param metadata the metadata associated with the variable
+  ///
+  void Add(const std::string label, const Metadata &metadata);
+
+  ///
+  /// Allocate and add a variable<T> to the container
+  ///
+  /// This function will eventually look at the metadata flags to
+  /// identify the size of the first dimension based on the
+  /// topological location.  Dimensions will be taken from the metadata.
+  ///
+  /// @param labelArray the array of names of variables
+  /// @param metadata the metadata associated with the variable
+  ///
+  void Add(const std::vector<std::string> labelArray, const Metadata &metadata);
+
+  void Add(std::shared_ptr<Swarm> swarm) {
+    swarmVector_.push_back(swarm);
+    swarmMap_[swarm->label()] = swarm;
+  }
+
+  ///
+  /// Get a swarm from the container
+  /// @param label the name of the swarm
+  /// @return the Swarm if found or throw exception
+  Swarm& Get(std::string label) {
+    if (swarmMap_.count(label) == 0) {
+      throw std::invalid_argument (std::string("\n") +
+                                   std::string(label) +
+                                   std::string(" swarm not found in Get()\n") );
+    }
+    return *swarmMap_[label];
+  }
+
+  Swarm& Get(const int index) {
+    return *(swarmVector_[index]);
+  }
+
+  int Index(const std::string& label) {
+    for (int i = 0; i < swarmVector_.size(); i++) {
+      if (! swarmVector_[i]->label().compare(label)) return i;
+    }
+    return -1;
+  }
+
+  /// Gets an array of real variables from container.
+  /// @param names is the variables we want
+  /// @param indexCount a map of names to std::pair<index,count> for each name
+  /// @param sparse_ids if specified is list of sparse ids we are interested in.  Note
+  ///        that non-sparse variables specified are aliased in as is.
+  int GetSwarms(const std::vector<std::string>& names,
+                std::vector<Swarm>& sRet,
+                std::map<std::string,std::pair<int,int>>& indexCount);
+
+  const SwarmVector &GetSwarmVector() const { return swarmVector_; }
+  const SwarmMap &GetSwarmMap() const { return swarmMap_; }
+
+  ///
+  /// Remove a variable from the container or throw exception if not
+  /// found.
+  /// @param label the name of the variable to be deleted
+  void Remove(const std::string label);
+
+  // Temporary functions till we implement a *real* iterator
+
+  /// Print list of labels in container
+  void Print();
+
+  // return number of stored arrays
+  int Size() {return swarmVector_.size();}
+
+  // Element accessor functions
+  std::vector<std::shared_ptr<Swarm>>& allSwarms() {
+    return swarmVector_;
+  }
+
+  // Communication routines
+  void SetupPersistentMPI();
+  void SetBoundaries();
+  void SendBoundaryBuffers();
+  void ReceiveAndSetBoundariesWithWait();
+  bool ReceiveBoundaryBuffers();
+  void StartReceiving(BoundaryCommSubset phase);
+  void ClearBoundary(BoundaryCommSubset phase);
+
+  bool operator==(const SwarmContainer &cmp) {
+    // Test that labels of swarms are the same
+    std::vector<std::string> my_keys;
+    std::vector<std::string> cmp_keys;
+    for (auto &s : swarmMap_) {
+      my_keys.push_back(s.first);
+    }
+    for (auto &s : cmp.GetSwarmMap()) {
+      cmp_keys.push_back(s.first);
+    }
+    return my_keys == cmp_keys;
+  }
+
+ private:
+  int debug=0;
+
+  SwarmVector swarmVector_ = {};
+  SwarmMap swarmMap_ = {};
+};
+
+} // namespace parthenon
+#endif // INTERFACE_SWARM_CONTAINER_HPP_
diff --git a/src/interface/update.cpp b/src/interface/update.cpp
index 318b693..ebcc81e 100644
--- a/src/interface/update.cpp
+++ b/src/interface/update.cpp
@@ -96,6 +96,35 @@ TaskStatus FluxDivergence(Container<Real> &in, Container<Real> &dudt_cont) {
   return TaskStatus::complete;
 }
 
+//void TransportSwarm(Swarm &in, const Real dt, Swarm &out) {
+TaskStatus TransportSwarm(Swarm &in, Swarm &out) {
+  // TODO BRR put const real dt back in the args list
+  const Real dt = 0.1;
+  int nmax_active = in.get_nmax_active();
+
+  ParticleVariable<Real> &x_in = in.GetReal("x");
+  ParticleVariable<Real> &y_in = in.GetReal("y");
+  ParticleVariable<Real> &z_in = in.GetReal("z");
+  ParticleVariable<Real> &x_out = out.GetReal("x");
+  ParticleVariable<Real> &y_out = out.GetReal("y");
+  ParticleVariable<Real> &z_out = out.GetReal("z");
+  //ParticleVariable<Real> vx = in.Get("vx");
+  //ParticleVariable<Real> vy = in.Get("vy");
+  //ParticleVariable<Real> vz = in.Get("vz");
+  //ParticleVariable<int> mask = in.Get("mask");
+  double vx = 1.;
+  double vy = 1.;
+  double vz = 1.;
+
+  for (int n = 0; n < nmax_active; n++) {
+    x_out(n) = x_in(n) + vx*dt;
+    y_out(n) = y_in(n) + vy*dt;
+    z_out(n) = z_in(n) + vz*dt;
+  }
+
+  return TaskStatus::complete;
+}
+
 void UpdateContainer(Container<Real> &in, Container<Real> &dudt_cont, const Real dt,
                      Container<Real> &out) {
   MeshBlock *pmb = in.pmy_block;
diff --git a/src/interface/update.hpp b/src/interface/update.hpp
index baccfbf..f8a1a1b 100644
--- a/src/interface/update.hpp
+++ b/src/interface/update.hpp
@@ -24,6 +24,7 @@ namespace Update {
 TaskStatus FluxDivergence(Container<Real> &in, Container<Real> &dudt_cont);
 void UpdateContainer(Container<Real> &in, Container<Real> &dudt_cont, const Real dt,
                      Container<Real> &out);
+TaskStatus TransportSwarm(Swarm &in, Swarm &out);
 void AverageContainers(Container<Real> &c1, Container<Real> &c2, const Real wgt1);
 Real EstimateTimestep(Container<Real> &rc);
 
diff --git a/src/interface/variable.cpp b/src/interface/variable.cpp
index f2e16c4..e02dfb7 100644
--- a/src/interface/variable.cpp
+++ b/src/interface/variable.cpp
@@ -155,6 +155,24 @@ std::string EdgeVariable<T>::info() {
   return s;
 }
 
+template <typename T>
+std::string ParticleVariable<T>::info() {
+  char tmp[100] = "";
+
+  // first add label
+  std::string s = this->label();
+  s.resize(20, '.');
+
+  // now append size
+  snprintf(tmp, sizeof(tmp), "%d", data.GetDim(1));
+  s += std::string(tmp);
+
+  // now append flag
+  s += " : " + this->metadata().MaskAsString();
+
+  return s;
+}
+
 template class CellVariable<Real>;
 template class FaceVariable<Real>;
 template class EdgeVariable<Real>;
diff --git a/src/interface/variable.hpp b/src/interface/variable.hpp
index de8c388..a94b0bb 100644
--- a/src/interface/variable.hpp
+++ b/src/interface/variable.hpp
@@ -199,6 +199,39 @@ class EdgeVariable {
   std::string label_;
 };
 
+template <typename T>
+class ParticleVariable {
+ public:
+  /// Initialize a particle variable
+  ParticleVariable(const std::string label, const int npool, const Metadata &metadata)
+      : data(label, npool),
+        npool_(npool), m_(metadata), label_(label) {}
+
+  // accessors
+  template <class... Args>
+  KOKKOS_FORCEINLINE_FUNCTION auto &operator()(Args... args) {
+    return data(std::forward<Args>(args)...);
+  }
+
+  ///< retrieve metadata for variable
+  const Metadata metadata() const { return m_; }
+
+  bool IsSet(const MetadataFlag bit) const { return m_.IsSet(bit); }
+
+  ///< retrieve label for variable
+  std::string label() { return label_; }
+
+  /// return information string
+  std::string info();
+
+  ParArrayND<T> data;
+
+ private:
+  int npool_;
+  Metadata m_;
+  std::string label_;
+};
+
 template <typename T>
 using CellVariableVector = std::vector<std::shared_ptr<CellVariable<T>>>;
 template <typename T>
@@ -209,6 +242,11 @@ using MapToCellVars = std::map<std::string, std::shared_ptr<CellVariable<T>>>;
 template <typename T>
 using MapToFace = std::map<std::string, std::shared_ptr<FaceVariable<T>>>;
 
+template <typename T>
+using ParticleVariableVector = std::vector<std::shared_ptr<ParticleVariable<T>>>;
+template <typename T>
+using MapToParticle = std::map<std::string, std::shared_ptr<ParticleVariable<T>>>;
+
 } // namespace parthenon
 
 #endif // INTERFACE_VARIABLE_HPP_
diff --git a/src/mesh/meshblock.cpp b/src/mesh/meshblock.cpp
index 4a2cc9c..1e05a64 100644
--- a/src/mesh/meshblock.cpp
+++ b/src/mesh/meshblock.cpp
@@ -85,6 +85,7 @@ MeshBlock::MeshBlock(int igid, int ilid, LogicalLocation iloc, RegionSize input_
   }
 
   Container<Real> &real_container = real_containers.Get();
+  SwarmContainer &swarm_container = real_containers.GetSwarm();
 
   // Set the block pointer for the containers
   real_container.setBlock(this);
@@ -137,6 +138,19 @@ MeshBlock::MeshBlock(int igid, int ilid, LogicalLocation iloc, RegionSize input_
     }
   }
 
+  // Add swarm properties data
+  for (int i = 0; i < properties.size(); i++) {
+    StateDescriptor &state = properties[i]->State();
+    for (auto const &q : state.AllSwarms()) {
+      swarm_container.Add(q.first, q.second);
+    }
+  }
+  for (auto const &pkg : packages) {
+    for (auto const &q: pkg.second->AllSwarms()) {
+      swarm_container.Add(q.first, q.second);
+    }
+  }
+
   // TODO(jdolence): Should these loops be moved to Variable creation
   ContainerIterator<Real> ci(real_container, {Metadata::Independent});
   int nindependent = ci.vars.size();