subgrid is functional

31 files changed, 2412 insertions, 1613 deletions

diff --git a/LICENSE b/LICENSE
deleted file mode 100644
index a9d1f99..0000000
--- a/LICENSE
+++ /dev/null
@@ -1,25 +0,0 @@
-BSD 2-Clause License
-
-Copyright (c) 2021, Ilhan Özgen Xian
-All rights reserved.
-
-Redistribution and use in source and binary forms, with or without
-modification, are permitted provided that the following conditions are met:
-
-1. Redistributions of source code must retain the above copyright notice, this
-   list of conditions and the following disclaimer.
-
-2. Redistributions in binary form must reproduce the above copyright notice,
-   this list of conditions and the following disclaimer in the documentation
-   and/or other materials provided with the distribution.
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
-FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
-SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README b/README
new file mode 100644
index 0000000..cdbc310
--- /dev/null
+++ b/README
@@ -0,0 +1,63 @@
+
+
+                               /|  |\
+                                //\\
+
+
+			      tracksuit
+			      =========
+
+Contents
+
+1. About
+2. Dependencies
+3. Installation
+4. Usage
+5. Acknowledgment
+
+
+				About
+				-----
+
+/tracksuit/ is a subgrid-based particle tracker for integrated
+hydrology simulations.  it is specifically targeting the output of the
+hydrology solver amanzi/ats.
+
+
+			     Dependencies
+			     ------------
+
++ /tracksuit/ uses syntax that compiles for c++ version 14 or higher.
++ it requires hdf5 library with c++ support.
++ it uses openmp for shared memory parallelisation.
+
+
+			     Installation
+			     ------------
+
+/tracksuit/ ships with a makefile. adjust the path pointing to the
+hdf5 library. in the source directory, running
+
+$ make; make clean
+
+will produce the executable `ptracker'.
+
+
+				Usage
+				-----
+
+the program is run with
+
+$ ./ptracker -[a|c] <path/to/option/file>
+
+-c (c)ompiles the mesh, creating the output grid.h5
+-a (a)dvects particles on the subgrid
+
+			   Acknowledgments
+			   ---------------
+
+This work is supported as part of the Interoperable Design of
+Extreme-scale Application Software - Watersheds project
+(IDEAS-Watersheds) funded by the U.S. Department of Energy, Office of
+Science, Office of Biological and Environmental Research under Award
+Number DE-AC02-05CH11231.
diff --git a/README.md b/README.md
deleted file mode 100644
index 4c3b452..0000000
--- a/README.md
+++ /dev/null
@@ -1,24 +0,0 @@
-# ægent provocateur
-
-```
-
-   /  /^    agent provocateur
-  o--o      =================
-
-```
-
-_agent provocateur_ is a HydroAgent realisation that uses
-{[Amanzi/ATS](https://amanzi.github.io)} model results for the
-environment variables.
-
-## Dependencies
-
-- HDF5
-- MPI
-
-## References
-
-- S.M. Reaney, 2007. The use of agent based modelling techniques in
-  hydrology: determining the spatial and temporal origin of channel
-  flow in semi-arid catchments, _Earth Surface Processes and
-  Landforms_, 33, 317-327.
diff --git a/agent.cpp b/agent.cpp
deleted file mode 100644
index 2a5a183..0000000
--- a/agent.cpp
+++ /dev/null
@@ -1,266 +0,0 @@
-// -*- mode: c++ -*-
-
-#include "agent.hpp"
-
-// ------------------------------------------------------------------------
-// -- PUBLIC METHODS
-// ------------------------------------------------------------------------
-
-void Agent::advance(Mesh* mesh, double rel_dt_surface, double rel_dt_subsurface, size_t step)
-{
-
-  switch(domain[step-1]) {
-  case SURFACE_FLAG:
-    advance_surface(mesh, rel_dt_surface, step);
-    break;
-  case SUBSURFACE_FLAG:
-    advance_subsurface(mesh, rel_dt_subsurface, step);
-    break;
-  default:
-    set_path(path[step-1][0], path[step-1][1], path[step-1][2], step);
-    time[step] = time[step-1];
-    domain[step] = OUTSIDE_FLAG;
-    break;
-  }
-  
-}
-
-// ------------------------------------------------------------------------
-// GETTERS AND SETTERS
-// ------------------------------------------------------------------------
-
-void Agent::init(int n_iteration)
-{
-  time.resize(n_iteration);
-  domain.resize(n_iteration);
-  path.resize(n_iteration);
-  for (int i = 0; i < n_iteration; i ++)
-    path[i].resize(DIM);
-}
-
-void Agent::set_cell_index(int value)
-{
-  cell_index = value;
-}
-
-void Agent::set_domain(int value, size_t step)
-{
-  domain[step] = value;
-}
-
-double Agent::get_time(size_t step)
-{
-  return time[step];
-}
-
-void Agent::set_path(double x, double y, double z, size_t step)
-{
-  path[step][0] = x;
-  path[step][1] = y;
-  path[step][2] = z;
-}
-
-double Agent::get_timestepsize()
-{
-  return time_step_size;
-}
-
-void Agent::set_timestepsize(double value)
-{
-  time_step_size = value;
-}
-
-double Agent::get_path(size_t step, size_t index)
-{
-  return path[step][index];
-}
-
-int Agent::get_domain(size_t step)
-{
-  return domain[step];
-}
-
-// ------------------------------------------------------------------------
-// -- PRIVATE METHODS
-// ------------------------------------------------------------------------
-
-void Agent::advance_surface(Mesh* mesh, double rel_dt, size_t step)
-{
-
-#if CLASSIC_SCHEME
-  
-  Cell *location = mesh->get_cell(cell_index, SURFACE_FLAG);
-
-  std::vector<double> velocity(DIM);
-  std::vector<double> velocity_slope(DIM);
-  std::vector<double> coordinates(DIM-1);
-  std::vector<double> coordinates_neighbor(DIM-1);
-  std::vector<double> distance(DIM-1);
-  
-  domain[step] = SURFACE_FLAG; // set domain for this step, may change later if agent infiltrates
-
-  for (size_t i = 0; i < DIM; i ++) {
-    velocity[i] = location->get_velocity(i);
-    velocity_slope[i] = rel_dt * (location->get_velocity_next(i) - velocity[i]);
-    velocity[i] = velocity[i] + velocity_slope[i];
-  }
-
-  // CFL criteria
-  double length = location->get_characteristic_length();
-  double velocity_norm = sqrt(velocity[0]*velocity[0] + velocity[1]*velocity[1]);
-
-  double max_time_step = time_step_size; // time step to next time level
-  time_step_size = std::min(CFL * length / (velocity_norm + EPS), MAXTIMESTEP);
-  if (max_time_step > 0)
-    time_step_size = std::min(time_step_size, max_time_step);
-
-  time[step] = time[step-1] + time_step_size;
-
-  // check if agent infiltrates
-  double h  = location->get_water_depth();
-  int flag  = INFILTRATE_FALSE;
-  
-  if (h > EPS) { // only wet cell can infiltrate
-    double qe = std::min(velocity[2], 0.0);
-    double pi = fabs(qe) * time_step_size / h;
-    double trial = (double) (rand()%100) / 100.0;
-    if (trial < pi)
-      flag = INFILTRATE_TRUE;
-  }
-
-  if (flag == INFILTRATE_FALSE) {
-    
-    // advect the agent horizontally, which means the last dimension
-    // (z-direction) is omitted
-    for (size_t i = 0; i < DIM-1; i ++) {
-      coordinates[i] = path[step-1][i] + time_step_size * velocity[i];
-      distance[i] = location->get_coordinate(i) - coordinates[i];
-    }
-
-    // horizontal distance to the centroid of current cell
-    double distance_norm = sqrt(distance[0]*distance[0]+distance[1]*distance[1]);
-
-    // check if cell has changed
-    size_t n_cells = mesh->get_n_cells(SURFACE_FLAG);
-    for (size_t i = 0; i < n_cells; i ++) {
-
-      for (size_t j = 0; j < DIM-1; j ++) {
-	coordinates_neighbor[j] = mesh->get_cell(i, SURFACE_FLAG)->get_coordinate(j);
-	distance[j] = coordinates_neighbor[j] - coordinates[j];
-      }
-
-      double distance_norm_ = sqrt(distance[0]*distance[0]+distance[1]*distance[1]);
-    
-      if (distance_norm > distance_norm_) {
-	cell_index = i;
-	distance_norm = distance_norm_;
-      }
-    
-    }
-
-    // check if out of the domain: if the agent travelled thrice the
-    // maximum length of the cell and still can't find a neighboring
-    // cell, it must be advected out of the domain
-    if (distance_norm > 15.0 * location->get_maximum_length())
-      domain[step] = OUTSIDE_FLAG;
-  
-  } else {
-    // infiltration detected
-    domain[step] = SUBSURFACE_FLAG;
-    coordinates[0] = path[step-1][0];
-    coordinates[1] = path[step-1][1];
-    cell_index = mesh->get_cell(cell_index, SURFACE_FLAG)->get_coupled_neighbor();
-  }
-
-  set_path(coordinates[0], coordinates[1], location->get_coordinate(2), step);
-
-#endif
-  
-}
-
-
-void Agent::advance_subsurface(Mesh* mesh, double rel_dt, size_t step)
-{
-
-#if CLASSIC_SCHEME
-  
-  Cell *location = mesh->get_cell(cell_index, SUBSURFACE_FLAG);
-  
-  std::vector<double> velocity(DIM);
-  std::vector<double> velocity_slope(DIM);
-  std::vector<double> coordinates(DIM);
-  std::vector<double> coordinates_neighbor(DIM);
-  std::vector<double> distance(DIM);
-
-  domain[step] = SUBSURFACE_FLAG; // set domain for this step, may change later if agent infiltrates
-  
-  for (size_t i = 0; i < DIM; i ++) {
-    velocity[i] = location->get_velocity(i);
-    velocity_slope[i] = rel_dt * (location->get_velocity_next(i) - velocity[i]);
-    velocity[i] = velocity[i] + velocity_slope[i];
-  }
-
-  // CFL criteria
-  double length  = location->get_characteristic_length();
-  double velocity_norm = sqrt(velocity[0] * velocity[0] + velocity[1] * velocity[1] + velocity[2] * velocity[2]);
-
-  double max_time_step = time_step_size;
-  time_step_size = std::min(CFL * length / (velocity_norm + EPS), MAXTIMESTEP);
-  if (max_time_step > 0)
-    time_step_size = std::min(time_step_size, max_time_step);
-
-  time[step] = time[step-1] + time_step_size;
-
-  // advect the agent
-  for (size_t i = 0; i < DIM; i ++) {
-    coordinates[i] = path[step-1][i] + time_step_size * velocity[i];
-    distance[i] = location->get_coordinate(i) - coordinates[i];
-  }
-  
-  set_path(coordinates[0], coordinates[1], coordinates[2], step);
-
-  // distance to the centroid of current cell
-  double distance_norm = sqrt(distance[0] * distance[0] + distance[1] * distance[1] + distance[2] * distance[2]);
-  
-  // check if cell has changed
-  size_t n_cells = mesh->get_n_cells(SUBSURFACE_FLAG);
-  for (size_t i = 0; i < n_cells; i ++) {
-
-    for (size_t j = 0; j < DIM; j ++) {
-      coordinates_neighbor[j] = mesh->get_cell(i, SUBSURFACE_FLAG)->get_coordinate(j);
-      distance[j] = coordinates_neighbor[j] - coordinates[j];
-    }
-
-    double distance_norm_ = sqrt(distance[0] * distance[0] + distance[1] * distance[1] + distance[2] * distance[2]);
-
-    if (distance_norm > distance_norm_) {
-      cell_index = i;
-      distance_norm = distance_norm_;
-    }
-    
-  }
-
-  // check if out of the domain: if the agent travelled thrice the
-  // maximum length of the cell and still can't find a neighboring
-  // cell, it must be advected out of the domain
-  //  if (distance_norm > 15.0 * location->get_maximum_length()) {
-  //    domain[step] = OUTSIDE_FLAG;
-  //  }
-
-  // check if exfiltrates
-  int neighbor_index = mesh->get_cell(cell_index, SUBSURFACE_FLAG)->get_coupled_neighbor();
-  double z_surface;
-  if (neighbor_index > -1) {
-    // border to surface detected
-    z_surface = mesh->get_cell(neighbor_index, SURFACE_FLAG)->get_coordinate(2);
-    if (coordinates[2] >= z_surface) {
-      // exfiltration detected
-      cell_index = (int) neighbor_index;
-      domain[step] = SURFACE_FLAG;
-      set_path(coordinates[0], coordinates[1], z_surface, step);
-    }
-  }
-  
-#endif
-  
-}
diff --git a/agent.hpp b/agent.hpp
deleted file mode 100644
index c21e06a..0000000
--- a/agent.hpp
+++ /dev/null
@@ -1,43 +0,0 @@
-// -*- mode: c++ -*-
-
-#ifndef AGENT_HPP
-#define AGENT_HPP
-
-#include "header.hpp"
-#include "cell.hpp"
-#include "mesh.hpp"
-
-class Agent
-{
-
-public:
-
-  //  Agent();
-
-  void init(int n_iteration);
-  void set_cell_index(int value);
-  void set_domain(int value, size_t step);
-  void set_path(double x, double y, double z, size_t step);
-  double get_path(size_t step, size_t index);
-  int get_domain(size_t step);
-  double get_timestepsize();
-  void set_timestepsize(double value);
-  double get_time(size_t step);
-
-  void advance(Mesh* mesh, double rel_dt_surface, double rel_dt_subsurface, size_t step);
-
-private:
-
-  int cell_index;
-  double time_step_size;
-
-  std::vector<int> domain;
-  std::vector<double> time;
-  std::vector<std::vector<double>> path;
-
-  void advance_surface(Mesh *mesh, double rel_dt, size_t step);
-  void advance_subsurface(Mesh *mesh, double rel_dt, size_t step);
-  
-};
-
-#endif
diff --git a/cell.cpp b/cell.cpp
deleted file mode 100644
index 0f18478..0000000
--- a/cell.cpp
+++ /dev/null
@@ -1,149 +0,0 @@
-// -*- mode: c++ -*-
-
-#include "cell.hpp"
-
-// ------------------------------------------------------------------------
-// constructors
-// ------------------------------------------------------------------------
-
-Cell::Cell()
-{
-  
-  cell_id = -1;
-  coupled_neighbor = -1;
-  characteristic_length = -1.0;
-  maximum_length = -1.0;
-  
-  water_depth = 0.0;
-  water_depth_next = 0.0;
-
-  velocity.resize(DIM);
-  velocity_next.resize(DIM);
-  coordinates.resize(DIM);
-
-  neighbors.resize(N_NEIGHBORS);
-  
-  for (size_t i = 0; i < DIM; i ++) {
-    coordinates[i] = 0.0;
-    velocity[i] = 0.0;
-    velocity_next[i] = 0.0;
-  }
-
-  for (size_t i = 0; i < N_NEIGHBORS; i ++)
-    neighbors[i] = -1;
-  
-}
-
-Cell::Cell(double x, double y, double z)
-{
-
-  cell_id = -1;
-  coupled_neighbor = -1;
-  characteristic_length = -1.0;
-  maximum_length = -1.0;
-
-  water_depth = 0.0;
-  water_depth_next = 0.0;
-
-  velocity.resize(DIM);
-  velocity_next.resize(DIM);
-  coordinates.resize(DIM);
-  neighbors.resize(N_NEIGHBORS);
-
-  coordinates[0] = x;
-  coordinates[1] = y;
-  coordinates[2] = z;
-  
-  for (size_t i = 0; i < DIM; i ++) {
-    velocity[i] = 0.0;
-    velocity_next[i] = 0.0;
-  }
-
-  for (size_t i = 0; i < N_NEIGHBORS; i ++)
-    neighbors[i] = -1;
-  
-}
-
-// ------------------------------------------------------------------------
-// methods
-// ------------------------------------------------------------------------
-
-hsize_t Cell::get_neighbor(int index) {
-  return neighbors[index];
-}
-
-void Cell::set_neighbor(int index, hsize_t value) {
-  neighbors[index] = value;
-}
-
-void Cell::set_cell_id(int value) {
-  cell_id = value;
-}
-
-int Cell::get_cell_id() {
-  return cell_id;
-}
-
-void Cell::set_coupled_neighbor(hsize_t value) {
-  coupled_neighbor = value;
-}
-
-hsize_t Cell::get_coupled_neighbor() {
-  return coupled_neighbor;
-}
-  
-void Cell::set_water_depth(double value) {
-  water_depth = value;
-}
-
-void Cell::set_water_depth_next(double value) {
-  water_depth_next = value;
-}
-
-void Cell::set_velocity(int index, double value) {
-  velocity[index] = value;
-}
-
-void Cell::set_velocity_next(int index, double value) {
-  velocity_next[index] = value;
-}
-
-double Cell::get_coordinate(int index) {
-  return coordinates[index];
-}
-
-void Cell::set_coordinate(int index, double value) {
-  coordinates[index] = value;
-}
-
-double Cell::get_water_depth() {
-  return water_depth;
-}
-
-double Cell::get_water_depth_next() {
-  return water_depth_next;
-}
-
-double Cell::get_velocity(int index) {
-  return velocity[index];
-}
-
-double Cell::get_velocity_next(int index) {
-  return velocity_next[index];
-}
-
-double Cell::get_characteristic_length() {
-  return characteristic_length;
-}
-
-void Cell::set_characteristic_length(double value) {
-  characteristic_length = value;
-}
-
-double Cell::get_maximum_length() {
-  return maximum_length;
-}
-
-void Cell::set_maximum_length(double value) {
-  maximum_length = value;
-}
diff --git a/cell.hpp b/cell.hpp
deleted file mode 100644
index 37d4d0e..0000000
--- a/cell.hpp
+++ /dev/null
@@ -1,67 +0,0 @@
-// -*- mode: c++ -*-
-
-#ifndef CELL_HPP
-#define CELL_HPP
-
-#include "header.hpp"
-
-class Cell
-{
-  
-public:
-
-  Cell();
-  Cell(double x, double y, double z);
-
-  hsize_t get_neighbor(int index);
-  void set_neighbor(int index, hsize_t value);
-
-  void set_cell_id(int value);
-  int get_cell_id();
-  
-  void set_coupled_neighbor(hsize_t value);
-  hsize_t get_coupled_neighbor();
-
-  void set_water_depth(double value);
-  void set_velocity(int index, double value);
-  void set_water_depth_next(double value);
-  void set_velocity_next(int index, double value);
-
-  double get_water_depth();
-  double get_velocity(int index);
-  double get_water_depth_next();
-  double get_velocity_next(int index);
-
-  double get_characteristic_length();
-  void set_characteristic_length(double value);
-
-  double get_maximum_length();
-  void set_maximum_length(double value);
-  
-  double get_coordinate(int index);
-  void set_coordinate(int index, double value);
-
-private:
-
-  int cell_id;
-  int coupled_neighbor;		        // id of the neighbor cell that's in the
-					// other domain, for example the surface
-					// neighbor of a subsurface cell
-  
-  double characteristic_length;		// length used to compute CFL number
-  double maximum_length;                // length to check if out of domain
-  double water_depth;			// unused if subsurface cell
-  double water_depth_next;
-  
-  std::vector<double> velocity;		// either surface velocity or
-  std::vector<double> velocity_next;     // subsurface velocity
-  
-  std::vector<double> coordinates;
-  std::vector<hsize_t> neighbors;	// default to number of
-					// 3D neighbors, last 2
-					// entries unused if
-					// surface cell
-  
-};
-
-#endif
diff --git a/column.hpp b/column.hpp
new file mode 100644
index 0000000..f612458
--- /dev/null
+++ b/column.hpp
@@ -0,0 +1,61 @@
+// -*- mode: c++ -*-
+
+#ifndef COLUMN_HPP
+#define COLUMN_HPP
+
+#include "header.hpp"
+
+// Content
+// =======
+//
+// 1. Column
+// ---
+// 1.1 sort
+
+// ------------------------------------------------------------------------
+// Column
+// ------------------------------------------------------------------------
+// The column is a class with a top surface cell and its stacked
+// subsurface cells. The subsurface cells are sorted from top to bottom.
+// ------------------------------------------------------------------------
+class Column
+{
+
+public:
+
+    std::vector<point> vertices;
+    point centroid;
+
+    std::vector<int> stack_index; // stacked subsurface cell indices
+    std::vector<double> stack_lower_horizon; // minimum elevation of
+    // stacked subsurface cells
+    std::vector<double> stack_upper_horizon; // maximum elevation of
+    // stacked subsurface cells
+
+    double area;
+    //  double columndepth;
+
+    void sort()
+    {
+
+        // sort indices and upper horizons
+        std::vector<std::pair<int,double>> zipped;
+        zip(stack_index, stack_upper_horizon, zipped);
+        std::sort(std::begin(zipped), std::end(zipped),
+        [&](const auto& a, const auto& b) {
+            return a.second > b.second;
+        });
+        unzip(stack_index, stack_upper_horizon, zipped);
+        zipped.clear();
+        // sort the lower horizons
+        std::sort(std::begin(stack_lower_horizon),
+                  std::end(stack_lower_horizon),
+        [&](const auto& a, const auto& b) {
+            return a > b;
+        });
+
+    }
+
+};
+
+#endif
diff --git a/environment.hpp b/environment.hpp
new file mode 100644
index 0000000..fb2eec6
--- /dev/null
+++ b/environment.hpp
@@ -0,0 +1,350 @@
+// -*- mode: c++ -*-
+
+#ifndef ENVIRONMENT_HPP
+#define ENVIRONMENT_HPP
+
+#include "header.hpp"
+#include "parser.hpp"
+
+// Content
+// =======
+//
+// 1. Environment
+//
+
+// ------------------------------------------------------------------------
+// Environment
+// ------------------------------------------------------------------------
+// The environment reads and stores hydrological values.  For the
+// surface, it stores infiltration and evapotranspiration fluxes,
+// water depth, and velocities.  For the subsurface it stores
+// velocities and water content.
+//
+// The environment variables can be updated by reading in a different
+// file at a different time step.
+// ------------------------------------------------------------------------
+class Environment
+{
+
+public:
+
+    int nx;
+    int ny;
+    int nz;
+
+    double xll_corner;
+    double yhl_corner;
+
+    index3d indices;   // holds the subgrid to mesh map
+    grid3d elevations; // holds the topographic elevation map
+
+    std::vector<double> depths; // holds the vertical resolution
+
+    // surface variables /these have the length of the cell # in the surface mesh/
+    std::vector<double> exchange_flux;
+    std::vector<double> evaporation_flux;
+    std::vector<double> cell_area;
+    std::vector<double> water_depth;
+    std::vector<double> velocity_x;
+    std::vector<double> velocity_y;
+
+    // subsurface variables /these have the length of the cell # in the subsurface mesh/
+    std::vector<double> darcy_velocity_x;
+    std::vector<double> darcy_velocity_y;
+    std::vector<double> darcy_velocity_z;
+
+    Environment()
+    {
+
+        std::cout << "[env] empty environment initialised" << std::endl;
+
+    }
+
+    Environment(Parser &value)
+    {
+        parser = value;
+        timestep = 0;
+
+        std::cout << "[env] building the environment" << std::endl;
+
+        read_grid();
+        read_environment(timestep);
+
+        std::cout << std::endl;
+        std::cout << "[env] initialised 0th environment" << std::endl;
+
+    }
+
+    void update()
+    {
+
+        timestep = timestep + 1;
+        read_environment(timestep);
+
+        std::cout << "[env] updated to " << timestep << "th environment" << std::endl;
+
+    }
+
+private:
+
+    Parser parser;
+    int timestep;
+
+    // ------------------------------------------------------------------------
+    // generic read function for dataset
+    // ------------------------------------------------------------------------
+    void read_variable(H5::H5File fp,
+                       std::string groupname,
+                       std::string datasetname,
+                       std::vector<double> &values)
+    {
+
+        H5::Group group         = fp.openGroup(groupname.c_str());
+        H5::DataSet dataset     = group.openDataSet(datasetname.c_str());
+        H5::DataType dtype      = dataset.getDataType();
+        H5::DataSpace dataspace = dataset.getSpace();
+
+        int rank = dataspace.getSimpleExtentNdims();
+
+        hsize_t *dims = new hsize_t[rank];
+        dataspace.getSimpleExtentDims(dims);
+
+        int dim_total = 1;
+        for (int i = 0; i < rank; i ++)
+            dim_total *= dims[i];
+
+        double *data = new double[dim_total];
+        dataset.read(data, H5::PredType::IEEE_F64LE);
+
+        values.resize(dim_total);
+        for (int i = 0; i < dim_total; i ++)
+            values[i] = data[i];
+
+        delete []data;
+        delete []dims;
+
+    }
+
+    // ------------------------------------------------------------------------
+    // read in the grid
+    // ------------------------------------------------------------------------
+    void read_grid()
+    {
+
+        H5::H5File fp_grid = H5::H5File("./grid.h5", H5F_ACC_RDONLY);
+
+        // read dimensions
+        H5::DataSet dset = fp_grid.openDataSet("dimensions");
+        H5::DataType dtype = dset.getDataType();
+        H5::DataSpace dspace = dset.getSpace();
+
+        int rank = dspace.getSimpleExtentNdims();
+        hsize_t* dims = new hsize_t[rank];
+        dspace.getSimpleExtentDims(dims);
+
+        int dim_total = 1;
+        for (int i = 0; i < rank; i ++)
+            dim_total *= dims[i];
+
+        int *dimensions = new int[dim_total];
+        dset.read(dimensions, H5::PredType::NATIVE_INT);
+
+        nx = dimensions[0];
+        ny = dimensions[1];
+        nz = dimensions[2];
+
+        delete []dimensions;
+        delete []dims;
+
+        // allocate memory
+
+        indices.resize(nx);
+        elevations.resize(nx);
+        for (int i = 0; i < nx; i ++) {
+            indices[i].resize(ny);
+            elevations[i].resize(ny);
+            for (int j = 0; j < ny; j ++) {
+                indices[i][j].resize(nz);
+                elevations[i][j].resize(nz);
+            }
+        }
+
+        // read depths
+        H5::DataSet dset_depths = fp_grid.openDataSet("depths");
+        H5::DataType dtype_depths = dset_depths.getDataType();
+        H5::DataSpace dspace_depths = dset_depths.getSpace();
+
+        rank = dspace_depths.getSimpleExtentNdims();
+        dims = new hsize_t[rank];
+        dspace_depths.getSimpleExtentDims(dims);
+
+        dim_total = 1;
+        for (int i = 0; i < rank; ++i)
+            dim_total *= dims[i];
+
+        double *data_depths = new double[dim_total];
+        dset_depths.read(data_depths, H5::PredType::NATIVE_DOUBLE);
+
+        depths.resize(dim_total);
+
+        for (int i = 0; i < dim_total; ++i)
+            depths[i] = data_depths[i];
+
+        // read indices
+        H5::DataSet dset_indices = fp_grid.openDataSet("indices");
+        H5::DataType dtype_indices = dset_indices.getDataType();
+        H5::DataSpace dspace_indices = dset_indices.getSpace();
+
+        rank = dspace_indices.getSimpleExtentNdims();
+        dims = new hsize_t[rank];
+        dspace_indices.getSimpleExtentDims(dims);
+
+        dim_total = 1;
+        for (int i = 0; i < rank; i ++)
+            dim_total *= dims[i];
+
+        int *data_indices = new int[dim_total];
+        dset_indices.read(data_indices, H5::PredType::NATIVE_INT);
+
+        for (int i = 0; i < nx; i ++) {
+            for (int j = 0; j < ny; j ++) {
+                for (int k = 0; k < nz; k ++) {
+                    indices[i][j][k] = data_indices[i + nx * (j + ny * k)];
+                }
+            }
+        }
+
+        delete []data_depths;
+        delete []data_indices;
+        delete []dims;
+
+        // read elevations
+        H5::DataSet dset_elevations = fp_grid.openDataSet("elevations");
+        H5::DataType dtype_elevations = dset_elevations.getDataType();
+        H5::DataSpace dspace_elevations = dset_elevations.getSpace();
+
+        rank = dspace_elevations.getSimpleExtentNdims();
+        dims = new hsize_t[rank];
+        dspace_elevations.getSimpleExtentDims(dims);
+
+        dim_total = 1;
+        for (int i = 0; i < rank; i ++)
+            dim_total *= dims[i];
+
+        double *data_elevations = new double[dim_total];
+        dset_elevations.read(data_elevations, H5::PredType::NATIVE_DOUBLE);
+
+        for (int i = 0; i < nx; i ++) {
+            for (int j = 0; j < ny; j ++) {
+                for (int k = 0; k < nz; k ++) {
+                    elevations[i][j][k] = data_elevations[i + nx * (j + ny * k)];
+                }
+            }
+        }
+
+        delete []data_elevations;
+        delete []dims;
+
+        // read corners
+        H5::DataSet dset_corner = fp_grid.openDataSet("upperleft");
+        H5::DataType dtype_corner = dset_corner.getDataType();
+        H5::DataSpace dspace_corner = dset_corner.getSpace();
+
+        rank = dspace_corner.getSimpleExtentNdims();
+        dims = new hsize_t[rank];
+        dspace_corner.getSimpleExtentDims(dims);
+
+        dim_total = 1;
+        for (int i = 0; i < rank; i ++)
+            dim_total *= dims[i];
+
+        double *data_corner = new double[dim_total];
+        dset_corner.read(data_corner, H5::PredType::NATIVE_DOUBLE);
+
+        xll_corner = data_corner[0];
+        yhl_corner = data_corner[1];
+
+        delete []data_corner;
+        delete []dims;
+
+        fp_grid.close();
+
+    }
+
+    // ------------------------------------------------------------------------
+    // read in environment variables
+    // ------------------------------------------------------------------------
+    void read_environment(int timestep)
+    {
+
+        //
+        // we start with surface environment variables
+        //
+        H5::H5File fp_surface = H5::H5File(parser.surface_data_file.c_str(),
+                                           H5F_ACC_RDONLY);
+
+        // exchange flux
+        read_variable(fp_surface,
+                      SURFACE_EXCHANGE_FLUX,
+                      parser.surface_data_group[timestep],
+                      exchange_flux);
+        // cell area
+        // todo (this is redundant to do at each time step)
+        read_variable(fp_surface,
+                      SURFACE_AREA,
+                      parser.surface_data_group[timestep],
+                      cell_area);
+        // water depth
+        read_variable(fp_surface,
+                      SURFACE_WATER_DEPTH,
+                      parser.surface_data_group[timestep],
+                      water_depth);
+        // velocity in x-direction
+        read_variable(fp_surface,
+                      SURFACE_VELOCITY_X,
+                      parser.surface_data_group[timestep],
+                      velocity_x);
+        // velocity in y-direction
+        read_variable(fp_surface,
+                      SURFACE_VELOCITY_Y,
+                      parser.surface_data_group[timestep],
+                      velocity_y);
+        // todo (use transpiration)
+        read_variable(fp_surface,
+                      SURFACE_EVAPOTRANSPIRATION_FLUX,
+                      parser.surface_data_group[timestep],
+                      evaporation_flux);
+
+        fp_surface.close();
+
+        //
+        // now we read in subsurface
+        //
+        H5::H5File fp_subsurface = H5::H5File(parser.subsurface_data_file.c_str(),
+                                              H5F_ACC_RDONLY);
+
+        // darcy velocity in x-direction
+        read_variable(fp_subsurface,
+                      DARCY_VELOCITY_X,
+                      parser.subsurface_data_group[timestep],
+                      darcy_velocity_x);
+
+        // darcy velocity in y-direction
+        read_variable(fp_subsurface,
+                      DARCY_VELOCITY_Y,
+                      parser.subsurface_data_group[timestep],
+                      darcy_velocity_y);
+
+        // darcy velocity in z-direction
+        read_variable(fp_subsurface,
+                      DARCY_VELOCITY_Z,
+                      parser.subsurface_data_group[timestep],
+                      darcy_velocity_z);
+
+        fp_subsurface.close();
+
+    }
+
+};
+
+#endif
diff --git a/header.hpp b/header.hpp
index 9976419..0575269 100644
--- a/header.hpp
+++ b/header.hpp
@@ -5,29 +5,21 @@
 
 #include <iostream>
 #include <vector>
+#include <array>
 #include <string>    // string class
 #include <algorithm> // min
 #include <float.h>   // DBL_MAX
-#include <math.h>    // fabs, sqrt
-#include <mpi.h>
+#include <math.h>    // fabs, sqrt, ceil
 #include <fstream>   // ifstream
 #include <sstream>   // stringstream
 #include <string.h>  // strcmp
+#include <random>    // random number generator
+#include <omp.h>
 
 #include <H5Cpp.h>
 #include <H5File.h>
 
-#define CLASSIC_SCHEME 1
-
-#define MAXTIMESTEP 864000.0 // (s) ten days
-#define EPS 1.0e-10
-
-#define CFL 10.0
-
 #define DIM 3
-#define N_NEIGHBORS 5
-#define N_VERTICES 6
-#define N_EDGES 9
 
 #define SURFACE_NCOL 4
 #define SUBSURFACE_NCOL 7
@@ -35,16 +27,35 @@
 #define INFILTRATE_TRUE 1
 #define INFILTRATE_FALSE 2
 
-// #define SURFACE_DATA_GROUP_NAME_H5 "11841"
-// #define SUBSURFACE_DATA_GROUP_NAME_H5 "11841"
-
 #define MESH_MIXED_ELEMENTS "/0/Mesh/MixedElements"
 #define MESH_NODES "/0/Mesh/Nodes"
 
-#define AREADIV 0.01
+#define SURFACE_EXCHANGE_FLUX "surface-surface_subsurface_flux.cell.0"
+#define SURFACE_WATER_DEPTH "surface-ponded_depth.cell.0"
+#define SURFACE_VELOCITY_X "surface-velocity.cell.0"
+#define SURFACE_VELOCITY_Y "surface-velocity.cell.1"
+#define SURFACE_EVAPOTRANSPIRATION_FLUX "surface-potential_transpiration.cell.0"
+#define SURFACE_AREA "surface-cell_volume.cell.0"
+
+#define DARCY_VELOCITY_X "darcy_velocity.cell.0"
+#define DARCY_VELOCITY_Y "darcy_velocity.cell.1"
+#define DARCY_VELOCITY_Z "darcy_velocity.cell.2"
 
 #define SUBSURFACE_FLAG 1
 #define SURFACE_FLAG 2
-#define OUTSIDE_FLAG 3
+#define OUTLET_FLAG 3
+#define EVAPORATION_FLAG 4
+
+#define WATERDEPTH_DRY 1.0e-8
+#define ZERO_VELOCITY 1.0e-10
+#define DT_MAX 432000.0
+
+#define MOLARVOLUME 1.8e-5
+
+#define CFL 1.0
+
+typedef std::array<double, DIM> point;
+typedef std::vector<std::vector<std::vector<int>>> index3d;
+typedef std::vector<std::vector<std::vector<double>>> grid3d;
 
 #endif
diff --git a/input/groups.txt b/input/groups.txt
new file mode 100644
index 0000000..beb3810
--- /dev/null
+++ b/input/groups.txt
@@ -0,0 +1,108 @@
+0
+51
+68
+88
+106
+116
+123
+128
+133
+159
+191
+215
+228
+232
+247
+253
+257
+265
+273
+279
+283
+290
+294
+298
+302
+306
+313
+319
+323
+327
+331
+334
+338
+344
+348
+352
+356
+359
+363
+368
+374
+379
+387
+391
+395
+398
+400
+402
+406
+410
+415
+418
+422
+425
+429
+437
+444
+448
+452
+456
+461
+463
+466
+470
+474
+477
+484
+489
+493
+497
+503
+508
+513
+520
+524
+528
+532
+537
+547
+573
+588
+597
+613
+627
+637
+647
+656
+672
+683
+690
+694
+701
+706
+716
+733
+754
+787
+797
+814
+836
+852
+857
+873
+899
+934
+958
+990
+998
diff --git a/input/options.input b/input/options.input
index ab297f6..f48010d 100644
--- a/input/options.input
+++ b/input/options.input
@@ -1,16 +1,18 @@
 // input file for agent provocateur
 
-n_agent         : 100						// Number of agents
-fuzz            : 0.1						// Fuzzy equivalence constant, used for geometry
+n_agent         : 1000						// Number of agents
+fuzz            : 0.01						// Fuzzy equivalence constant, used for geometry
 seed            : 108						// Random seed
 n_iteration     : 1000						// Number of iteration steps
+subgrid_dx      : 5.0						// Subgrid cell size
 
-surfacemesh	: ./001-run/simrun-000004/ats_vis_surface_mesh.h5	// Surface mesh file
-surfacedata	: ./001-run/simrun-000004/ats_vis_surface_data.h5	// Surface data file
-surfacetime	: 0.0 172800.0 432000.0 518400.0 691200.0		// Time levels
-surfacegroup	: 0 106 116 123 128					// Group names for surface
+surfacemesh	: ./input/test/ats_vis_surface_mesh.h5		// Surface mesh file
+surfacedata	: ./input/test/ats_vis_surface_data.h5		// Surface data file
+surfacetime     : ./input/times.txt				// Time levels for surface 
+surfacegroup	: ./input/groups.txt				// Group names for surface
+
+subsurfacemesh  : ./input/test/ats_vis_mesh.h5			// Subsurface mesh file
+subsurfacedata  : ./input/test/ats_vis_data.h5			// Subsurface data file
+subsurfacetime  : ./input/subtimes.txt			 	// Time levels for subsurface
+subsurfacegroup : ./input/subgroups.txt				// Group names for subsurface
 
-subsurfacemesh  : ./001-run/simrun-000004/ats_vis_mesh.h5		// Subsurface mesh file
-subsurfacedata  : ./001-run/simrun-000004/ats_vis_data.h5		// Subsurface data file
-subsurfacetime  : 0.0 432000.0 864000.0 1296000.0 1728000.0		// Time levels
-subsurfacegroup : 0 116 130 191	230					// Group names for subsurface
diff --git a/input/subgroups.txt b/input/subgroups.txt
new file mode 100644
index 0000000..a5f28cb
--- /dev/null
+++ b/input/subgroups.txt
@@ -0,0 +1,43 @@
+0
+77
+116
+130
+191
+230
+253
+269
+283
+296
+306
+321
+331
+340
+352
+361
+374
+389
+398
+404
+415
+424
+437
+450
+461
+468
+477
+491
+503
+516
+528
+541
+588
+623
+647
+679
+694
+709
+754
+802
+852
+887
+958
diff --git a/input/subtimes.txt b/input/subtimes.txt
new file mode 100644
index 0000000..45adcab
--- /dev/null
+++ b/input/subtimes.txt
@@ -0,0 +1,43 @@
+0
+118276
+236550
+354826
+473100
+591376
+709651
+827926
+946201
+1.06447e+06
+1.18276e+06
+1.30103e+06
+1.4193e+06
+1.53757e+06
+1655856
+1.77413e+06
+1.8924e+06
+2.01067e+06
+2.12896e+06
+2.24723e+06
+2.3655e+06
+2.48378e+06
+2.60206e+06
+2.72033e+06
+2.8386e+06
+2.95688e+06
+3.07516e+06
+3.19343e+06
+3.3117e+06
+3.42998e+06
+3.54826e+06
+3.66653e+06
+3.7848e+06
+3.90308e+06
+4.02136e+06
+4.13963e+06
+4.2579e+06
+4.37618e+06
+4.49446e+06
+4.61273e+06
+4.731e+06
+4.84928e+06
+4.96756e+06
diff --git a/input/times.txt b/input/times.txt
new file mode 100644
index 0000000..45adcab
--- /dev/null
+++ b/input/times.txt
@@ -0,0 +1,43 @@
+0
+118276
+236550
+354826
+473100
+591376
+709651
+827926
+946201
+1.06447e+06
+1.18276e+06
+1.30103e+06
+1.4193e+06
+1.53757e+06
+1655856
+1.77413e+06
+1.8924e+06
+2.01067e+06
+2.12896e+06
+2.24723e+06
+2.3655e+06
+2.48378e+06
+2.60206e+06
+2.72033e+06
+2.8386e+06
+2.95688e+06
+3.07516e+06
+3.19343e+06
+3.3117e+06
+3.42998e+06
+3.54826e+06
+3.66653e+06
+3.7848e+06
+3.90308e+06
+4.02136e+06
+4.13963e+06
+4.2579e+06
+4.37618e+06
+4.49446e+06
+4.61273e+06
+4.731e+06
+4.84928e+06
+4.96756e+06
diff --git a/main.cpp b/main.cpp
index 62d491a..45d3cb1 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,284 +1,75 @@
 // -*- mode: c++ -*-
 
-#include "mesh.hpp"
-#include "agent.hpp"
 #include "parser.hpp"
+#include "mesh.hpp"
+#include "environment.hpp"
+#include "particle.hpp"
 
 int main(int argc, char* argv[])
 {
 
-  int rank;
-  int size;
-
-  MPI_Init(&argc, &argv);
-  MPI_Comm_size(MPI_COMM_WORLD, &size);
-  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-
-  if (rank == 0) {
-    
-    std::cout << std::endl;
-    std::cout << "       /  /^    agent provocateur" << std::endl;
-    std::cout << "      o--o      =================" << std::endl;
-    std::cout << std::endl;
-
-    std::cout << "[~>] Read in input data" << std::endl;
-
-  }
-
-  Parser parser;
-  parser.read_options("input/options.input");
+    if (argc < 3) {
+        std::cerr << "[err] usage: " << argv[0] << " -[c|a] <fname>.input" << std::endl;
+        std::cerr << "[err] c -> compile mesh | a -> advect particles" << std::endl;
+        return EXIT_FAILURE;
+    }
 
-  if (rank == 0) {
+    // ------------------------------------------------------------------------
+    // handle command line arguments
+    // ------------------------------------------------------------------------
 
-    std::cout << "[~~>] Number of agents: " << parser.n_agent << std::endl;
-    std::cout << "[~~>] Number of iterations: " << parser.n_iteration << std::endl;
-    std::cout << "[~~>] Random seed: " << parser.seed << std::endl;
-    std::cout << "[~~>] Geometric fuzz: " << parser.fuzz << std::endl;
+    std::string flag = argv[1];
 
-    std::cout << "[~~>] Surface mesh file: " << parser.surface_mesh_file << std::endl;
-    std::cout << "[~~>] Surface data file: " << parser.surface_data_file << std::endl;
-    std::cout << "[~~>] Surface time: ";
-    for (uint i = 0; i < parser.surface_time.size(); i ++)
-      std::cout << parser.surface_time[i] << " ";
-    std::cout <<std::endl;
-    std::cout << "[~~>] Surface group: ";
-    for (uint i = 0; i < parser.surface_data_group.size(); i ++)
-      std::cout << parser.surface_data_group[i] << " ";
-    std::cout <<std::endl;
+    if (flag != "-c" && flag != "-a") {
+        std::cerr << "[err] flag not understood: options are -c and -a" << std::endl;
+        return EXIT_FAILURE;
+    }
 
-    std::cout << "[~~>] Subsurface mesh file: " << parser.subsurface_mesh_file << std::endl;
-    std::cout << "[~~>] Subsurface data file: " << parser.subsurface_data_file << std::endl;
-    std::cout << "[~~>] Subsurface time: ";
-    for (uint i = 0; i < parser.subsurface_time.size(); i ++)
-      std::cout << parser.subsurface_time[i] << " ";
-    std::cout <<std::endl;
-    std::cout << "[~~>] Subsurface group: ";
-    for (uint i = 0; i < parser.subsurface_data_group.size(); i ++)
-      std::cout << parser.subsurface_data_group[i] << " ";
-    std::cout << std::endl;
+    // ------------------------------------------------------------------------
+    // read configuration file
+    // ------------------------------------------------------------------------
     std::cout << std::endl;
-    
-    std::cout << "[~>] Building the subsurface mesh" << std::endl;
-    
-  }
-
-  Mesh mesh = Mesh();
-
-  mesh.set_fuzz(parser.fuzz);
-  
-  H5::H5File fp_subsurface = H5::H5File(parser.subsurface_mesh_file.c_str(), H5F_ACC_RDONLY);
-  H5::H5File dp_subsurface = H5::H5File(parser.subsurface_data_file.c_str(), H5F_ACC_RDONLY);
-  H5::H5File fp_surface = H5::H5File(parser.surface_mesh_file.c_str(), H5F_ACC_RDONLY);
-  H5::H5File dp_surface = H5::H5File(parser.surface_data_file.c_str(), H5F_ACC_RDONLY);
-  
-  mesh.build_subsurface_mesh(fp_subsurface);
-  fp_subsurface.close();
-
-  if (rank == 0) {
-    std::cout << "[~~>] Number of cells: " << mesh.get_n_cells(SUBSURFACE_FLAG) << std::endl;
+    std::cout << "      /|  |\\   tracksuit" << std::endl;
+    std::cout << "       //\\\\    =========" << std::endl;
     std::cout << std::endl;
-    
-    std::cout << "[~>] Building the surface mesh (includes coupling with subsurface and may take a while)" << std::endl;
-  }
-  
-  mesh.build_surface_mesh(fp_surface, dp_surface, parser.surface_data_group[0]);
-  fp_surface.close();
 
-  if (rank == 0) {
-
-    int len_surf = mesh.get_n_cells(SURFACE_FLAG);
-
-    std::cout << "[~~>] Number of cells: " << len_surf << std::endl;
-    std::cout << std::endl;
-    
-    FILE *mesh_fp = fopen("mesh.csv", "w");
-    fprintf(mesh_fp, "x (m), y (m), z (m), index\n");
-    Cell *c;
-    
-    for (int n = 0; n < len_surf; n ++) {
-      c = mesh.get_cell(n, SURFACE_FLAG);
-      fprintf(mesh_fp, "%f, %f, %f, %d\n",
-	      c->get_coordinate(0),
-	      c->get_coordinate(1),
-	      c->get_coordinate(2), n);
+    Parser parser;
+    parser.read_options(argv[2]);
+    parser.print_options();
+
+    // ------------------------------------------------------------------------
+    // compile the mesh to produce a subgrid
+    // ------------------------------------------------------------------------
+
+    if (flag == "-c") {
+        Mesh mesh = Mesh(parser);
+        mesh.compile();
+        // mesh.subgrid();
+        // mesh.map();
+        mesh.dump();
+        std::cout << "[   ] meshing completed, enjoy the subgrid" << std::endl;
+        std::cout << std::endl;
     }
-  
-    fclose(mesh_fp);
-
-    std::cout << "[~>] Reading the subsurface environment" << std::endl;
-  
-  }
 
-  // Place agents in the domain and advance N_ITERATION steps
+    // ------------------------------------------------------------------------
+    // start the particle tracking
+    // ------------------------------------------------------------------------
 
-  // Bootstrap 
-  srand(parser.seed + rank);
-  int len = mesh.get_n_cells(SUBSURFACE_FLAG);
-  int index;
+    if (flag == "-a") {
 
-  std::cout << "[~>] Releasing " << parser.n_agent << " agents" << std::endl;
+        ParticleManager pmanager = ParticleManager(parser);
 
-  size_t partition = parser.n_agent / size;
-  size_t low_bound = rank * partition;
-  size_t high_bound = (rank + 1) * partition;
-  if (rank == size - 1) // top processor picks up extras
-    high_bound = parser.n_agent;
-  
-  for (size_t i = low_bound; i < high_bound; i ++) {
-
-    // for each agent, start with a fresh environment
-    H5::H5File _dp_subsurface = H5::H5File(parser.subsurface_data_file.c_str(), H5F_ACC_RDONLY);
-	  
-    mesh.build_environment(_dp_subsurface,
-			   parser.subsurface_data_group[0],
-			   parser.subsurface_data_group[1],
-			   SUBSURFACE_FLAG);
-    _dp_subsurface.close();
-
-    H5::H5File _dp_surface = H5::H5File(parser.surface_data_file.c_str(), H5F_ACC_RDONLY);
-	  
-    mesh.build_environment(_dp_surface,
-			   parser.surface_data_group[0],
-			   parser.surface_data_group[1],
-			   SURFACE_FLAG);
-    _dp_surface.close();
-
-    // create agent
-
-    Agent agent;
-    agent.init(parser.n_iteration);
-    
-    index = rand()%len;
-    agent.set_cell_index(index);
-    agent.set_path(mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(0),
-		   mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(1),
-		   mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(2), 0);
-  
-    agent.set_domain(SUBSURFACE_FLAG, 0);
-    agent.set_domain(SUBSURFACE_FLAG, 1);
-
-    agent.set_path(mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(0),
-		   mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(1),
-		   mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(2), 1);
-
-    int step;
-    int surface_time_counter = 0;
-    int subsurface_time_counter = 0;
-    int n_levels_subsurface = parser.subsurface_time.size();
-    int n_levels_surface    = parser.surface_time.size();
-
-    // if the size of the vectors is larger than 1, there are multiple
-    // time levels and the environment is dynamic.  we increase the
-    // counter to the next time step to linearly interpolate between
-    // values
-    if (n_levels_subsurface > 1)
-      subsurface_time_counter ++;
-    if (n_levels_surface > 1)
-      surface_time_counter ++;
-    
-    double time_horizon_surface    = parser.surface_time[surface_time_counter];
-    double time_horizon_subsurface = parser.subsurface_time[subsurface_time_counter];
-    
-    char output_fname[1024];
-
-    // move agent
-    for (step = 1; step < parser.n_iteration; step ++) {
-
-      double rel_dt_surface = (agent.get_time(step) - parser.surface_time[surface_time_counter-1]) / (parser.surface_time[surface_time_counter] - parser.surface_time[surface_time_counter-1]);
-      double rel_dt_subsurface = (agent.get_time(step) - parser.subsurface_time[subsurface_time_counter-1]) / (parser.subsurface_time[subsurface_time_counter] - parser.subsurface_time[subsurface_time_counter-1]);
-
-      rel_dt_surface    = std::max(0.0, rel_dt_surface);
-      rel_dt_subsurface = std::max(0.0, rel_dt_subsurface);
-      
-      agent.advance(&mesh, rel_dt_surface, rel_dt_subsurface, step);
-
-      // update environment if necessary
-      double current_time = agent.get_time(step);
-      double dt_max       = time_horizon_subsurface - current_time;
-      agent.set_timestepsize(dt_max);
-
-      if ((step - 1) % (parser.n_iteration / 5) == 0 && i % (parser.n_agent / 5) == 0) {
-	std::cout << std::endl;
-	std::cout << "Agent #" << i << std::endl;
-	std::cout << "Current time: " << current_time << "/" << time_horizon_subsurface << " s";
-	std::cout << std::endl;
-	std::cout << "Current counter: " << subsurface_time_counter << "/" << n_levels_subsurface - 1 << std::endl;
-	std::cout << std::endl;
-      }
-
-      // check the subsurface
-      if (subsurface_time_counter < n_levels_subsurface - 1) {
-	
-	// once the simulation time exceeds the time frame of the
-	// environment, we assume steady state
-
-	if (current_time >= time_horizon_subsurface) {
-
-	  if (step % (parser.n_iteration / 5) == 0 && i % (parser.n_agent / 5) == 0) {
-	    std::cout << "  Advance to group " << parser.subsurface_data_group[subsurface_time_counter] << std::endl;
-	    std::cout << std::endl;
-	  }
-	  
-	  time_horizon_subsurface = parser.subsurface_time[subsurface_time_counter+1];
-	  
-	  H5::H5File _dp = H5::H5File(parser.subsurface_data_file.c_str(), H5F_ACC_RDONLY);
-	  mesh.build_environment(_dp,
-				 parser.subsurface_data_group[subsurface_time_counter],
-				 parser.subsurface_data_group[subsurface_time_counter+1],
-				 SUBSURFACE_FLAG);
-	  _dp.close();
-	  
-	  subsurface_time_counter ++;
-	  
-	}
-	
-      }
-
-      if (surface_time_counter < n_levels_surface - 1) {
-
-	if (current_time >= time_horizon_surface) {
-
-	  if (step % (parser.n_iteration / 5) == 0 && i % (parser.n_agent / 5) == 0) {
-	    std::cout << "  Advance to group " << parser.surface_data_group[surface_time_counter] << std::endl;
-	    std::cout << std::endl;
-	  }
-	  
-	  time_horizon_surface = parser.surface_time[surface_time_counter+1];
-	  
-	  H5::H5File _dp = H5::H5File(parser.surface_data_file.c_str(), H5F_ACC_RDONLY);
-	  mesh.build_environment(_dp,
-				 parser.surface_data_group[surface_time_counter],
-				 parser.surface_data_group[surface_time_counter+1],
-				 SURFACE_FLAG);
-	  _dp.close();
-	  
-	  surface_time_counter ++;
-	  
-	}
-	
-      }
-      
-    }
+        for (size_t t = 1; t < parser.surface_time.size(); ++t) {
+            pmanager.advect(parser.surface_time[t]);
+            pmanager.environment.update();
+        }
 
-    // write out results
-    sprintf(output_fname, "results/agent-%06lu.csv", i);
-    FILE *result_fp = fopen(output_fname, "w");
+        for (size_t i = 0; i < pmanager.particles.size(); ++i) {
+            pmanager.particles[i].dump(i);
+        }
 
-    fprintf(result_fp, "time (s), x (m), y (m), z (m), domain\n");
-    for (step = 0; step < parser.n_iteration; step ++) {
-      fprintf(result_fp, "%f, %f, %f, %f, %d\n",
-	      agent.get_time(step),
-	      agent.get_path(step, 0),
-	      agent.get_path(step, 1),
-	      agent.get_path(step, 2),
-	      agent.get_domain(step));
     }
-    fclose(result_fp);
-      
-  }
 
-  MPI_Finalize();
-  
-  return MPI_SUCCESS;
+    return EXIT_SUCCESS;
 
 }
diff --git a/makefile b/makefile
index ea6e6aa..c0ff9e3 100644
--- a/makefile
+++ b/makefile
@@ -1,12 +1,12 @@
-PRG=agent-provocateur
+PRG=ptracker
 
-#CC=g++-11
-CC=mpic++
-CXXFLAGS=-Wall -Wextra -pedantic -O0
+CC=g++-11
+#CC=mpic++
+CXXFLAGS=-Wall -Wextra -pedantic -O0 -fopenmp
+#LFLAGS=-I/home/iozgen/Documents/software/hdf5/include -L/home/iozgen/Documents/software/hdf5/lib
 LFLAGS=-I/opt/hdf5/include -L/opt/hdf5/lib
-SRC=parser.cpp cell.cpp mesh.cpp main.cpp agent.cpp
+SRC=main.cpp
 LIB=-lm -lhdf5 -lhdf5_cpp
-
 OBJ=$(SRC:.cpp=.o)
 
 all:$(PRG)
@@ -20,14 +20,20 @@ $(PRG): $(OBJ)
 	$(CC) $(LFLAGS) $(CXXFLAGS) -c $< -o $@
 
 run: $(PRG)
-	./$(PRG)
+	./$(PRG) input/options.input
+
+tag: $(SRC)
+	etags *.hpp *.cpp
+
+indent: $(SRC)
+	astyle --style=kr *.hpp *.cpp
 
 clean:
-	rm -f $(PRG) $(OBJ)
+	rm -f $(OBJ) *.orig
 
 cleanall:
 	make clean
-	rm -f *.dat; rm -f *.png
+	rm -f *.dat; rm -f *.png *.csv TAGS $(PRG)
 
-.PHONY: all clean cleanall
+.PHONY: all tag clean cleanall
 
diff --git a/mesh.cpp b/mesh.cpp
deleted file mode 100644
index 406dfba..0000000
--- a/mesh.cpp
+++ /dev/null
@@ -1,544 +0,0 @@
-// -*- mode: c++ -*-
-
-#include "mesh.hpp"
-
-
-// ------------------------------------------------------------------------
-// -- PUBLIC METHODS
-// ------------------------------------------------------------------------
-
-void Mesh::build_surface_mesh(H5::H5File fp, H5::H5File ep, std::string datasetname)
-{
-
-  get_mesh_metadata(fp, SURFACE_NCOL);
-
-  handle_surface_node_data(fp);
-
-  build_surface_cells(ep, datasetname);
-  
-}
-
-void Mesh::build_subsurface_mesh(H5::H5File fp)
-{
-
-  get_mesh_metadata(fp, SUBSURFACE_NCOL);
-
-  handle_subsurface_node_data(fp);
-
-  build_subsurface_cells();
-  
-}
-
-void Mesh::build_environment(H5::H5File ep,
-			     std::string datasetname,
-			     std::string datasetname_next,
-			     int flag)
-{
-
-  std::string groupname;
-  
-  // Read in velocities in x direction
-    
-  if (flag == SUBSURFACE_FLAG) {
-    groupname = "darcy_velocity.cell.0";
-  } else if (flag == SURFACE_FLAG) {
-    groupname = "surface-velocity.cell.0";
-  }
-
-  read_velocity(ep, groupname, datasetname, datasetname_next, flag, 0);
-
-  // Read in velocities in y direction
-
-  if (flag == SUBSURFACE_FLAG) {
-    groupname = "darcy_velocity.cell.1";
-  } else if (flag == SURFACE_FLAG) {
-    groupname = "surface-velocity.cell.1";
-  }
-
-  read_velocity(ep, groupname, datasetname, datasetname_next, flag, 1);
-
-  // Read velocities in z direction
-  
-  if (flag == SUBSURFACE_FLAG) {
-    groupname = "darcy_velocity.cell.2";
-  } else if (flag == SURFACE_FLAG) {
-    groupname = "surface-surface_subsurface_flux.cell.0";
-  }
-
-  read_velocity(ep, groupname, datasetname, datasetname_next, flag, 2);
-
-  
-  // Surface environment needs adjustments to the vertical velocity
-  // and needs to read in water depth
-  
-  if (flag == SURFACE_FLAG) {
-    
-    // Rescale the exchange flux from mole to m/s
-    for (int i = 0; i < n_cells_surface; i ++) {
-      surface_cells[i].set_velocity(2, surface_cells[i].get_velocity(2) * 1.8e-5 * AREADIV);
-      surface_cells[i].set_velocity_next(2, surface_cells[i].get_velocity_next(2) * 1.8e-5 * AREADIV);
-    }
-
-    // Read in water depth
-    read_water_depth(ep, datasetname, datasetname_next);
-    
-  }
-  
-}
-
-// ------------------------------------------------------------------------
-// GETTERS AND SETTERS
-// ------------------------------------------------------------------------
-
-Mesh::Mesh()
-{
-  n_cells_surface = 0;
-  n_cells_subsurface = 0;
-  high_coordinates.resize(DIM);
-  low_coordinates.resize(DIM);
-}
-
-int Mesh::get_n_cells(int flag)
-{
-  int n_cells;
-  
-  if (flag == SUBSURFACE_FLAG) {
-    n_cells = n_cells_subsurface;
-  } else {
-    n_cells = n_cells_surface;
-  }
-  
-  return n_cells;
-}
-
-Cell* Mesh::get_cell(int index, int flag)
-{
-
-  if (flag == SUBSURFACE_FLAG) {
-    return &(subsurface_cells[index]);
-  } else {
-    return &(surface_cells[index]);
-  } 
-}
-
-void Mesh::set_high_coordinate(int index, double value)
-{
-  high_coordinates[index] = value;
-}
-
-double Mesh::get_high_coordinate(int index)
-{
-  return high_coordinates[index];
-}
-
-void Mesh::set_low_coordinate(int index, double value)
-{
-  low_coordinates[index] = value;
-}
-
-double Mesh::get_low_coordinate(int index)
-{
-  return low_coordinates[index];
-}
-
-void Mesh::set_fuzz(double value)
-{
-  fuzz = value;
-}
-
-// ------------------------------------------------------------------------
-// -- PRIVATE METHODS
-// ------------------------------------------------------------------------
-
-void Mesh::handle_surface_node_data(H5::H5File fp)
-{
-
-  H5::DataSet dataset_n = fp.openDataSet(MESH_NODES);
-  H5::DataType datatype_n = dataset_n.getDataType();
-  H5::DataSpace dataspace_n = dataset_n.getSpace();
-
-  int rank = dataspace_n.getSimpleExtentNdims();
-
-  hsize_t *dims_n = new hsize_t[rank];
-  dataspace_n.getSimpleExtentDims(dims_n);
-
-#if VERBOSE
-  for (int i = 0; i < rank; i ++)
-    std::cout << "[~~~>] Dimension in " << i << " of nodes: " << dims_n[i] << std::endl;
-#endif
-  
-  double *data_n = new double[dims_n[0]*dims_n[1]];
-  dataset_n.read(data_n, H5::PredType::IEEE_F64LE);
-  
-  nodes.resize(dims_n[0]);
-
-  for (size_t i = 0; i < dims_n[0]; i ++) {
-    nodes[i].resize(dims_n[1]);
-    for (size_t j = 0; j < dims_n[1]; j ++)
-      nodes[i][j] = data_n[i * dims_n[1] + j];
-  }
-  
-}
-
-void Mesh::build_surface_cells(H5::H5File ep, std::string datasetname)
-{
-
-  std::vector<double> coords1(DIM-1);
-  std::vector<double> coords2(DIM-1);
-  std::vector<double> coords3(DIM-1);
-  std::vector<double> coords_avg(DIM-1);
-
-  read_elevation(ep, datasetname);
-
-  for (int i = 0; i < n_cells_surface; i ++) {
-
-    int n1 = element_node_indices[i][1];
-    int n2 = element_node_indices[i][2];
-    int n3 = element_node_indices[i][3];
-
-    std::vector<double> edge_length(DIM);
-
-    for (size_t k = 0; k < DIM; k ++)
-      edge_length[k] = DBL_MAX;
-
-
-    double smallest_distance;
-    double largest_distance;
-    
-    for (size_t j = 0; j < DIM-1; j ++) {
-
-      coords1[j] = nodes[n1][j];
-      coords2[j] = nodes[n2][j];
-      coords3[j] = nodes[n3][j];
-
-      coords_avg[j] = (coords1[j] + coords2[j] + coords3[j]) / 3.0;
-
-      edge_length[0] = fabs(coords1[j] - coords_avg[j]);
-      edge_length[1] = fabs(coords2[j] - coords_avg[j]);
-      edge_length[2] = fabs(coords3[j] - coords_avg[j]);
-
-      smallest_distance = edge_length[0];
-      largest_distance  = edge_length[0];
-
-      for (size_t k = 1; k < DIM; k ++) {
-	smallest_distance = std::min(smallest_distance, edge_length[k]);
-	largest_distance  = std::max(largest_distance, edge_length[k]);
-      }
-      
-    }
-
-    surface_cells[i].set_characteristic_length(smallest_distance);
-    surface_cells[i].set_maximum_length(largest_distance);
-
-    for (size_t k = 0; k < DIM-1; k ++)
-      surface_cells[i].set_coordinate(k, coords_avg[k]);
-    
-    surface_cells[i].set_cell_id((int) i);
-
-    // find the coupled cell in the subsurface
-    std::vector<double> distance(DIM);
-    int coupled_neighbor = -1; // store current coupled neighbor here
-    double dz = DBL_MAX;       
-
-    for (int k = 0; k < n_cells_subsurface; k ++) {
-
-      for (size_t j = 0; j < DIM; j ++)
-	distance[j] = subsurface_cells[k].get_coordinate(j) - surface_cells[i].get_coordinate(j);
-
-      // Because the mesh is a stacked mesh, the only candidates for the subsurface
-      // neighbor are the cells located in the same column.  These cells will have
-      // the same centroid in the xy-plane as the surface cell.  The problem reduces
-      // to finding the nearest cell in the vertical direction in the column.
-      
-      double xy_distance = sqrt(distance[0]*distance[0] + distance[1]*distance[1]);
-
-      if (xy_distance <= fuzz) {
-	dz = std::min(dz, fabs(distance[2]));
-	coupled_neighbor = k;
-	if (dz <= 2.0 * subsurface_cells[k].get_characteristic_length())
-	  break;
-      }
-      
-    }
-
-    surface_cells[i].set_coupled_neighbor(coupled_neighbor);
-    subsurface_cells[coupled_neighbor].set_coupled_neighbor(i);
-    
-  }
-
-}
-
-void Mesh::handle_subsurface_node_data(H5::H5File fp)
-{
-  
-  H5::DataSet dataset_n = fp.openDataSet(MESH_NODES);
-  H5::DataType datatype_n = dataset_n.getDataType();
-  H5::DataSpace dataspace_n = dataset_n.getSpace();
-
-  int rank = dataspace_n.getSimpleExtentNdims();
-
-  hsize_t *dims_n = new hsize_t[rank];
-  dataspace_n.getSimpleExtentDims(dims_n);
-
-#if VERBOSE
-  for (int i = 0; i < rank; i ++)
-    std::cout << "[~~~>] Dimension in " << i << " of nodes: " << dims_n[i] << std::endl;
-#endif
-  
-  double *data_n = new double[dims_n[0]*dims_n[1]];
-  dataset_n.read(data_n, H5::PredType::IEEE_F64LE);
-
-  nodes.resize(dims_n[0]);
-
-  for (size_t i = 0; i < dims_n[0]; i ++) {
-    nodes[i].resize(dims_n[1]);
-    for (size_t j = 0; j < dims_n[1]; j ++)
-      nodes[i][j] = data_n[i * dims_n[1] + j];
-  }
-  
-}
-
-void Mesh::build_subsurface_cells()
-{
-
-  std::vector<double> coords1(DIM);
-  std::vector<double> coords2(DIM);
-  std::vector<double> coords3(DIM);
-  std::vector<double> coords4(DIM);
-  std::vector<double> coords5(DIM);
-  std::vector<double> coords6(DIM);
-
-  std::vector<double> coords_avg(DIM);
-
-  for (size_t i = 0; i < DIM; i ++) {
-    set_high_coordinate(i, DBL_MIN);
-    set_low_coordinate(i, DBL_MAX);
-  }
-
-  double smallest_distance;
-  double largest_distance;
-    
-  for (int i = 0; i < n_cells_subsurface; i ++) {
-
-    int n1 = element_node_indices[i][1];
-    int n2 = element_node_indices[i][2];
-    int n3 = element_node_indices[i][3];
-    int n4 = element_node_indices[i][4];
-    int n5 = element_node_indices[i][5];
-    int n6 = element_node_indices[i][6];
-
-    std::vector<double> edge_length(N_VERTICES);
-
-    for (size_t k = 0; k < N_VERTICES; k ++)
-      edge_length[k] = DBL_MAX;
-    
-    for (size_t j = 0; j < DIM; j++) {
-      
-      coords1[j] = nodes[n1][j];
-      coords2[j] = nodes[n2][j];
-      coords3[j] = nodes[n3][j];
-      coords4[j] = nodes[n4][j];
-      coords5[j] = nodes[n5][j];
-      coords6[j] = nodes[n6][j];
-
-      coords_avg[j] = (coords1[j] + coords2[j] + coords3[j] + coords4[j] + coords5[j] + coords6[j]) / 6.0;
-
-      // the characteristic length is the radius of the insphere of
-      // the wedge
-      edge_length[0] = fabs(coords1[j] - coords_avg[j]);
-      edge_length[1] = fabs(coords2[j] - coords_avg[j]);
-      edge_length[2] = fabs(coords3[j] - coords_avg[j]);
-      edge_length[3] = fabs(coords4[j] - coords_avg[j]);
-      edge_length[4] = fabs(coords5[j] - coords_avg[j]);
-      edge_length[5] = fabs(coords6[j] - coords_avg[j]);
-
-      smallest_distance = edge_length[0];
-      largest_distance  = edge_length[0];
-
-      for (size_t k = 1; k < N_VERTICES; k ++) {
-      smallest_distance = std::min(smallest_distance, edge_length[k]);
-      largest_distance  = std::max(largest_distance, edge_length[k]);
-      }
-      
-    }
-
-    subsurface_cells[i].set_characteristic_length(smallest_distance);
-    subsurface_cells[i].set_maximum_length(largest_distance);
-    
-    for (size_t k = 0; k < DIM; k ++) {
-
-      // update lowest and highest corner of mesh
-      high_coordinates[k] = std::max(coords_avg[k], high_coordinates[k]);
-      low_coordinates[k] = std::min(coords_avg[k], low_coordinates[k]);
-
-      subsurface_cells[i].set_coordinate(k, coords_avg[k]);
-      
-    }
-    
-    subsurface_cells[i].set_cell_id((int) i);
-    
-  }
-  
-}
-
-void Mesh::get_mesh_metadata(H5::H5File fp, int flag)
-{
-
-  H5::DataSet dataset_e = fp.openDataSet(MESH_MIXED_ELEMENTS);
-  H5::DataType datatype_e = dataset_e.getDataType();
-  H5::DataSpace dataspace_e = dataset_e.getSpace();
-
-  int rank = dataspace_e.getSimpleExtentNdims();
-
-  hsize_t *dims_e = new hsize_t[rank];
-  dataspace_e.getSimpleExtentDims(dims_e);
-
-#if VERBOSE
-  for (int i = 0; i < rank; i ++)
-      std::cout << "[~~~>] Dimension in " << i << " of element node map: " << dims_e[i] << std::endl;
-#endif
-
-  int *data_e = new int[dims_e[0]*dims_e[1]];
-  dataset_e.read(data_e, H5::PredType::STD_I32LE);
-
-  size_t nx_e = flag;
-  size_t ny_e = dims_e[0] * dims_e[1] / nx_e;
-  
-  element_node_indices.resize(ny_e);
-
-  for (size_t i = 0; i < ny_e; i ++) {
-    element_node_indices[i].resize(nx_e);
-    for (size_t j = 0; j < nx_e; j ++)
-      element_node_indices[i][j] = data_e[i * nx_e + j];
-  }
-  
-  if (flag == SUBSURFACE_NCOL) {
-    n_cells_subsurface = ny_e;
-    subsurface_cells.resize(ny_e);
-  } else if (flag == SURFACE_NCOL) {
-    n_cells_surface = ny_e;
-    surface_cells.resize(ny_e);
-  }
-  
-}
-
-void Mesh::read_velocity(H5::H5File fp,
-			 std::string groupname,
-			 std::string datasetname,
-			 std::string datasetname_next,
-			 int flag,
-			 int index)
-{
-
-  // read in the current time step data
-  
-  H5::Group group		= fp.openGroup(groupname.c_str());
-  H5::DataSet dataset		= group.openDataSet(datasetname.c_str());
-  H5::DataType datatype		= dataset.getDataType();
-  H5::DataSpace dataspace	= dataset.getSpace();
-  
-  int rank = dataspace.getSimpleExtentNdims();
-
-  hsize_t *dims = new hsize_t[rank];
-  dataspace.getSimpleExtentDims(dims);
-  
-  double *data = new double[dims[0]*dims[1]];
-  dataset.read(data, H5::PredType::IEEE_F64LE);
-
-  if (flag == SUBSURFACE_FLAG) {
-    for (int i = 0; i < n_cells_subsurface; i ++)
-      subsurface_cells[i].set_velocity(index, data[i]);
-  } else if (flag == SURFACE_FLAG) {
-    for (int i = 0; i < n_cells_surface; i ++)
-      surface_cells[i].set_velocity(index, data[i]);
-  }
-
-  // read in the next time step data for linear interpolation
-  
-  H5::DataSet dataset_next	= group.openDataSet(datasetname_next.c_str());
-  H5::DataType datatype_next	= dataset_next.getDataType();
-  H5::DataSpace dataspace_next	= dataset_next.getSpace();
-
-  rank = dataspace_next.getSimpleExtentNdims();
-  
-  dims = new hsize_t[rank];
-  dataspace_next.getSimpleExtentDims(dims);
-
-  data = new double[dims[0]*dims[1]];
-  dataset_next.read(data, H5::PredType::IEEE_F64LE);
-
-  if (flag == SUBSURFACE_FLAG) {
-    for (int i = 0; i < n_cells_subsurface; i ++)
-      subsurface_cells[i].set_velocity_next(index, data[i]);
-  } else if (flag == SURFACE_FLAG) {
-    for (int i = 0; i < n_cells_surface; i ++)
-      surface_cells[i].set_velocity_next(index, data[i]);
-  }
-
-}
-
-void Mesh::read_water_depth(H5::H5File fp,
-			    std::string datasetname,
-			    std::string datasetname_next)
-{
-
-  // read in current time step
-  
-  H5::Group group		= fp.openGroup("surface-ponded_depth.cell.0");
-  H5::DataSet dataset		= group.openDataSet(datasetname.c_str());
-  H5::DataType datatype		= dataset.getDataType();
-  H5::DataSpace dataspace	= dataset.getSpace();
-  
-  int rank = dataspace.getSimpleExtentNdims();
-
-  hsize_t *dims = new hsize_t[rank];
-  dataspace.getSimpleExtentDims(dims);
-  
-  double *data = new double[dims[0]*dims[1]];
-  dataset.read(data, H5::PredType::IEEE_F64LE);
-
-  for (int i = 0; i < n_cells_surface; i ++) {
-    surface_cells[i].set_water_depth(data[i]);
-  }
-
-  // read in next time step
-
-  H5::DataSet dataset_next	= group.openDataSet(datasetname_next.c_str());
-  H5::DataType datatype_next	= dataset_next.getDataType();
-  H5::DataSpace dataspace_next	= dataset_next.getSpace();
-  
-  rank = dataspace_next.getSimpleExtentNdims();
-
-  dims = new hsize_t[rank];
-  dataspace_next.getSimpleExtentDims(dims);
-  
-  data = new double[dims[0]*dims[1]];
-  dataset_next.read(data, H5::PredType::IEEE_F64LE);
-
-  for (int i = 0; i < n_cells_surface; i ++) {
-    surface_cells[i].set_water_depth_next(data[i]);
-  }
-  
-}
-
-void Mesh::read_elevation(H5::H5File fp, std::string datasetname)
-{
-
-  H5::Group group		= fp.openGroup("surface-elevation.cell.0");
-  H5::DataSet dataset		= group.openDataSet(datasetname.c_str());
-  H5::DataType datatype		= dataset.getDataType();
-  H5::DataSpace dataspace	= dataset.getSpace();
-  
-  int rank = dataspace.getSimpleExtentNdims();
-
-  hsize_t *dims = new hsize_t[rank];
-  dataspace.getSimpleExtentDims(dims);
-  
-  double *data = new double[dims[0]*dims[1]];
-  dataset.read(data, H5::PredType::IEEE_F64LE);
-  
-  for (int i = 0; i < n_cells_surface; i ++) {
-    surface_cells[i].set_coordinate(2, data[i]);
-  }
-  
-}
diff --git a/mesh.hpp b/mesh.hpp
index c22cb31..358b2b4 100644
--- a/mesh.hpp
+++ b/mesh.hpp
@@ -4,67 +4,669 @@
 #define MESH_HPP
 
 #include "header.hpp"
-#include "cell.hpp"
+#include "parser.hpp"
+#include "subgrid.hpp"
+#include "tooling.hpp"
+#include "column.hpp"
 
+// Content
+// =======
+//
+// 1. Mesh
+// ---
+// 1.1 compile
+// 1.2 dump
+// 1.3 initialise_subgrid
+// 1.4 map_surface
+// 1.5 map_subsurface
+// ---
+// 1.6 get_bounding_box
+// 1.7 read_node_indices
+// 1.8 read_nodes
+// 1.9 initialise_surface
+// 1.10 initialise_subsurface
+//
+
+// ------------------------------------------------------------------------
+// Mesh
+// ------------------------------------------------------------------------
+// The mesh is a collection of columns. The columns needn't to know
+// their neighbours.
+// ------------------------------------------------------------------------
 class Mesh
 {
 
 public:
 
-  Mesh();
+    std::vector<Column> columns;
+
+    Mesh(Parser &value)
+    {
+        parser = value;
+    }
+
+    // ---
+    // re-build the mesh and the subgrid from what little
+    // information is contained in the ats output files
+    // ---
+    void compile()
+    {
+
+        std::cout << "[msh] compiling the mesh" << std::endl;
+
+        H5::H5File fp_surface    = H5::H5File(parser.surface_mesh_file.c_str(),
+                                              H5F_ACC_RDONLY);
+        H5::H5File fp_subsurface = H5::H5File(parser.subsurface_mesh_file.c_str(),
+                                              H5F_ACC_RDONLY);
+
+
+        //
+
+        std::cout << "[msh] surface mesh indexing" << std::endl;
+
+        std::vector<std::vector<int>> node_indices;
+        read_node_indices(fp_surface, node_indices);
+        n_vertices = node_indices[0].size();
+
+        std::vector<point> nodes;
+        read_nodes(fp_surface, nodes);
+
+        initialise_surface(node_indices, nodes);
+
+        //
+
+
+        //
+
+        std::cout << "[msh] subsurface mesh indexing" << std::endl;
+
+        for (size_t i = 0; i < node_indices.size(); i ++)
+            node_indices[i].clear();
+        read_node_indices(fp_subsurface, node_indices);
+
+        nodes.clear();
+        read_nodes(fp_subsurface, nodes);
+
+        //
+
+
+        std::cout << "[msh] bounding box: {";
+        std::cout << low_coordinates[0] << ", " << low_coordinates[1];
+        std::cout << "} {";
+        std::cout << high_coordinates[0] << ", " << high_coordinates[1];
+        std::cout << "}" << std::endl;
+
+
+        //
+
+        int nz = probe_column_depth(node_indices, nodes);
+        initialise_subgrid(nz);
+
+        map_surface();
+
+        initialise_subsurface(node_indices, nodes);
+
+        map_subsurface();
+
+        std::cout << std::endl;
+
+    }
+
+    // ---
+    // dump the subgrid with all topology and geometry information
+    // ---
+    void dump()
+    {
+
+        std::cout << "[msh] dumping the subgrid" << std::endl;
+
+        int nx = sgrid.indices.size();
+        int ny = sgrid.indices[0].size();
+        int nz = sgrid.indices[0][0].size();
+
+        int *subgrid_indices_array = new int[nx * ny * nz];
+        double *subgrid_elevations_array = new double[nx * ny * nz];
+        double *depths_array = new double[nz-1];
+
+        for (int i = 0; i < nx; ++i) {
+            for (int j = 0; j < ny; ++j) {
+                for (int k = 0; k < nz; ++k) {
+                    subgrid_indices_array[i + nx * (j + ny * k)] = sgrid.indices[i][j][k];
+                    subgrid_elevations_array[i + nx * (j + ny * k)] = sgrid.elevations[i][j][k];
+                }
+            }
+        }
+
+        for (int i = 0; i < nz-1; ++i) {
+            depths_array[i] = sgrid.dzs[i];
+        }
+
+        // open a hdf5 file and write stuff to file
+        H5::H5File file("grid.h5", H5F_ACC_TRUNC);
+
+        hsize_t rank = 1;
+        hsize_t dims[1];
+
+        // ---
+        // write cellsize
+        // ---
+        double cellsize[1] = {parser.subgrid_dx};
+        dims[0] = 1;
+
+        H5::DataSpace dspace_cellsize(rank, dims);
+        H5::FloatType dtype_cellsize(H5::PredType::NATIVE_DOUBLE);
+        dtype_cellsize.setOrder(H5T_ORDER_LE);
+
+        H5::DataSet dset_cellsize = file.createDataSet("cellsize", dtype_cellsize, dspace_cellsize);
+        dset_cellsize.write(cellsize, H5::PredType::NATIVE_DOUBLE);
+
+        // ---
+        // write the lower left corner
+        // ---
+        double upperleft[2] = {low_coordinates[0], high_coordinates[1]};
+        dims[0] = 2;
+
+        H5::DataSpace dspace_upperleft(rank, dims);
+        H5::FloatType dtype_upperleft(H5::PredType::NATIVE_DOUBLE);
+        dtype_upperleft.setOrder(H5T_ORDER_LE);
+
+        H5::DataSet dset_upperleft = file.createDataSet("upperleft", dtype_upperleft, dspace_upperleft);
+        dset_upperleft.write(upperleft, H5::PredType::NATIVE_DOUBLE);
+
+        // ---
+        // write depths
+        // ---
+        dims[0] = nz-1;
+
+        H5::DataSpace dspace_depths(rank, dims);
+        H5::FloatType dtype_depths(H5::PredType::NATIVE_DOUBLE);
+        dtype_depths.setOrder(H5T_ORDER_LE);
+
+        H5::DataSet dset_depths = file.createDataSet("depths", dtype_depths, dspace_depths);
+        dset_depths.write(depths_array, H5::PredType::NATIVE_DOUBLE);
+
+        // ---
+        // write dimensions
+        // ---
+        int dimensions[3] = {nx, ny, nz};
+        dims[0] = 3;
+
+        H5::DataSpace dspace_dims(rank, dims);
+        H5::IntType dtype_dims(H5::PredType::NATIVE_INT);
+        dtype_dims.setOrder(H5T_ORDER_LE);
+
+        H5::DataSet dset_dims = file.createDataSet("dimensions", dtype_dims, dspace_dims);
+        dset_dims.write(dimensions, H5::PredType::NATIVE_INT);
+
+        dims[0] = nx * ny * nz;
+
+        // ---
+        // write indices
+        // ---
+        H5::DataSpace dspace_indices(rank, dims);
+        H5::IntType dtype_indices(H5::PredType::NATIVE_INT);
+        dtype_indices.setOrder(H5T_ORDER_LE);
+
+        H5::DataSet dset_indices = file.createDataSet("indices", dtype_indices, dspace_indices);
+        dset_indices.write(subgrid_indices_array, H5::PredType::NATIVE_INT);
+
+        // ---
+        // write elevations
+        // ---
+        H5::DataSpace dspace_elevations(rank, dims);
+        H5::FloatType dtype_elevations(H5::PredType::NATIVE_DOUBLE);
+        dtype_elevations.setOrder(H5T_ORDER_LE);
+
+        H5::DataSet dset_elevations = file.createDataSet("elevations", dtype_elevations, dspace_elevations);
+        dset_elevations.write(subgrid_elevations_array, H5::PredType::NATIVE_DOUBLE);
+
+        file.close();
+
+        // free memory
+
+        delete []subgrid_indices_array;
+        delete []subgrid_elevations_array;
+        delete []depths_array;
+
+        std::cout << "[msh] subgrid dumped to grid.h5" << std::endl;
+
+    }
+
+    // ---
+    // initialise the subgrid, compute its extent and how many grid
+    // points are needed
+    // ---
+    void initialise_subgrid(int nz)
+    {
+
+        std::vector<double> depths(nz);
+        std::cout << "[msh] initialising the subgrid" << std::endl;
+
+        for (int i = 0; i < nz; i ++)
+            depths[i] = columns[0].stack_upper_horizon[i] - columns[0].stack_lower_horizon[i];
+
+        sgrid = Subgrid(parser, low_coordinates, high_coordinates, depths);
+
+    }
+
+    void get_bounding_box(std::vector<point> &vertices, double &xmin, double &xmax, double &ymin, double &ymax)
+    {
+
+        xmin = DBL_MAX;
+        ymin = DBL_MAX;
+        xmax = DBL_MIN;
+        ymax = DBL_MIN;
+
+        for (size_t k = 0; k < vertices.size(); k ++) {
+
+            xmin = std::min(vertices[k][0], xmin);
+            xmax = std::max(vertices[k][0], xmax);
+
+            ymin = std::min(vertices[k][1], ymin);
+            ymax = std::max(vertices[k][1], ymax);
+
+        }
+
+    }
+
+    // ---
+    // assign each subgrid subsurface point the index of the
+    // subsurface cell containing it
+    // ---
+    void map_subsurface()
+    {
+
+        size_t i;
+        std::cout << "[msh] mapping subsurface cells to subgrid...";
+
+        double offset_x = low_coordinates[0];
+        double offset_y = high_coordinates[1];
+        double gridsize = parser.subgrid_dx;
+
+        for (i = 0; i < columns.size(); ++ i) {
+
+            // get bounding box
+            double xmin = DBL_MAX;
+            double ymin = DBL_MAX;
+            double xmax = DBL_MIN;
+            double ymax = DBL_MIN;
+
+            get_bounding_box(columns[i].vertices, xmin, xmax, ymin, ymax);
+
+            int colmin = (int) floor((xmin - offset_x) / gridsize);
+            int rowmin = (int) floor((offset_y - ymax) / gridsize);
+            int colmax = (int) floor((xmax - offset_x) / gridsize);
+            int rowmax = (int) floor((offset_y - ymin) / gridsize);
+
+            // check limits
+            colmax = std::min(colmax, sgrid.nx - 1);
+            rowmax = std::min(rowmax, sgrid.ny - 1);
+
+            for (int ix = colmin; ix <= colmax; ++ ix) {
+                for (int iy = rowmin; iy <= rowmax; ++ iy) {
+
+                    point p1 { offset_x + 0.5 * gridsize + ix * gridsize,
+                               offset_y - 0.5 * gridsize - iy * gridsize, 0.0 };
+
+                    if (is_inside2d(p1, columns[i].centroid, columns[i].vertices)) {
+
+                        for (size_t iz = 1; iz < columns[i].stack_index.size(); ++ iz) {
+
+                            // map subsurface
+                            sgrid.indices[ix][iy][iz] = columns[i].stack_index[iz];
+                            sgrid.elevations[ix][iy][iz] = 0.5 * (columns[i].stack_lower_horizon[iz] + columns[i].stack_upper_horizon[iz]);
+
+                        }
+
+                    }
+
+                }
+            }
+
+        }
+
+        std::cout << " done" << std::endl;
+
+    }
 
-  void build_surface_mesh(H5::H5File fp, H5::H5File ep, std::string datasetname);
-  void build_subsurface_mesh(H5::H5File fp);
+    // ---
+    // assign each subgrid surface point the index of the surface cell
+    // containing it
+    // ---
+    void map_surface()
+    {
 
-  void build_environment(H5::H5File ep,
-			 std::string datasetname,
-			 std::string datasetname_next,
-			 int flag);
-  
-  //  void build_subsurface_environment(H5::H5File ep);
+        // the strategy is to assign the top layer of the subgrid to the
+        // corresponding cell in a 2d manner.  the underlying cells can be
+        // mapped using indices and offset computation.
 
-  Cell* get_cell(int index, int flag);
+        size_t i;
+        std::cout << "[msh] mapping surface cells to subgrid...";
 
-  void set_high_coordinate(int index, double value);
-  double get_high_coordinate(int index);
-  void set_low_coordinate(int index, double value);
-  double get_low_coordinate(int index);
-  void set_fuzz(double value);
+        // index [0,0] corresponds to upper left corner of the domain.
+        // we compute the offset in x and y direction accordingly.
+        double offset_x = low_coordinates[0];
+        double offset_y = high_coordinates[1];
+        double gridsize = parser.subgrid_dx;
 
-  int get_n_cells(int flag);
+        #pragma omp parallel
+        {
+            #pragma omp for private(i)
+            for (i = 0; i < columns.size(); ++ i) {
+
+                // get bounding box
+                double xmin = DBL_MAX;
+                double ymin = DBL_MAX;
+                double xmax = DBL_MIN;
+                double ymax = DBL_MIN;
+
+                get_bounding_box(columns[i].vertices, xmin, xmax, ymin, ymax);
+
+                int colmin = (int) floor((xmin - offset_x) / gridsize);
+                int rowmin = (int) floor((offset_y - ymax) / gridsize);
+                int colmax = (int) floor((xmax - offset_x) / gridsize);
+                int rowmax = (int) floor((offset_y - ymin) / gridsize);
+
+                // check limits
+                colmax = std::min(colmax, sgrid.nx - 1);
+                rowmax = std::min(rowmax, sgrid.ny - 1);
+
+                for (int ix = colmin; ix <= colmax; ++ ix) {
+                    for (int iy = rowmin; iy <= rowmax; ++ iy) {
+
+                        point p1 { offset_x + 0.5 * gridsize + ix * gridsize, offset_y - 0.5 * gridsize - iy * gridsize, 0.0 };
+
+                        if (is_inside2d(p1, columns[i].centroid, columns[i].vertices)) {
+
+                            // map surface
+                            sgrid.indices[ix][iy][0] = i;
+                            sgrid.elevations[ix][iy][0] = columns[i].centroid[2];
+
+                        }
+
+                    }
+                }
+
+            }
+
+        } // end of parallel region
+
+        std::cout << " done" << std::endl;
+
+    }
 
 private:
 
-  int n_cells_subsurface;
-  int n_cells_surface;
-
-  double fuzz;
-  
-  std::vector<Cell> subsurface_cells;
-  std::vector<Cell> surface_cells;
-
-  std::vector<double> high_coordinates;
-  std::vector<double> low_coordinates;
-
-  std::vector<std::vector<double>> nodes;
-  std::vector<std::vector<int>> element_node_indices;
-  
-  void handle_subsurface_node_data(H5::H5File fp);
-  void build_subsurface_cells();
-
-  void handle_surface_node_data(H5::H5File fp);
-  void build_surface_cells(H5::H5File fp, std::string datasetname);
-
-  void get_mesh_metadata(H5::H5File fp, int ncol);
-  void read_velocity(H5::H5File fp,
-		     std::string groupname,
-		     std::string datasetname,
-		     std::string datasetname_next,
-		     int flag,
-		     int index);
-  void read_water_depth(H5::H5File fp, std::string datasetname, std::string datasetname_next);
-  void read_elevation(H5::H5File fp, std::string datasetname);
-  
+    Parser parser;
+    Subgrid sgrid;
+    point high_coordinates;
+    point low_coordinates;
+
+    int n_columns;
+    int n_vertices;
+
+    // ---
+    // read the node indices from hdf5 file
+    // ---
+    void read_node_indices(H5::H5File fp, std::vector<std::vector<int>> &node_indices)
+    {
+
+        H5::DataSet dset = fp.openDataSet(MESH_MIXED_ELEMENTS);
+        H5::DataType dtype = dset.getDataType();
+        H5::DataSpace dspace = dset.getSpace();
+
+        int rank = dspace.getSimpleExtentNdims();
+        hsize_t* dims = new hsize_t[rank];
+        dspace.getSimpleExtentDims(dims);
+
+        int n_dims = 1;
+        for (int i = 0; i < rank; i++)
+            n_dims *= dims[i];
+
+        int* data = new int[n_dims];
+        dset.read(data, H5::PredType::STD_I32LE);
+
+        // columns and vertices change from surface/subsurface
+        size_t ncols = data[0];
+        size_t nvers = 0;
+        if (ncols == 4) {
+            nvers = 3;
+        } else if (ncols == 8) {
+            ncols = 7;
+            nvers = 6;
+        }
+        size_t nrows = n_dims / ncols;
+        //
+
+        std::cout << "[msh] number of vertices per cell: " << nvers << std::endl;
+
+        node_indices.resize(nrows);
+
+        for (size_t i = 0; i < nrows; i ++) {
+            node_indices[i].resize(nvers);
+            for (size_t j = 0; j < nvers; j ++)
+                node_indices[i][j] = data[i * ncols + j + 1];
+        }
+
+        n_columns = nrows;
+
+    }
+
+    // ---
+    // read in the coordinates from the hdf5 file
+    // ---
+    void read_nodes(H5::H5File fp, std::vector<point> &nodes)
+    {
+
+        H5::DataSet dset = fp.openDataSet(MESH_NODES);
+        H5::DataType dtype = dset.getDataType();
+        H5::DataSpace dspace = dset.getSpace();
+
+        int rank = dspace.getSimpleExtentNdims();
+
+        hsize_t *dims = new hsize_t[rank];
+        dspace.getSimpleExtentDims(dims);
+
+        int n_dims = 1;
+        for (int i = 0; i < rank; i++)
+            n_dims *= dims[i];
+
+        double *data = new double[n_dims];
+        dset.read(data, H5::PredType::IEEE_F64LE);
+
+        nodes.resize(dims[0]);
+
+        for (size_t i = 0; i < dims[0]; i ++) {
+            for (size_t j = 0; j < dims[1]; j ++)
+                nodes[i][j] = data[i * dims[1] + j];
+        }
+
+    }
+
+    // ---
+    // initialise the top cell (vertices and centroids)
+    // ---
+    void initialise_surface(std::vector<std::vector<int>> &node_indices, std::vector<point> &nodes)
+    {
+
+        columns.resize(node_indices.size());
+
+        for (int i = 0; i < DIM; i ++) {
+            high_coordinates[i] = DBL_MIN;
+            low_coordinates[i] = DBL_MAX;
+        }
+
+        for (size_t i = 0; i < node_indices.size(); i ++) {
+
+            columns[i].vertices.resize(n_vertices);
+
+            for (int j = 0; j < n_vertices; j ++) {
+
+                int index = node_indices[i][j];
+
+                for (size_t k = 0; k < DIM; k ++) {
+
+                    double val = nodes[index][k];
+
+                    // bounding box for the entire domain
+                    high_coordinates[k] = std::max(val, high_coordinates[k]);
+                    low_coordinates[k]  = std::min(val, low_coordinates[k]);
+                    //
+
+                    columns[i].vertices[j][k] = val;
+                    columns[i].centroid[k] += val / (double) n_vertices;
+
+                }
+            }
+        }
+
+    }
+
+    int probe_column_depth(std::vector<std::vector<int>> &node_indices, std::vector<point> &nodes)
+    {
+
+        int maxstack = 0;
+
+        for (size_t i = 0; i < node_indices.size(); ++ i) {
+
+            point centroid;
+            for (size_t j = 0; j < DIM; ++ j) {
+                centroid[j] = 0.0;
+            }
+
+            size_t index_count = node_indices[i].size();
+
+            for (size_t j = 0; j < index_count; ++ j) {
+                int index = node_indices[i][j];
+                for (size_t k = 0; k < DIM; ++ k) {
+                    double val = nodes[index][k];
+                    centroid[k] += val / (double) index_count;
+                }
+            }
+
+            double distance = get_distance2d(columns[0].centroid, centroid);
+            if (distance <= parser.fuzz) {
+
+                columns[0].stack_index.push_back(i);
+
+                // determine horizons
+                point p = nodes[node_indices[i][0]];
+                double zmax = p[2];
+                double zmin;
+
+                for (size_t k = 1; k < node_indices[i].size(); ++k) {
+                    point q = nodes[node_indices[i][k]];
+                    if (get_distance2d(p, q) < parser.fuzz) {
+                        if (zmax < q[2]) {
+                            zmin = zmax;
+                            zmax = q[2];
+                        } else {
+                            zmin = q[2];
+                        }
+                        break;
+                    }
+                }
+
+                columns[0].stack_upper_horizon.push_back(zmax);
+                columns[0].stack_lower_horizon.push_back(zmin);
+
+                maxstack += 1;
+            }
+
+        }
+
+        columns[0].sort();
+
+        return maxstack;
+
+    }
+
+    // ---
+    // initialise the stack, calculating the elevations
+    // and node indices in the stack
+    // ---
+    void initialise_subsurface(std::vector<std::vector<int>> &node_indices, std::vector<point> &nodes)
+    {
+
+        size_t i;
+
+        double offset_x = low_coordinates[0];
+        double offset_y = high_coordinates[1];
+        double gridsize = parser.subgrid_dx;
+
+        #pragma omp parallel
+        {
+            #pragma omp for private(i)
+            for (i = 0; i < node_indices.size(); ++ i) {
+
+                // compute the centroid
+                point centroid = {0.0, 0.0, 0.0};
+                size_t index_count = node_indices[i].size();
+
+                // vertical bounds
+                double zmax = DBL_MIN;
+                double zmin = DBL_MAX;
+
+                for (size_t j = 0; j < index_count; ++ j) {
+
+                    int index = node_indices[i][j];
+
+                    double ztmp = nodes[index][2];
+                    zmax = std::max(zmax, ztmp);
+                    zmin = std::min(zmin, ztmp);
+
+                    for (size_t k = 0; k < DIM; ++ k) {
+                        double val = nodes[index][k];
+                        centroid[k] += val / (double) index_count;
+                    }
+
+                }
+
+                //
+
+                // convert into indices and add to stack of column `index`
+                int ix = (int) floor((centroid[0] - offset_x) / gridsize);
+                int iy = (int) floor((offset_y - centroid[1]) / gridsize);
+                int index = sgrid.indices[ix][iy][0];
+
+                if (index < 0) {
+
+                    printf("[msh] process %d: cannot find a column for cell %ld (%d, %d), centroid {%.12f %.12f}\n",
+                           omp_get_thread_num(),
+                           i, ix, iy, centroid[0], centroid[1]);
+                    exit(EXIT_FAILURE);
+
+                }
+
+                if (index > 0) {
+                    columns[index].stack_index.push_back(i);
+                    columns[index].stack_upper_horizon.push_back(zmax);
+                    columns[index].stack_lower_horizon.push_back(zmin);
+                }
+                //
+
+            }
+
+            // sanity check and sorting
+            #pragma omp for private(i)
+            for (i = 1; i < columns.size(); ++ i) {
+
+                if ((int) columns[i].stack_index.size() != sgrid.nz) {
+                    printf("[msh] process %d: incomplete column %ld with %ld/%d subsurface cells\n",
+                           omp_get_thread_num(),
+                           i,
+                           columns[i].stack_index.size(), sgrid.nz);
+                    exit(EXIT_FAILURE);
+                }
+
+                columns[i].sort();
+            }
+
+
+        }
+
+    }
+
 };
 
 #endif
diff --git a/parser.cpp b/parser.cpp
deleted file mode 100644
index 39527ed..0000000
--- a/parser.cpp
+++ /dev/null
@@ -1,114 +0,0 @@
-#include "header.hpp"
-#include "parser.hpp"
-
-void ParserLine::parse()
-{
-  
-  key.clear();
-  value.clear();
-  
-  if (!line.empty() && line.find("//", 0) != 0) {
-
-    uint splitloc = line.find(":", 0); // find the colon
-    key = line.substr(0, splitloc);    // store the key and value strings
-    key.erase(std::remove(key.begin(), key.end(), ' '), key.end()); // remove spaces and tabs
-    key.erase(std::remove(key.begin(), key.end(), '\t'), key.end());
-    std::string val = line.substr(splitloc+1, line.length()-splitloc);
-
-    size_t splitter = val.find("//", 0);
-    std::string strloc;
-    if (splitter != std::string::npos) {
-      strloc = val.substr(0, splitter);
-    } else {
-      strloc = val;
-    }
-
-    value.clear();
-    value.str(strloc);
-    
-  }
-  
-}
-
-
-void Parser::read_options(std::string fname)
-{
-
-  std::ifstream fInStream(fname);
-  std::string line;
-  ParserLine pline;
-  
-  // initialize all read-in values to -999
-  fuzz    = -999;
-  n_agent = -999;
-  n_iteration = -999;
-
-  if (fInStream.is_open()) {
-    
-    while (std::getline(fInStream, line)) {
-
-      pline.line = line;
-      pline.parse();
-
-      if (!pline.key.empty()) {
-	
-	// Match the key and store the value
-	if (!strcmp("surfacemesh", pline.key.c_str())) {
-	  pline.value >> surface_mesh_file;
-	}
-	else if (!strcmp("surfacedata", pline.key.c_str())) {
-	  pline.value >> surface_data_file;
-	}
-	else if (!strcmp("surfacetime", pline.key.c_str())) {
-	  double time;
-	  while (pline.value >> time) {
-	    surface_time.push_back(time);
-	  }
-	}
-	else if (!strcmp("surfacegroup", pline.key.c_str())) {
-	  std::string group;
-	  while (pline.value >> group) {
-	    surface_data_group.push_back(group);
-	  }
-	}
-	else if (!strcmp("subsurfacemesh", pline.key.c_str())) {
-	  pline.value >> subsurface_mesh_file;
-	}
-	else if (!strcmp("subsurfacedata", pline.key.c_str())) {
-	  pline.value >> subsurface_data_file;
-	}
-	else if (!strcmp("subsurfacetime", pline.key.c_str())) {
-	  double time;
-	  while (pline.value >> time) {
-	    subsurface_time.push_back(time);
-	  }
-	}
-	else if (!strcmp("subsurfacegroup", pline.key.c_str())) {
-	  std::string group;
-	  while (pline.value >> group) {
-	    subsurface_data_group.push_back(group);
-	  }
-	}
-	else if (!strcmp("n_agent", pline.key.c_str())) {
-	  pline.value >> n_agent;
-	}
-	else if (!strcmp("fuzz", pline.key.c_str())) {
-	  pline.value >> fuzz;
-	}
-	else if (!strcmp("seed", pline.key.c_str())) {
-	  pline.value >> seed;
-	}
-	else if (!strcmp("n_iteration", pline.key.c_str())) {
-	  pline.value >> n_iteration;
-	}
-	else {
-	  std::cerr << "[!!] Unrecognised option key" << std::endl;
-	}
-	
-      }
-      
-    }
-    
-  }
-  
-}
diff --git a/parser.hpp b/parser.hpp
index 7ffda0d..fa4e3ab 100644
--- a/parser.hpp
+++ b/parser.hpp
@@ -10,17 +10,57 @@
 
 #include "header.hpp"
 
+// Content
+// =======
+//
+// 1. ParserLine
+// ---
+// 1.1 parse
+//
+// 2. Parser
+// ---
+// 2.1 print_options
+// 2.2 read_options
+//
+
 class ParserLine
 {
 
 public:
 
-  std::string line;
-  std::string key;
-  std::stringstream value;
+    std::string line;
+    std::string key;
+    std::stringstream value;
+
+    void parse()
+    {
+
+        key.clear();
+        value.clear();
+
+        if (!line.empty() && line.find("//", 0) != 0) {
+
+            uint splitloc = line.find(":", 0); // find the colon
+            key = line.substr(0, splitloc);    // store the key and value strings
+            key.erase(std::remove(key.begin(), key.end(), ' '), key.end()); // remove spaces and tabs
+            key.erase(std::remove(key.begin(), key.end(), '\t'), key.end());
+            std::string val = line.substr(splitloc+1, line.length()-splitloc);
+
+            size_t splitter = val.find("//", 0);
+            std::string strloc;
+            if (splitter != std::string::npos) {
+                strloc = val.substr(0, splitter);
+            } else {
+                strloc = val;
+            }
+
+            value.clear();
+            value.str(strloc);
+
+        }
+
+    }
 
-  void parse();
-  
 };
 
 class Parser
@@ -28,28 +68,189 @@ class Parser
 
 public:
 
-  std::string surface_mesh_file;
-  std::string surface_data_file;
-  std::string subsurface_mesh_file;
-  std::string subsurface_data_file;
-
-  std::vector<std::string> surface_data_group;
-  std::vector<double> surface_time;
-  
-  std::vector<std::string> subsurface_data_group;
-  std::vector<double> subsurface_time;
-
-  int n_agent;
-  int n_iteration;
-  int seed;
-  
-  int surface_np;
-  int subsurface_np;
-  
-  double fuzz;
-
-  void read_options(std::string fname);
-  
+    std::string surface_mesh_file;
+    std::string surface_data_file;
+    std::string subsurface_mesh_file;
+    std::string subsurface_data_file;
+
+    std::vector<std::string> surface_data_group;
+    std::vector<double> surface_time;
+
+    std::vector<std::string> subsurface_data_group;
+    std::vector<double> subsurface_time;
+
+    int n_agent;
+    int n_iteration;
+    int seed;
+
+    double start_time = 0.0;
+    double end_time   = 1.0;
+
+    double fuzz;
+    double subgrid_dx;
+
+    void print_options()
+    {
+
+        std::cout << "[prs] number of agents: " << n_agent << std::endl;
+        std::cout << "[prs] number of iterations: " << n_iteration << std::endl;
+        std::cout << "[prs] random seed: " << seed << std::endl;
+        std::cout << "[prs] geometric fuzz: " << fuzz << std::endl;
+        std::cout << "[prs] subgrid spacing: " << subgrid_dx << std::endl;
+
+        std::cout << "[prs] surface mesh file: " << surface_mesh_file << std::endl;
+        std::cout << "[prs] surface data file: " << surface_data_file << std::endl;
+        std::cout << "[prs] surface time: ";
+        for (uint i = 0; i < surface_time.size(); i ++)
+            std::cout << surface_time[i] << " ";
+        std::cout <<std::endl;
+        std::cout << "[prs] surface group: ";
+        for (uint i = 0; i < surface_data_group.size(); i ++)
+            std::cout << surface_data_group[i] << " ";
+        std::cout <<std::endl;
+
+        std::cout << "[prs] subsurface mesh file: " << subsurface_mesh_file << std::endl;
+        std::cout << "[prs] subsurface data file: " << subsurface_data_file << std::endl;
+        std::cout << "[prs] subsurface time: ";
+        for (uint i = 0; i < subsurface_time.size(); i ++)
+            std::cout << subsurface_time[i] << " ";
+        std::cout <<std::endl;
+        std::cout << "[prs] subsurface group: ";
+        for (uint i = 0; i < subsurface_data_group.size(); i ++)
+            std::cout << subsurface_data_group[i] << " ";
+        std::cout << std::endl;
+        std::cout << std::endl;
+
+    }
+
+    void read_options(std::string fname)
+    {
+
+        std::ifstream fInStream(fname);
+        std::string line;
+        ParserLine pline;
+
+        // initialize all read-in values to -999
+        fuzz    = -999;
+        n_agent = -999;
+        n_iteration = -999;
+
+        if (fInStream.is_open()) {
+
+            while (std::getline(fInStream, line)) {
+
+                pline.line = line;
+                pline.parse();
+
+                if (!pline.key.empty()) {
+
+                    // Match the key and store the value
+                    if (!strcmp("surfacemesh", pline.key.c_str())) {
+                        pline.value >> surface_mesh_file;
+                    } else if (!strcmp("surfacedata", pline.key.c_str())) {
+                        pline.value >> surface_data_file;
+                    } else if (!strcmp("surfacetime", pline.key.c_str())) {
+
+		      std::string timefile;
+		      pline.value >> timefile;
+
+		      std::ifstream in(timefile);
+		      double time;
+
+		      if (in.is_open()) {
+			
+			while (!in.eof()) {
+			  in >> time;
+			  surface_time.push_back(time);
+			}
+
+		      } else {
+			std::cerr << "[err] file " << timefile << " can't be opened" << std::endl;
+		      }
+		      
+                    } else if (!strcmp("surfacegroup", pline.key.c_str())) {
+
+		      std::string groupfile;
+		      pline.value >> groupfile;
+
+		      std::ifstream in(groupfile);
+		      std::string group;
+
+		      if (in.is_open()) {
+
+			while (!in.eof()) {
+			  in >> group;
+			  surface_data_group.push_back(group);
+			}
+
+		      } else {
+			std::cerr << "[err] file " << groupfile << " can't be opened" << std::endl;
+		      }
+
+                    } else if (!strcmp("subsurfacemesh", pline.key.c_str())) {
+                        pline.value >> subsurface_mesh_file;
+                    } else if (!strcmp("subsurfacedata", pline.key.c_str())) {
+                        pline.value >> subsurface_data_file;
+                    } else if (!strcmp("subsurfacetime", pline.key.c_str())) {
+
+		      std::string timefile;
+		      pline.value >> timefile;
+
+		      std::ifstream in(timefile);
+		      double time;
+
+		      if (in.is_open()) {
+		      
+			while (!in.eof()) {
+			  in >> time;
+			  subsurface_time.push_back(time);
+			}
+
+		      } else {
+			std::cerr << "[err] file " << timefile << " can't be opened" << std::endl;
+		      }
+		      
+                    } else if (!strcmp("subsurfacegroup", pline.key.c_str())) {
+
+		      std::string groupfile;
+		      pline.value >> groupfile;
+
+		      std::ifstream in(groupfile);
+		      std::string group;
+
+		      if (in.is_open()) {
+
+			while (!in.eof()) {
+			  in >> group;
+			  subsurface_data_group.push_back(group);
+			}
+
+		      } else {
+			std::cerr << "[err] file " << groupfile << " can't be opened" << std::endl;
+		      }
+		      
+                    } else if (!strcmp("n_agent", pline.key.c_str())) {
+                        pline.value >> n_agent;
+                    } else if (!strcmp("fuzz", pline.key.c_str())) {
+                        pline.value >> fuzz;
+                    } else if (!strcmp("subgrid_dx", pline.key.c_str())) {
+                        pline.value >> subgrid_dx;
+                    } else if (!strcmp("seed", pline.key.c_str())) {
+                        pline.value >> seed;
+                    } else if (!strcmp("n_iteration", pline.key.c_str())) {
+                        pline.value >> n_iteration;
+                    } else {
+                        std::cerr << "[err] unrecognised option key" << std::endl;
+                    }
+
+                }
+
+            }
+
+        }
+
+    }
+
 };
 
 #endif
diff --git a/particle.hpp b/particle.hpp
new file mode 100644
index 0000000..0f5026a
--- /dev/null
+++ b/particle.hpp
@@ -0,0 +1,436 @@
+// -*- mode: c++ -*-
+
+#ifndef PARTICLE_HPP
+#define PARTICLE_HPP
+
+#include "header.hpp"
+#include "parser.hpp"
+#include "environment.hpp"
+
+// ------------------------------------------------------------------------
+// Particle
+// ------------------------------------------------------------------------
+// The particle advects passively based on the state of the environment.
+// It can infiltrate, exfiltrate, and evaporate.  As it is advected, it
+// keeps track of its path and its time inside the domain.  At the end
+// of the particle-tracking run, it knows its origin and its faith.
+// ------------------------------------------------------------------------
+class Particle
+{
+
+public:
+
+    int cell_index;
+    double t, x, y, z, dx, dy, dz;
+
+    std::vector<int> domain;
+    std::vector<point> coordinates;
+    std::vector<double> time;
+
+    Particle()
+    {
+        t = x = y = z = 0.0;
+    }
+
+    Particle(double t0, double x0, double y0, double z0)
+    {
+        t = t0;
+        x = x0;
+        y = y0;
+        z = z0;
+    }
+
+    void advect(double vx, double vy, double vz, double dt)
+    {
+
+        vx = check_zero(vx);
+        vy = check_zero(vy);
+        vz = check_zero(vz);
+
+        dx = dy = dz = 0.0;
+
+        dx = vx * dt;
+        dy = vy * dt;
+        dz = vz * dt;
+
+        t += dt;
+        x += dx;
+        y += dy;
+        z += dz;
+
+    }
+
+    void dump(int index)
+    {
+
+        char output_fname[1024];
+        sprintf(output_fname, "output/p-%06d.csv", index);
+
+        FILE *fp;
+
+        fp = fopen(output_fname, "w");
+
+        fprintf(fp, "#time\tdomain\tx\ty\tz\n");
+
+        for (size_t i = 0; i < time.size(); i ++) {
+            fprintf(fp, "%f\t%d\t%f\t%f\t%f\n", time[i],
+                    domain[i],
+                    coordinates[i][0], coordinates[i][1],
+                    coordinates[i][2]);
+        }
+
+        fclose(fp);
+
+    }
+
+};
+
+// ------------------------------------------------------------------------
+// Particle Manager
+// ------------------------------------------------------------------------
+//
+// ------------------------------------------------------------------------
+class ParticleManager
+{
+
+public:
+
+    std::vector<Particle> particles;
+    Environment environment;
+
+    std::default_random_engine rng;
+    std::uniform_real_distribution<double> dist;
+
+    ParticleManager(Parser &value)
+    {
+
+        parser = value;
+        environment = Environment(parser);
+
+        double dx = parser.subgrid_dx;
+
+        std::uniform_int_distribution<int> xdist(0, environment.nx - 1);
+        std::uniform_int_distribution<int> ydist(0, environment.ny - 1);
+        std::uniform_int_distribution<int> zdist(1, environment.nz - 1);
+
+        rng = std::default_random_engine(parser.seed);
+        dist = std::uniform_real_distribution<double>(0.0, 1.0);
+
+        for (int i = 0; i < parser.n_agent; ++ i) {
+
+            Particle p;
+
+            p.cell_index = -9999;
+            int ix, iy, iz;
+
+            while (p.cell_index == -9999) {
+
+                ix = xdist(rng);
+                iy = ydist(rng);
+                iz = zdist(rng) * 0;
+
+                p.cell_index = environment.indices[ix][iy][iz];
+
+            }
+
+            // get the coordinates
+            p.x = environment.xll_corner + 0.5 * dx + ix * dx;
+            p.y = environment.yhl_corner - 0.5 * dx - iy * dx;
+            p.z = environment.elevations[ix][iy][iz];
+            p.t = parser.start_time;
+
+            if (iz > 0) {
+                // if the particle is not in the top layer, it's in the subsurface
+                p.domain.push_back(SUBSURFACE_FLAG);
+            } else {
+                // otherwise, it's at the surface
+                p.domain.push_back(SURFACE_FLAG);
+            }
+
+            p.coordinates.push_back({ p.x, p.y, p.z });
+            p.time.push_back(p.t);
+
+            particles.push_back(p);
+
+        }
+
+    }
+
+    void advect(double timehorizon)
+    {
+
+        std::cout << "[ptl] advecting particles" << std::endl;
+
+        size_t i;
+
+        #pragma omp parallel
+        {
+
+            #pragma omp for private(i)
+            for (i = 0; i < particles.size(); ++i) {
+
+                while (particles[i].time.back() < timehorizon) {
+
+                    int flag = particles[i].domain.back();
+
+                    switch(flag) {
+
+                    case SURFACE_FLAG:
+                        advecth(particles[i], timehorizon);
+                        break;
+
+                    case SUBSURFACE_FLAG:
+                        advectv(particles[i], timehorizon);
+                        break;
+
+                    default:
+                        double dt = std::min(DT_MAX, timehorizon - particles[i].time.back());
+                        particles[i].t += dt;
+                        particles[i].time.push_back(particles[i].t);
+                        particles[i].coordinates.push_back({particles[i].x, particles[i].y, particles[i].z});
+                        particles[i].domain.push_back(particles[i].domain.back());
+
+                    }
+
+                }
+
+            }
+
+        }
+
+    }
+
+    // move particle in the surface
+    void advecth(Particle &particle, double timehorizon)
+    {
+
+        double vx, vy, vmag, dt;
+        double pi;
+        double dx = parser.subgrid_dx;
+        bool is_inf;
+        int ix, iy, iz;
+
+        vx = environment.velocity_x[particle.cell_index];
+        vy = environment.velocity_y[particle.cell_index];
+
+        vmag = std::sqrt(vx * vx + vy * vy);
+
+        dt = CFL * dx / (vmag + 1.0e-10);
+        dt = std::min(dt, DT_MAX);
+        dt = std::min(dt, timehorizon - particle.time.back());
+
+        // check if infiltrates
+
+        is_inf = false;
+
+        double h = environment.water_depth[particle.cell_index];
+        double a = environment.cell_area[particle.cell_index];
+        double e = environment.exchange_flux[particle.cell_index] / a;
+
+        e = std::min(e, 0.0);
+        pi = std::abs(e) * dt / (h + WATERDEPTH_DRY);
+        pi = std::min(pi, 1.0);
+
+        // printf("%f %f %f\n", vx, vy, e);
+
+        if (dist(rng) < pi) {
+            is_inf = true;
+            //            printf("%f %f %f\n", vx, vy, e);
+        }
+
+        //
+
+        if (!is_inf) {
+            particle.advect(vx, vy, 0.0, dt);
+            iz = 0;
+        } else {
+            iz = 1;
+        }
+
+        // get the indices
+        ix = (int) floor((particle.x - environment.xll_corner) / dx);
+        iy = (int) floor((environment.yhl_corner - particle.y) / dx);
+
+        // update the particle
+        particle.time.push_back(particle.t);
+        particle.cell_index = environment.indices[ix][iy][iz];
+        particle.z = environment.elevations[ix][iy][iz];
+
+        int flag = SURFACE_FLAG;
+        if (particle.cell_index == -9999) {
+            flag = OUTLET_FLAG;
+            particle.z = particle.coordinates.back()[2];
+        } else if (iz > 0) {
+            flag = SUBSURFACE_FLAG;
+        }
+
+        particle.coordinates.push_back({particle.x,
+                                        particle.y,
+                                        particle.z});
+
+        particle.domain.push_back(flag);
+
+    }
+
+    // move particle in the subsurface
+    void advectv(Particle &particle, double timehorizon)
+    {
+
+        double vx, vy, vz, vmag, dt;
+        double dx = parser.subgrid_dx;
+        int ix, iy, iz;
+
+        //
+        // get the indices to determine the cell size
+        //
+
+        ix = (int) floor((particle.x - environment.xll_corner) / dx);
+        iy = (int) floor((environment.yhl_corner - particle.y) / dx);
+
+        double upper_horizon = environment.elevations[ix][iy][0];
+
+        for (size_t i = 0; i < environment.depths.size(); ++i) {
+
+            double dz = environment.depths[i];
+            upper_horizon -= dz;
+
+            if (particle.z > upper_horizon) {
+                iz = i + 1;
+                break;
+            }
+
+        }
+
+        //
+
+        vx = environment.darcy_velocity_x[particle.cell_index];
+        vy = environment.darcy_velocity_y[particle.cell_index];
+        vz = environment.darcy_velocity_z[particle.cell_index];
+        vmag = std::sqrt(vx * vx + vy * vy + vz * vz);
+
+        // cell size
+
+        dx = std::min(dx, environment.depths[iz-1]);
+        dt = CFL * dx / (vmag + 1.0e-10);
+        dt = std::min(dt, DT_MAX);
+        dt = std::min(dt, timehorizon - particle.time.back());
+
+        //
+
+        particle.advect(vx, vy, vz, dt);
+
+        dx = parser.subgrid_dx;
+
+        // get the new indices
+        ix = (int) floor((particle.x - environment.xll_corner) / dx);
+        iy = (int) floor((environment.yhl_corner - particle.y) / dx);
+
+        upper_horizon = environment.elevations[ix][iy][0];
+
+        if (particle.z > upper_horizon) {
+
+            iz = 0; // resurfaced (exfiltration)
+
+        } else {
+
+            for (size_t i = 0; i < environment.depths.size(); ++i) {
+
+                double dz = environment.depths[i];
+                upper_horizon -= dz;
+
+                if (particle.z > upper_horizon) {
+                    iz = i + 1;
+                    break;
+                }
+
+            }
+
+        }
+
+        // update the particle
+        int flag = SUBSURFACE_FLAG;
+
+        particle.time.push_back(particle.t);
+        particle.cell_index = environment.indices[ix][iy][iz];
+
+        if (iz == 0) {
+            particle.z = environment.elevations[ix][iy][iz];
+            flag = SURFACE_FLAG;
+        } else if (particle.cell_index == -9999) {
+            flag = OUTLET_FLAG;
+            particle.z = particle.coordinates.back()[2];
+        }
+
+        particle.coordinates.push_back({particle.x,
+                                        particle.y,
+                                        particle.z});
+
+        particle.domain.push_back(flag);
+
+    }
+
+
+    // // move particle in the subsurface
+    // void advectv()
+    // {
+
+    //     std::cout << "[ptl] advecting particles" << std::endl;
+
+    //     size_t i;
+    //     double dx = parser.subgrid_dx;
+
+    //     #pragma omp parallel
+    //     {
+
+    //         #pragma omp for private(i)
+    //         for (i = 0; i < particles.size(); ++i) {
+
+    //             while (particles[i].time.back() < timehorizon) {
+
+    //                 int flag = particles[i].domain.back();
+
+    //                 double vx, vy, vz, vmag, dt;
+
+    //                 vx = environment.velocity_x[particles[i].cell_index];
+    //                 vy = environment.velocity_y[particles[i].cell_index];
+
+    //                 vmag = std::sqrt(vx * vx + vy * vy);
+
+    //                 dt = CFL * dx / (vmag + 1.0e-10);
+    //                 dt = std::min(dt, DT_MAX);
+    //                 dt = std::min(dt, timehorizon - particles[i].time.back());
+
+    //                 particles[i].advect(vx, vy, 0.0, dt);
+
+    //                 // get the indices
+    //                 int ix, iy, iz;
+    //                 ix = (int) floor((particles[i].x - environment.xll_corner) / dx);
+    //                 iy = (int) floor((environment.yhl_corner - particles[i].y) / dx);
+    //                 iz = 0;
+
+    //                 // update the particle
+    //                 particles[i].time.push_back(particles[i].t);
+    //                 particles[i].cell_index = environment.indices[ix][iy][iz];
+    //                 particles[i].z = environment.elevations[ix][iy][iz];
+    //                 particles[i].coordinates.push_back({particles[i].x,
+    //                                                     particles[i].y,
+    //                                                     particles[i].z});
+    //                 particles[i].domain.push_back(flag);
+
+    //             }
+
+    //         }
+
+    //     }
+
+    // }
+
+
+
+private:
+
+    Parser parser;
+
+};
+
+#endif
diff --git a/subgrid.hpp b/subgrid.hpp
new file mode 100644
index 0000000..f050212
--- /dev/null
+++ b/subgrid.hpp
@@ -0,0 +1,99 @@
+// -*- mode: c++ style: k&r -*-
+
+#ifndef SUBGRID_HPP
+#define SUBGRID_HPP
+
+#include "header.hpp"
+
+// Content
+// =======
+//
+// 1. Subgrid
+// ---
+// 1.1 init
+//
+
+class Subgrid
+{
+
+public:
+
+    int nx;
+    int ny;
+    int nz;
+
+    index3d indices;   // three-dimensional array holding cell indices
+    grid3d elevations; // three-dimensional array holding elevations
+
+    std::vector<double> dzs;
+
+    Subgrid()
+    {
+        for (size_t i = 0; i < DIM; i ++) {
+            min[i] = 0.0;
+            max[i] = 0.0;
+        }
+    }
+
+    Subgrid(Parser &p_value, point &min_value, point &max_value, std::vector<double> &depths)
+    {
+        parser = p_value;
+
+        for (size_t i = 0; i < DIM; i ++) {
+            min[i] = min_value[i];
+            max[i] = max_value[i];
+        }
+
+        size_t ndepths = depths.size();
+        dzs.resize(ndepths);
+        for (size_t i = 0; i < ndepths; i ++)
+            dzs[i] = depths[i];
+
+        init();
+
+    }
+
+private:
+
+    Parser parser;
+    point min;
+    point max;
+
+    void init()
+    {
+
+        nx = (int) ceil((max[0] - min[0]) / parser.subgrid_dx);
+        ny = (int) ceil((max[1] - min[1]) / parser.subgrid_dx);
+        nz = dzs.size();
+
+#if VERBOSE
+        std::cout << "[sgd] vertical_res: { ";
+        for (size_t i = 0; i < dzs.size(); ++i) {
+            std::cout << dzs[i] << " ";
+        }
+        std::cout << "}" << std::endl;
+#endif
+
+        // initialise the grid and cell map
+        indices.resize(nx);
+        elevations.resize(nx);
+        for (int ix = 0; ix < nx; ++ix) {
+            indices[ix].resize(ny);
+            elevations[ix].resize(ny);
+            for (int iy = 0; iy < ny; ++iy) {
+                indices[ix][iy].resize(nz + 1);
+                elevations[ix][iy].resize(nz + 1);
+                std::fill(indices[ix][iy].begin(), indices[ix][iy].end(), -9999);
+                std::fill(elevations[ix][iy].begin(), elevations[ix][iy].end(), -9999.0);
+            }
+        }
+
+        std::cout << "[sgd] subgrid";
+        std::cout << " {" << indices.size() << ", " << indices[0].size() << ", " << indices[0][0].size() << "} ";
+        std::cout << "initialised with " << nx * ny * nz << " points" << std::endl;
+
+    }
+
+};
+
+#endif
diff --git a/tooling.hpp b/tooling.hpp
new file mode 100644
index 0000000..fc3aea2
--- /dev/null
+++ b/tooling.hpp
@@ -0,0 +1,102 @@
+// -*- mode: c++ -*-
+
+#ifndef TOOLING_HPP
+#define TOOLING_HPP
+
+#include <vector>
+#include <random>
+#include <algorithm>
+
+#define FUZZ 1.0e-12
+
+// Content
+// =======
+//
+// 1. zip
+// 2. unzip
+// 3. get_distance2d
+// 4. is_intersecting2d
+// 5. is_inside2d
+// 6. sgn
+//
+
+template <typename A, typename B>
+void zip(const std::vector<A> &a, const std::vector<B> &b, std::vector<std::pair<A,B>> &zipped)
+{
+    for (size_t i = 0; i < a.size(); ++i)
+        zipped.push_back(std::make_pair(a[i], b[i]));
+}
+
+template <typename A, typename B>
+void unzip(std::vector<A> &a, std::vector<B> &b, const std::vector<std::pair<A,B>> &zipped)
+{
+    for (size_t i = 0; i < a.size(); ++i) {
+        a[i] = zipped[i].first;
+        b[i] = zipped[i].second;
+    }
+}
+
+double get_distance2d(const point &p1, const point &p2)
+{
+    double distance = 0.0;
+    for (size_t i = 0; i < 2; i ++)
+        distance += (p1[i] - p2[i]) * (p1[i] - p2[i]);
+    return sqrt(distance);
+}
+
+bool is_intersecting2d(const point &p1, const point &p2, const point &p3, const point &p4)
+{
+
+    // values to determine t
+    double dx13 = p1[0] - p3[0];
+    double dx34 = p3[0] - p4[0];
+    double dx12 = p1[0] - p2[0];
+
+    double dy34 = p3[1] - p4[1];
+    double dy13 = p1[1] - p3[1];
+    double dy12 = p1[1] - p2[1];
+
+    // additional values to determine u
+    double dx21 = p2[0] - p1[0];
+    double dy21 = p2[1] - p1[1];
+
+    double det = (dx12 * dy34 - dy12 * dx34);
+    double t = (dx13 * dy34 - dy13 * dx34) / det;
+    double u = (dx21 * dy13 - dy21 * dx13) / det;
+
+    if ((t >= -FUZZ && t <= (1.0 + FUZZ)) && (u >= -FUZZ && u <= (1.0 + FUZZ)))
+        return true;
+
+    return false;
+
+}
+
+// draw a line from point of interest to the centroid of the polygon
+// if the line intersects any of the edges, the point is outside of the polygon
+// ! assumes convex polygons
+bool is_inside2d(const point &p, const point &centroid, const std::vector<point> &polygon)
+{
+
+    if (is_intersecting2d(p, centroid, polygon[polygon.size()-1], polygon[0]))
+        return false;
+
+    for (size_t i = 0; i < polygon.size() - 1; i ++) {
+        if (is_intersecting2d(p, centroid, polygon[i], polygon[i+1]))
+            return false;
+    }
+
+    return true;
+
+}
+
+template <typename T> int sgn(T val)
+{
+    return (T(0) < val) - (val < T(0));
+}
+
+double check_zero(double val)
+{
+    return (std::abs(val) <= ZERO_VELOCITY) ? 0.0 : val;
+}
+
+#endif
diff --git a/tools/extract-groups.sh b/tools/extract-groups.sh
new file mode 100755
index 0000000..5923f17
--- /dev/null
+++ b/tools/extract-groups.sh
@@ -0,0 +1,5 @@
+#!/usr/bin/env bash
+
+groupname=$2
+
+h5ls -r $1 | grep "${groupname}" | sed 's/\/.*\///' | sed 's/ .*Dataset.*$//' | sed '/Group$/d' | sort -n
diff --git a/tools/extract-times.sh b/tools/extract-times.sh
new file mode 100755
index 0000000..880e0d0
--- /dev/null
+++ b/tools/extract-times.sh
@@ -0,0 +1,10 @@
+#!/usr/bin/env bash
+
+groupname=$2
+
+IFS=$'\n' read -d '' -r -a groups
+
+for group in "${groups[@]}"
+do
+    h5dump -a "${groupname}/${group}/Time" $1 | grep "(0)" | sed 's/^ .*(0): //'
+done
diff --git a/tools/extract.sh b/tools/extract.sh
new file mode 100644
index 0000000..7222f83
--- /dev/null
+++ b/tools/extract.sh
@@ -0,0 +1,8 @@
+#!/usr/bin/env bash
+
+filename=$1
+groupname=$2
+
+./extract-groups.sh ${filename} ${groupname} > _tmpgroups.txt
+cat _tmpgroups.txt | ./extract-times.sh ${filename} ${groupname} > _tmptimes.txt
+cat _tmptimes.txt | awk '{ print $1 * 360000.0 * 24 }'
diff --git a/tools/merge-results.sh b/tools/merge-results.sh
deleted file mode 100755
index 648017f..0000000
--- a/tools/merge-results.sh
+++ /dev/null
@@ -1,11 +0,0 @@
-#!/usr/bin/env bash
-
-echo "[ok] remove existing aggregation file"
-rm -f ../results/agent-total.csv
-echo "[ok] cat everything into agent-total.csv"
-cat ../results/*.csv > ../results/agent-total.csv
-echo "[ok] cleaning up the file"
-gsed -i '/^time/d' ../results/agent-total.csv
-echo "[ok] inserting header"
-gsed -i '1i time (s), x (m), y (m), z (m), domain' ../results/agent-total.csv
-echo "[ok] done"
diff --git a/tools/merger.sh b/tools/merger.sh
new file mode 100644
index 0000000..a97a95d
--- /dev/null
+++ b/tools/merger.sh
@@ -0,0 +1,6 @@
+#!/usr/bin/bash
+
+for file in output/*.csv
+do
+    cat ${file} >> merge.txt
+done
diff --git a/tools/plot-paths-integrated.gp b/tools/plot-paths-integrated.gp
new file mode 100644
index 0000000..8fe4a06
--- /dev/null
+++ b/tools/plot-paths-integrated.gp
@@ -0,0 +1,27 @@
+#!/usr/bin/env gnuplot
+
+set term png font "hack"
+set output "output/figure.png"
+
+set border linecolor rgbcolor "white"
+set key textcolor rgbcolor "white"
+
+set obj 1 rectangle behind from screen 0,0 to screen 1,1
+set obj 1 fillstyle solid 1.0 fillcolor rgbcolor "black"
+
+set title "ptracker" textcolor rgbcolor "white"
+set xlabel "(km)" textcolor rgbcolor "white"
+set ylabel "(km)" textcolor rgbcolor "white"
+set zlabel "z (km)" textcolor rgbcolor "white"
+
+set xtics 1
+set ytics 1
+
+set size ratio -1
+
+splot \
+	for [i = 100:999] "output/p-000".i.".csv" \
+	u ($3/1000):($4/1000):($5/1000)\
+	with lines\
+	notitle
+
diff --git a/tools/plot-paths.gp b/tools/plot-paths.gp
new file mode 100644
index 0000000..6bc2fb1
--- /dev/null
+++ b/tools/plot-paths.gp
@@ -0,0 +1,25 @@
+#!/usr/bin/env gnuplot
+
+set term png font "hack"
+set output "output/figure.png"
+
+set border linecolor rgbcolor "white"
+set key textcolor rgbcolor "white"
+
+set obj 1 rectangle behind from screen 0,0 to screen 1,1
+set obj 1 fillstyle solid 1.0 fillcolor rgbcolor "black"
+
+set title "ptracker" textcolor rgbcolor "white"
+set xlabel "easting (km)" textcolor rgbcolor "white"
+set ylabel "northing (km)" textcolor rgbcolor "white"
+set xtics 1
+set ytics 1
+
+set size ratio -1
+
+plot \
+	for [i = 100:999] "output/p-000".i.".csv" \
+	u ($3/1000):($4/1000)\
+	with lines\
+	notitle
+