diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2dd13529b..4a9843ff2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -344,9 +344,6 @@ add_library(libopenmc SHARED
   src/photon_physics.F90
   src/physics_common.F90
   src/physics.F90
-  src/plot.F90
-  src/plot_header.F90
-  src/progress_header.F90
   src/pugixml/pugixml_f.F90
   src/random_lcg.F90
   src/reaction_header.F90
@@ -414,6 +411,7 @@ add_library(libopenmc SHARED
   src/physics_mg.cpp
   src/plot.cpp
   src/position.cpp
+  src/progress_bar.cpp
   src/pugixml/pugixml_c.cpp
   src/random_lcg.cpp
   src/reaction.cpp
diff --git a/include/openmc/output.h b/include/openmc/output.h
index e396c07d8..6d6d6ba0f 100644
--- a/include/openmc/output.h
+++ b/include/openmc/output.h
@@ -4,6 +4,7 @@
 #ifndef OPENMC_OUTPUT_H
 #define OPENMC_OUTPUT_H
 
+#include <string>
 
 namespace openmc {
 
@@ -16,6 +17,21 @@ namespace openmc {
 
 void header(const char* msg, int level);
 
+
+//==============================================================================
+//! Retrieve a time stamp
+//!
+//! \return current time stamp (format: "yyyy-mm-dd hh:mm:ss")
+//==============================================================================
+
+std::string time_stamp();
+
+//==============================================================================
+//! Display plot information
+//==============================================================================
+
+extern "C" void print_plot();
+
 //==============================================================================
 //! Display information regarding cell overlap checking.
 //==============================================================================
diff --git a/include/openmc/plot.h b/include/openmc/plot.h
index 2048189f7..64c092ce8 100644
--- a/include/openmc/plot.h
+++ b/include/openmc/plot.h
@@ -1,15 +1,181 @@
 #ifndef OPENMC_PLOT_H
 #define OPENMC_PLOT_H
 
+#include <unordered_map>
+#include <sstream>
+
+#include "xtensor/xarray.hpp"
+
 #include "hdf5.h"
+#include "openmc/position.h"
+#include "openmc/constants.h"
+#include "openmc/particle.h"
+#include "openmc/xml_interface.h"
 
 namespace openmc {
 
-extern "C" void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace,
+//===============================================================================
+// Global variables
+//===============================================================================
+
+extern int PLOT_LEVEL_LOWEST; //!< lower bound on plot universe level
+
+extern std::unordered_map<int, int> plot_map; //!< map of plot ids to index
+
+extern "C" int32_t n_plots; //!< number of plots in openmc run
+
+class Plot;
+extern std::vector<Plot> plots; //!< Plot instance container
+
+//===============================================================================
+// RGBColor holds color information for plotted objects
+//===============================================================================
+
+struct RGBColor {
+  //Constructors
+  RGBColor() : red(0), green(0), blue(0) { };
+  RGBColor(const int v[3]) : red(v[0]), green(v[1]), blue(v[2]) { };
+  RGBColor(int r, int g, int b) : red(r), green(g), blue(b) { };
+
+  RGBColor(const std::vector<int> &v) {
+    if (v.size() != 3) {
+      throw std::out_of_range("Incorrect vector size for RGBColor.");
+    }
+    red = v[0];
+    green = v[1];
+    blue = v[2];
+  }
+
+  // Members
+  uint8_t red, green, blue;
+};
+
+typedef xt::xtensor<RGBColor, 2> ImageData;
+
+enum class PlotType {
+  slice = 1,
+  voxel = 2
+};
+
+enum class PlotBasis {
+  xy = 1,
+  xz = 2,
+  yz = 3
+};
+
+enum class PlotColorBy {
+  cells = 1,
+  mats = 2
+};
+
+//===============================================================================
+// Plot class
+//===============================================================================
+
+class Plot
+{
+
+public:
+  // Constructor
+  Plot(pugi::xml_node plot);
+
+  // Methods
+private:
+  void set_id(pugi::xml_node plot_node);
+  void set_type(pugi::xml_node plot_node);
+  void set_output_path(pugi::xml_node plot_node);
+  void set_bg_color(pugi::xml_node plot_node);
+  void set_basis(pugi::xml_node plot_node);
+  void set_origin(pugi::xml_node plot_node);
+  void set_width(pugi::xml_node plot_node);
+  void set_universe(pugi::xml_node plot_node);
+  void set_default_colors(pugi::xml_node plot_node);
+  void set_user_colors(pugi::xml_node plot_node);
+  void set_meshlines(pugi::xml_node plot_node);
+  void set_mask(pugi::xml_node plot_node);
+
+  // Members
+public:
+  int id_; //!< Plot ID
+  PlotType type_; //!< Plot type (Slice/Voxel)
+  PlotColorBy color_by_; //!< Plot coloring (cell/material)
+  Position origin_; //!< Plot origin in geometry
+  Position width_; //!< Plot width in geometry
+  PlotBasis basis_; //!< Plot basis (XY/XZ/YZ)
+  std::array<int, 3> pixels_; //!< Plot size in pixels
+  int meshlines_width_; //!< Width of lines added to the plot
+  int level_; //!< Plot universe level
+  int index_meshlines_mesh_; //!< Index of the mesh to draw on the plot
+  RGBColor meshlines_color_; //!< Color of meshlines on the plot
+  RGBColor not_found_; //!< Plot background color
+  std::vector<RGBColor> colors_; //!< Plot colors
+  std::string path_plot_; //!< Plot output filename
+};
+
+//===============================================================================
+// Non-member functions
+//===============================================================================
+
+//! Add mesh lines to image data of a plot object
+//! \param[in] plot object
+//! \param[out] image data associated with the plot object
+void draw_mesh_lines(Plot pl, ImageData& data);
+
+//! Write a ppm image to file using a plot object's image data
+//! \param[in] plot object
+//! \param[out] image data associated with the plot object
+void output_ppm(Plot pl,
+                const ImageData& data);
+
+//! Get the rgb color for a given particle position in a plot
+//! \param[in] particle with position for current pixel
+//! \param[in] plot object
+//! \param[out] rgb color
+//! \param[out] cell or material id for particle position
+void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id);
+
+//! Initialize a voxel file
+//! \param[in] id of an open hdf5 file
+//! \param[in] dimensions of the voxel file (dx, dy, dz)
+//! \param[out] dataspace pointer to voxel data
+//! \param[out] dataset pointer to voxesl data
+//! \param[out] pointer to memory space of voxel data
+void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace,
                            hid_t* dset, hid_t* memspace);
-extern "C" void voxel_write_slice(int x, hid_t dspace, hid_t dset,
+
+//! Write a section of the voxel data to hdf5
+//! \param[in] voxel slice
+//! \param[out] dataspace pointer to voxel data
+//! \param[out] dataset pointer to voxesl data
+//! \param[out] pointer to data to write
+void voxel_write_slice(int x, hid_t dspace, hid_t dset,
                                   hid_t memspace, void* buf);
-extern "C" void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace);
+//! Close voxel file entities
+//! \param[in] data space to close
+//! \param[in] dataset to close
+//! \param[in] memory space to close
+void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace);
+
+//===============================================================================
+// External functions
+//===============================================================================
+
+//! Read plots from plots.xml node
+//! \param[in] plot node of plots.xml
+extern "C" void read_plots(pugi::xml_node* plot_node);
+
+//! Create a ppm image for a plot object
+//! \param[in] plot object
+void create_ppm(Plot pl);
+
+//! Create an hdf5 voxel file for a plot object
+//! \param[in] plot object
+void create_voxel(Plot pl);
+
+//! Create a randomly generated RGB color
+//! \return RGBColor with random value
+RGBColor random_color();
+
 
 } // namespace openmc
 #endif // OPENMC_PLOT_H
diff --git a/include/openmc/progress_bar.h b/include/openmc/progress_bar.h
new file mode 100644
index 000000000..18a3eea47
--- /dev/null
+++ b/include/openmc/progress_bar.h
@@ -0,0 +1,22 @@
+#ifndef OPENMC_PROGRESSBAR_H
+#define OPENMC_PROGRESSBAR_H
+
+#include <string>
+
+class ProgressBar {
+
+public:  
+  // Constructor
+  ProgressBar();
+
+  void set_value(double val);
+  
+private:
+  std::string bar;
+  char bar_old[72] = "???% |                                                                |";
+  
+};
+
+
+#endif // OPENMC_PROGRESSBAR_H
+  
diff --git a/include/openmc/settings.h b/include/openmc/settings.h
index aa9f9c056..47ca0f078 100644
--- a/include/openmc/settings.h
+++ b/include/openmc/settings.h
@@ -61,7 +61,8 @@ extern std::string path_statepoint;       //!< path to a statepoint file
 
 extern "C" int32_t index_entropy_mesh;  //!< Index of entropy mesh in global mesh array
 extern "C" int32_t index_ufs_mesh;      //!< Index of UFS mesh in global mesh array
-
+extern "C" int32_t index_cmfd_mesh;     //!< Index of CMFD mesh in global mesh array
+ 
 extern "C" int32_t n_batches;      //!< number of (inactive+active) batches
 extern "C" int32_t n_inactive;     //!< number of inactive batches
 extern "C" int32_t gen_per_batch;  //!< number of generations per batch
diff --git a/include/openmc/string_functions.h b/include/openmc/string_functions.h
index a3da179a0..a653d1e77 100644
--- a/include/openmc/string_functions.h
+++ b/include/openmc/string_functions.h
@@ -14,5 +14,7 @@ char* strtrim(char* c_str);
 
 void to_lower(std::string& str);
 
+int word_count(std::string const& str);
+
 } // namespace openmc
 #endif // STRING_FUNCTIONS_H
diff --git a/include/openmc/xml_interface.h b/include/openmc/xml_interface.h
index 85db26faa..6ab8c148c 100644
--- a/include/openmc/xml_interface.h
+++ b/include/openmc/xml_interface.h
@@ -10,7 +10,6 @@
 #include "xtensor/xarray.hpp"
 #include "xtensor/xadapt.hpp"
 
-
 namespace openmc {
 
 inline bool
diff --git a/src/api.F90 b/src/api.F90
index c5ab06f59..bb99180cc 100644
--- a/src/api.F90
+++ b/src/api.F90
@@ -63,7 +63,6 @@ contains
     use geometry_header
     use material_header
     use photon_header
-    use plot_header
     use sab_header
     use settings
     use simulation_header
@@ -85,7 +84,6 @@ contains
     call free_memory_geometry()
     call free_memory_surfaces()
     call free_memory_material()
-    call free_memory_plot()
     call free_memory_volume()
     call free_memory_simulation()
     call free_memory_nuclide()
diff --git a/src/cmfd_execute.cpp b/src/cmfd_execute.cpp
index aa9907a5b..46371601b 100644
--- a/src/cmfd_execute.cpp
+++ b/src/cmfd_execute.cpp
@@ -9,10 +9,10 @@
 #include "openmc/mesh.h"
 #include "openmc/message_passing.h"
 #include "openmc/simulation.h"
+#include "openmc/settings.h"
 
 namespace openmc {
 
-extern "C" int index_cmfd_mesh;
 
 extern "C" void
 cmfd_populate_sourcecounts(int n_energy, const double* energies,
@@ -24,7 +24,7 @@ cmfd_populate_sourcecounts(int n_energy, const double* energies,
   openmc_source_bank(&source_bank, &n);
 
   // Get source counts in each mesh bin / energy bin
-  auto& m = meshes.at(index_cmfd_mesh);
+  auto& m = meshes.at(settings::index_cmfd_mesh);
   xt::xarray<double> counts = m->count_sites(simulation::work, source_bank, n_energy, energies, outside);
 
   // Copy data from the xarray into the source counts array
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 466bff9e7..74ca51223 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -23,9 +23,8 @@ module input_xml
   use mgxs_interface
   use nuclide_header
   use multipole_header
-  use output,           only: title, header, print_plot
+  use output,           only: title, header
   use photon_header
-  use plot_header
   use random_lcg,       only: prn
   use surface_header
   use set_header,       only: SetChar, SetInt
@@ -97,6 +96,14 @@ module input_xml
       integer(C_INT)                        :: n
     end function maximum_levels
 
+    subroutine read_plots(node_ptr) bind(C)
+      import C_PTR
+      type(C_PTR) :: node_ptr
+    end subroutine read_plots
+
+    subroutine print_plot() bind(C)
+    end subroutine print_plot
+
     subroutine set_particle_energy_bounds(particle, E_min, E_max) bind(C)
       import C_INT, C_DOUBLE
       integer(C_INT), value :: particle
@@ -2030,26 +2037,10 @@ contains
 
   subroutine read_plots_xml()
 
-    integer :: i, j
-    integer :: n_cols, col_id, n_comp, n_masks, n_meshlines
-    integer :: meshid
-    integer(C_INT) :: err, idx
-    integer, allocatable :: iarray(:)
     logical :: file_exists              ! does plots.xml file exist?
     character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml
-    character(MAX_LINE_LEN) :: temp_str
-    character(MAX_WORD_LEN) :: meshtype
-    type(ObjectPlot), pointer :: pl => null()
     type(XMLDocument) :: doc
     type(XMLNode) :: root
-    type(XMLNode) :: node_plot
-    type(XMLNode) :: node_col
-    type(XMLNode) :: node_mask
-    type(XMLNode) :: node_meshlines
-    type(XMLNode), allocatable :: node_plot_list(:)
-    type(XMLNode), allocatable :: node_col_list(:)
-    type(XMLNode), allocatable :: node_mask_list(:)
-    type(XMLNode), allocatable :: node_meshline_list(:)
 
     ! Check if plots.xml exists
     filename = trim(path_input) // "plots.xml"
@@ -2066,425 +2057,7 @@ contains
     call doc % load_file(filename)
     root = doc % document_element()
 
-    ! Get list pointer to XML <plot>
-    call get_node_list(root, "plot", node_plot_list)
-
-    ! Allocate plots array
-    n_plots = size(node_plot_list)
-    allocate(plots(n_plots))
-
-    READ_PLOTS: do i = 1, n_plots
-      pl => plots(i)
-
-      ! Get pointer to plot XML node
-      node_plot = node_plot_list(i)
-
-      ! Copy data into plots
-      if (check_for_node(node_plot, "id")) then
-        call get_node_value(node_plot, "id", pl % id)
-      else
-        call fatal_error("Must specify plot id in plots XML file.")
-      end if
-
-      ! Check to make sure 'id' hasn't been used
-      if (plot_dict % has(pl % id)) then
-        call fatal_error("Two or more plots use the same unique ID: " &
-             // to_str(pl % id))
-      end if
-
-      ! Copy plot type
-      temp_str = 'slice'
-      if (check_for_node(node_plot, "type")) &
-           call get_node_value(node_plot, "type", temp_str)
-      temp_str = to_lower(temp_str)
-      select case (trim(temp_str))
-      case ("slice")
-        pl % type = PLOT_TYPE_SLICE
-      case ("voxel")
-        pl % type = PLOT_TYPE_VOXEL
-      case default
-        call fatal_error("Unsupported plot type '" // trim(temp_str) &
-             // "' in plot " // trim(to_str(pl % id)))
-      end select
-
-      ! Set output file path
-      filename = "plot_" // trim(to_str(pl % id))
-      if (check_for_node(node_plot, "filename")) &
-           call get_node_value(node_plot, "filename", filename)
-      select case (pl % type)
-      case (PLOT_TYPE_SLICE)
-        pl % path_plot = trim(path_input) // trim(filename) // ".ppm"
-      case (PLOT_TYPE_VOXEL)
-        pl % path_plot = trim(path_input) // trim(filename) // ".h5"
-      end select
-
-      ! Copy plot pixel size
-      if (pl % type == PLOT_TYPE_SLICE) then
-        if (node_word_count(node_plot, "pixels") == 2) then
-          call get_node_array(node_plot, "pixels", pl % pixels(1:2))
-        else
-          call fatal_error("<pixels> must be length 2 in slice plot " &
-               // trim(to_str(pl % id)))
-        end if
-      else if (pl % type == PLOT_TYPE_VOXEL) then
-        if (node_word_count(node_plot, "pixels") == 3) then
-          call get_node_array(node_plot, "pixels", pl % pixels(1:3))
-        else
-          call fatal_error("<pixels> must be length 3 in voxel plot " &
-               // trim(to_str(pl % id)))
-        end if
-      end if
-
-      ! Copy plot background color
-      if (check_for_node(node_plot, "background")) then
-        if (pl % type == PLOT_TYPE_VOXEL) then
-          if (master) call warning("Background color ignored in voxel plot " &
-               // trim(to_str(pl % id)))
-        end if
-        if (node_word_count(node_plot, "background") == 3) then
-          call get_node_array(node_plot, "background", pl % not_found % rgb)
-        else
-          call fatal_error("Bad background RGB in plot " &
-               // trim(to_str(pl % id)))
-        end if
-      else
-        pl % not_found % rgb = (/ 255, 255, 255 /)
-      end if
-
-      ! Copy plot basis
-      if (pl % type == PLOT_TYPE_SLICE) then
-        temp_str = 'xy'
-        if (check_for_node(node_plot, "basis")) &
-             call get_node_value(node_plot, "basis", temp_str)
-        temp_str = to_lower(temp_str)
-        select case (trim(temp_str))
-        case ("xy")
-          pl % basis = PLOT_BASIS_XY
-        case ("xz")
-          pl % basis = PLOT_BASIS_XZ
-        case ("yz")
-          pl % basis = PLOT_BASIS_YZ
-        case default
-          call fatal_error("Unsupported plot basis '" // trim(temp_str) &
-               // "' in plot " // trim(to_str(pl % id)))
-        end select
-      end if
-
-      ! Copy plotting origin
-      if (node_word_count(node_plot, "origin") == 3) then
-        call get_node_array(node_plot, "origin", pl % origin)
-      else
-        call fatal_error("Origin must be length 3 in plot " &
-             // trim(to_str(pl % id)))
-      end if
-
-      ! Copy plotting width
-      if (pl % type == PLOT_TYPE_SLICE) then
-        if (node_word_count(node_plot, "width") == 2) then
-          call get_node_array(node_plot, "width", pl % width(1:2))
-        else
-          call fatal_error("<width> must be length 2 in slice plot " &
-               // trim(to_str(pl % id)))
-        end if
-      else if (pl % type == PLOT_TYPE_VOXEL) then
-        if (node_word_count(node_plot, "width") == 3) then
-          call get_node_array(node_plot, "width", pl % width(1:3))
-        else
-          call fatal_error("<width> must be length 3 in voxel plot " &
-               // trim(to_str(pl % id)))
-        end if
-      end if
-
-      ! Copy plot cell universe level
-      if (check_for_node(node_plot, "level")) then
-        call get_node_value(node_plot, "level", pl % level)
-
-        if (pl % level < 0) then
-          call fatal_error("Bad universe level in plot " &
-               // trim(to_str(pl % id)))
-        end if
-      else
-        pl % level = PLOT_LEVEL_LOWEST
-      end if
-
-      ! Copy plot color type and initialize all colors randomly
-      temp_str = "cell"
-      if (check_for_node(node_plot, "color_by")) &
-           call get_node_value(node_plot, "color_by", temp_str)
-      temp_str = to_lower(temp_str)
-      select case (trim(temp_str))
-      case ("cell")
-
-        pl % color_by = PLOT_COLOR_CELLS
-        allocate(pl % colors(n_cells))
-        do j = 1, n_cells
-          pl % colors(j) % rgb(1) = int(prn()*255)
-          pl % colors(j) % rgb(2) = int(prn()*255)
-          pl % colors(j) % rgb(3) = int(prn()*255)
-        end do
-
-      case ("material")
-
-        pl % color_by = PLOT_COLOR_MATS
-        allocate(pl % colors(n_materials))
-        do j = 1, n_materials
-          pl % colors(j) % rgb(1) = int(prn()*255)
-          pl % colors(j) % rgb(2) = int(prn()*255)
-          pl % colors(j) % rgb(3) = int(prn()*255)
-        end do
-
-      case default
-        call fatal_error("Unsupported plot color type '" // trim(temp_str) &
-             // "' in plot " // trim(to_str(pl % id)))
-      end select
-
-      ! Get the number of <color> nodes and get a list of them
-      call get_node_list(node_plot, "color", node_col_list)
-      n_cols = size(node_col_list)
-
-      ! Copy user specified colors
-      if (n_cols /= 0) then
-
-        if (pl % type == PLOT_TYPE_VOXEL) then
-          if (master) call warning("Color specifications ignored in voxel &
-               &plot " // trim(to_str(pl % id)))
-        end if
-
-        do j = 1, n_cols
-
-          ! Get pointer to color spec XML node
-          node_col = node_col_list(j)
-
-          ! Check and make sure 3 values are specified for RGB
-          if (node_word_count(node_col, "rgb") /= 3) then
-            call fatal_error("Bad RGB in plot " &
-                 // trim(to_str(pl % id)))
-          end if
-
-          ! Ensure that there is an id for this color specification
-          if (check_for_node(node_col, "id")) then
-            call get_node_value(node_col, "id", col_id)
-          else
-            call fatal_error("Must specify id for color specification in &
-                 &plot " // trim(to_str(pl % id)))
-          end if
-
-          ! Add RGB
-          if (pl % color_by == PLOT_COLOR_CELLS) then
-
-            if (cell_dict % has(col_id)) then
-              col_id = cell_dict % get(col_id)
-              call get_node_array(node_col, "rgb", pl % colors(col_id) % rgb)
-            else
-              call fatal_error("Could not find cell " // trim(to_str(col_id)) &
-                   // " specified in plot " // trim(to_str(pl % id)))
-            end if
-
-          else if (pl % color_by == PLOT_COLOR_MATS) then
-
-            if (material_dict % has(col_id)) then
-              col_id = material_dict % get(col_id)
-              call get_node_array(node_col, "rgb", pl % colors(col_id) % rgb)
-            else
-              call fatal_error("Could not find material " &
-                   // trim(to_str(col_id)) // " specified in plot " &
-                   // trim(to_str(pl % id)))
-            end if
-
-          end if
-        end do
-      end if
-
-      ! Deal with meshlines
-      call get_node_list(node_plot, "meshlines", node_meshline_list)
-      n_meshlines = size(node_meshline_list)
-      if (n_meshlines /= 0) then
-
-        if (pl % type == PLOT_TYPE_VOXEL) then
-          call warning("Meshlines ignored in voxel plot " &
-               // trim(to_str(pl % id)))
-        end if
-
-        select case(n_meshlines)
-          case (0)
-            ! Skip if no meshlines are specified
-          case (1)
-
-            ! Get pointer to meshlines
-            node_meshlines = node_meshline_list(1)
-
-            ! Check mesh type
-            if (check_for_node(node_meshlines, "meshtype")) then
-              call get_node_value(node_meshlines, "meshtype", meshtype)
-            else
-              call fatal_error("Must specify a meshtype for meshlines &
-                   &specification in plot " // trim(to_str(pl % id)))
-            end if
-
-            ! Ensure that there is a linewidth for this meshlines specification
-            if (check_for_node(node_meshlines, "linewidth")) then
-              call get_node_value(node_meshlines, "linewidth", &
-                   pl % meshlines_width)
-            else
-              call fatal_error("Must specify a linewidth for meshlines &
-                   &specification in plot " // trim(to_str(pl % id)))
-            end if
-
-            ! Check for color
-            if (check_for_node(node_meshlines, "color")) then
-
-              ! Check and make sure 3 values are specified for RGB
-              if (node_word_count(node_meshlines, "color") /= 3) then
-                call fatal_error("Bad RGB for meshlines color in plot " &
-                     // trim(to_str(pl % id)))
-              end if
-
-              call get_node_array(node_meshlines, "color", &
-                   pl % meshlines_color % rgb)
-            else
-
-              pl % meshlines_color % rgb = (/ 0, 0, 0 /)
-
-            end if
-
-            ! Set mesh based on type
-            select case (trim(meshtype))
-            case ('ufs')
-
-              if (index_ufs_mesh < 0) then
-                call fatal_error("No UFS mesh for meshlines on plot " &
-                     // trim(to_str(pl % id)))
-              end if
-
-              pl % index_meshlines_mesh = index_ufs_mesh
-
-            case ('cmfd')
-
-              if (.not. cmfd_run) then
-                call fatal_error("Need CMFD run to plot CMFD mesh for &
-                     &meshlines on plot " // trim(to_str(pl % id)))
-              end if
-
-              pl % index_meshlines_mesh = index_cmfd_mesh
-
-            case ('entropy')
-
-              if (index_entropy_mesh < 0) then
-                call fatal_error("No entropy mesh for meshlines on plot " &
-                     // trim(to_str(pl % id)))
-              end if
-
-              pl % index_meshlines_mesh = index_entropy_mesh
-
-            case ('tally')
-
-              ! Ensure that there is a mesh id if the type is tally
-              if (check_for_node(node_meshlines, "id")) then
-                call get_node_value(node_meshlines, "id", meshid)
-              else
-                call fatal_error("Must specify a mesh id for meshlines tally &
-                     &mesh specification in plot " // trim(to_str(pl % id)))
-              end if
-
-              ! Check if the specified tally mesh exists
-              err = openmc_get_mesh_index(meshid, idx)
-              if (err /= 0) then
-                call fatal_error("Could not find mesh " &
-                     // trim(to_str(meshid)) // " specified in meshlines for &
-                     &plot " // trim(to_str(pl % id)))
-              end if
-              pl % index_meshlines_mesh = idx
-
-            case default
-              call fatal_error("Invalid type for meshlines on plot " &
-                    // trim(to_str(pl % id)) // ": " // trim(meshtype))
-            end select
-
-          case default
-            call fatal_error("Mutliple meshlines specified in plot " &
-                 // trim(to_str(pl % id)))
-        end select
-
-      end if
-
-      ! Deal with masks
-      call get_node_list(node_plot, "mask", node_mask_list)
-      n_masks = size(node_mask_list)
-      if (n_masks /= 0) then
-
-        if (pl % type == PLOT_TYPE_VOXEL) then
-          if (master) call warning("Mask ignored in voxel plot " &
-               // trim(to_str(pl % id)))
-        end if
-
-        select case(n_masks)
-          case default
-            call fatal_error("Mutliple masks specified in plot " &
-                 // trim(to_str(pl % id)))
-          case (1)
-
-            ! Get pointer to mask
-            node_mask = node_mask_list(1)
-
-            ! Determine how many components there are and allocate
-            n_comp = 0
-            n_comp = node_word_count(node_mask, "components")
-            if (n_comp == 0) then
-              call fatal_error("Missing <components> in mask of plot " &
-                   // trim(to_str(pl % id)))
-            end if
-            allocate(iarray(n_comp))
-            call get_node_array(node_mask, "components", iarray)
-
-            ! First we need to change the user-specified identifiers to indices
-            ! in the cell and material arrays
-            do j=1, n_comp
-              col_id = iarray(j)
-
-              if (pl % color_by == PLOT_COLOR_CELLS) then
-
-                if (cell_dict % has(col_id)) then
-                  iarray(j) = cell_dict % get(col_id)
-                else
-                  call fatal_error("Could not find cell " &
-                       // trim(to_str(col_id)) // " specified in the mask in &
-                       &plot " // trim(to_str(pl % id)))
-                end if
-
-              else if (pl % color_by == PLOT_COLOR_MATS) then
-
-                if (material_dict % has(col_id)) then
-                  iarray(j) = material_dict % get(col_id)
-                else
-                  call fatal_error("Could not find material " &
-                       // trim(to_str(col_id)) // " specified in the mask in &
-                       &plot " // trim(to_str(pl % id)))
-                end if
-
-              end if
-            end do
-
-            ! Alter colors based on mask information
-            do j = 1, size(pl % colors)
-              if (.not. any(j == iarray)) then
-                if (check_for_node(node_mask, "background")) then
-                  call get_node_array(node_mask, "background", pl % colors(j) % rgb)
-                else
-                  pl % colors(j) % rgb(:) = [255, 255, 255]
-                end if
-              end if
-            end do
-
-            deallocate(iarray)
-
-        end select
-
-      end if
-
-      ! Add plot to dictionary
-      call plot_dict % set(pl % id, i)
-
-    end do READ_PLOTS
+    call read_plots(root % ptr)
 
     ! Close plots XML file
     call doc % clear()
diff --git a/src/output.F90 b/src/output.F90
index ef3d3259c..485d957a6 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -14,7 +14,6 @@ module output
   use mgxs_interface
   use nuclide_header
   use particle_header, only: LocalCoord, Particle
-  use plot_header
   use sab_header,      only: SAlphaBeta
   use settings
   use simulation_header
@@ -412,79 +411,6 @@ contains
 
   end subroutine print_batch_keff
 
-!===============================================================================
-! PRINT_PLOT displays selected options for plotting
-!===============================================================================
-
-  subroutine print_plot()
-
-    integer :: i ! loop index for plots
-    type(ObjectPlot), pointer :: pl
-
-    ! Display header for plotting
-    call header("PLOTTING SUMMARY", 5)
-
-    do i = 1, n_plots
-      pl => plots(i)
-
-      ! Plot id
-      write(ou,100) "Plot ID:", trim(to_str(pl % id))
-
-      ! Plot filename
-      write(ou,100) "Plot file:", trim(pl % path_plot)
-
-      ! Plot level
-      write(ou,100) "Universe depth:", trim(to_str(pl % level))
-
-      ! Plot type
-      if (pl % type == PLOT_TYPE_SLICE) then
-        write(ou,100) "Plot Type:", "Slice"
-      else if (pl % type == PLOT_TYPE_VOXEL) then
-        write(ou,100) "Plot Type:", "Voxel"
-      end if
-
-      ! Plot parameters
-      write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // &
-           " " // trim(to_str(pl % origin(2))) // " " // &
-           trim(to_str(pl % origin(3)))
-      if (pl % type == PLOT_TYPE_SLICE) then
-        write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
-             " " // trim(to_str(pl % width(2)))
-      else if (pl % type == PLOT_TYPE_VOXEL) then
-        write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
-             " " // trim(to_str(pl % width(2))) // &
-             " " // trim(to_str(pl % width(3)))
-      end if
-      if (pl % color_by == PLOT_COLOR_CELLS) then
-        write(ou,100) "Coloring:", "Cells"
-      else if (pl % color_by == PLOT_COLOR_MATS) then
-        write(ou,100) "Coloring:", "Materials"
-      end if
-      if (pl % type == PLOT_TYPE_SLICE) then
-        select case (pl % basis)
-        case (PLOT_BASIS_XY)
-          write(ou,100) "Basis:", "xy"
-        case (PLOT_BASIS_XZ)
-          write(ou,100) "Basis:", "xz"
-        case (PLOT_BASIS_YZ)
-          write(ou,100) "Basis:", "yz"
-        end select
-        write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // &
-             trim(to_str(pl % pixels(2)))
-      else if (pl % type == PLOT_TYPE_VOXEL) then
-        write(ou,100) "Voxels:", trim(to_str(pl % pixels(1))) // " " // &
-             trim(to_str(pl % pixels(2))) // " " // trim(to_str(pl % pixels(3)))
-      end if
-
-      write(ou,*)
-
-    end do
-
-    ! Format descriptor for columns
-100 format (1X,A,T25,A)
-
-  end subroutine print_plot
-
 !===============================================================================
 ! PRINT_RUNTIME displays the total time elapsed for the entire run, for
 ! initialization, for computation, and for intergeneration synchronization.
diff --git a/src/output.cpp b/src/output.cpp
index 6a3abeed6..79800e21a 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -5,13 +5,14 @@
 #include <iomanip>  // for setw
 #include <iostream>
 #include <sstream>
+#include <ctime>
 
 #include "openmc/cell.h"
 #include "openmc/geometry.h"
 #include "openmc/message_passing.h"
 #include "openmc/capi.h"
 #include "openmc/settings.h"
-
+#include "openmc/plot.h"
 
 namespace openmc {
 
@@ -35,13 +36,98 @@ header(const char* msg, int level) {
 
   // Print header based on verbosity level.
   if (settings::verbosity >= level) {
-    std::cout << out.str() << std::endl << std::endl;
+    std::cout << out.str() << "\n\n";
   }
 }
 
+std::string time_stamp()
+{
+  int base_year = 1990;
+  std::stringstream ts;
+  std::time_t t = std::time(0);   // get time now
+  std::tm* now = std::localtime(&t);
+  ts << now->tm_year + base_year << "-" << now->tm_mon
+     << "-" << now->tm_mday << " " << now->tm_hour
+     << ":" << now->tm_min << ":" << now->tm_sec;
+  return ts.str();
+}
+  
 //==============================================================================
 
-void print_overlap_check() {
+//===============================================================================
+// PRINT_PLOT displays selected options for plotting
+//===============================================================================
+
+void print_plot() {
+
+  header("PLOTTING SUMMARY", 5);
+
+  for (auto pl : plots) {
+    // Plot id
+    std::cout << "Plot ID: " << pl.id_ << "\n";
+    // Plot filename
+    std::cout << "Plot file: " << pl.path_plot_ << "\n";
+    // Plot level
+    std::cout << "Universe depth: " << pl.level_ << "\n";
+
+    // Plot type
+    if (PlotType::slice == pl.type_) {
+      std::cout << "Plot Type: Slice" << "\n";
+    } else if (PlotType::voxel == pl.type_) {
+      std::cout << "Plot Type: Voxel" << "\n";      
+    }
+
+    // Plot parameters
+    std::cout << "Origin: " << pl.origin_[0] << " "
+              << pl.origin_[1] << " "
+              << pl.origin_[2] << "\n";
+
+    if (PlotType::slice == pl.type_) {
+      std::cout << std::setprecision(4)
+                << "Width: "
+                << pl.width_[0] << " "
+                << pl.width_[1] << "\n";
+    } else if (PlotType::voxel == pl.type_) {
+      std::cout << std::setprecision(4)
+                << "Width: "
+                << pl.width_[0] << " "
+                << pl.width_[1] << " "
+                << pl.width_[2] << "\n";
+    }
+
+    if (PlotColorBy::cells == pl.color_by_) {
+      std::cout << "Coloring: Cells" << "\n";
+    } else if (PlotColorBy::mats == pl.color_by_) {
+      std::cout << "Coloring: Materials" << "\n";      
+    }
+    
+    if (PlotType::slice == pl.type_) {
+      switch(pl.basis_) {
+      case PlotBasis::xy:
+        std::cout <<  "Basis: XY" << "\n";
+        break;
+      case PlotBasis::xz:
+        std::cout <<  "Basis: XZ" << "\n";
+        break;
+      case PlotBasis::yz:
+        std::cout <<  "Basis: YZ" << "\n";
+        break;
+      }
+      std::cout << "Pixels: " << pl.pixels_[0] << " "
+                << pl.pixels_[1] << " " << "\n";
+    } else if (PlotType::voxel == pl.type_) {
+      std::cout << "Voxels: " << pl.pixels_[0] << " "
+                << pl.pixels_[1] << " "
+                << pl.pixels_[2] << "\n";
+    }
+
+    std::cout << "\n";
+    
+  }
+}
+  
+void
+print_overlap_check() {
 #ifdef OPENMC_MPI
   std::vector<int64_t> temp(overlap_check_count);
   int err = MPI_Reduce(temp.data(), overlap_check_count.data(),
diff --git a/src/plot.F90 b/src/plot.F90
deleted file mode 100644
index bfc8bc225..000000000
--- a/src/plot.F90
+++ /dev/null
@@ -1,456 +0,0 @@
-module plot
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants
-  use error,           only: fatal_error, write_message
-  use geometry,        only: find_cell, check_cell_overlap
-  use geometry_header, only: Cell, root_universe, cells
-  use hdf5_interface
-  use output,          only: time_stamp
-  use material_header, only: materials
-  use mesh_header,     only: meshes, RegularMesh
-  use particle_header
-  use plot_header
-  use progress_header, only: ProgressBar
-  use settings,        only: check_overlaps
-  use string,          only: to_str
-
-  implicit none
-  private
-
-  public :: openmc_plot_geometry
-
-  integer, parameter :: RED = 1
-  integer, parameter :: GREEN = 2
-  integer, parameter :: BLUE = 3
-
-contains
-
-!===============================================================================
-! RUN_PLOT controls the logic for making one or many plots
-!===============================================================================
-
-  function openmc_plot_geometry() result(err) bind(C)
-    integer(C_INT) :: err
-
-    integer :: i ! loop index for plots
-
-    do i = 1, n_plots
-      associate (pl => plots(i))
-        ! Display output message
-        call write_message("Processing plot " // trim(to_str(pl % id)) &
-             // ": " // trim(pl % path_plot) // " ...", 5)
-
-        if (pl % type == PLOT_TYPE_SLICE) then
-          ! create 2d image
-          call create_ppm(pl)
-        else if (pl % type == PLOT_TYPE_VOXEL) then
-          ! create dump for 3D silomesh utility script
-          call create_voxel(pl)
-        end if
-      end associate
-    end do
-
-    err = 0
-  end function openmc_plot_geometry
-
-!===============================================================================
-! POSITION_RGB computes the red/green/blue values for a given plot with the
-! current particle's position
-!===============================================================================
-
-  subroutine position_rgb(p, pl, rgb, id)
-    type(Particle), intent(inout) :: p
-    type(ObjectPlot), intent(in)  :: pl
-    integer, intent(out)          :: rgb(3)
-    integer, intent(out)          :: id
-
-    integer :: j
-    logical :: found_cell
-
-    p % n_coord = 1
-
-    call find_cell(p, found_cell)
-    j = p % n_coord
-    if (check_overlaps) call check_cell_overlap(p)
-
-    ! Set coordinate level if specified
-    if (pl % level >= 0) j = pl % level + 1
-
-    if (.not. found_cell) then
-      ! If no cell, revert to default color
-      rgb = pl % not_found % rgb
-      id = -1
-    else
-      if (pl % color_by == PLOT_COLOR_MATS) then
-        ! Assign color based on material
-        associate (c => cells(p % coord(j) % cell + 1))
-          if (c % type() == FILL_UNIVERSE) then
-            ! If we stopped on a middle universe level, treat as if not found
-            rgb = pl % not_found % rgb
-            id = -1
-          else if (p % material == MATERIAL_VOID) then
-            ! By default, color void cells white
-            rgb = 255
-            id = -1
-          else
-            rgb = pl % colors(p % material) % rgb
-            id = materials(p % material) % id()
-          end if
-        end associate
-      else if (pl % color_by == PLOT_COLOR_CELLS) then
-        ! Assign color based on cell
-        rgb = pl % colors(p % coord(j) % cell + 1) % rgb
-        id = cells(p % coord(j) % cell + 1) % id()
-      else
-        rgb = 0
-        id = -1
-      end if
-    end if
-
-  end subroutine position_rgb
-
-!===============================================================================
-! CREATE_PPM creates an image based on user input from a plots.xml <plot>
-! specification in the portable pixmap format (PPM)
-!===============================================================================
-
-  subroutine create_ppm(pl)
-    type(ObjectPlot), intent(in) :: pl
-
-    integer :: in_i
-    integer :: out_i
-    integer :: x, y      ! pixel location
-    integer :: rgb(3)    ! colors (red, green, blue) from 0-255
-    integer :: id
-    integer :: height, width
-    real(8) :: in_pixel
-    real(8) :: out_pixel
-    real(8) :: xyz(3)
-    integer, allocatable :: data(:,:,:)
-    type(Particle)    :: p
-
-    width = pl % pixels(1)
-    height = pl % pixels(2)
-
-    in_pixel  = pl % width(1)/dble(width)
-    out_pixel = pl % width(2)/dble(height)
-
-    ! Allocate and initialize results array
-    allocate(data(3, width, height))
-    data(:,:,:) = 0
-
-    if (pl % basis == PLOT_BASIS_XY) then
-      in_i  = 1
-      out_i = 2
-      xyz(1) = pl % origin(1) - pl % width(1) / TWO
-      xyz(2) = pl % origin(2) + pl % width(2) / TWO
-      xyz(3) = pl % origin(3)
-    else if (pl % basis == PLOT_BASIS_XZ) then
-      in_i  = 1
-      out_i = 3
-      xyz(1) = pl % origin(1) - pl % width(1) / TWO
-      xyz(2) = pl % origin(2)
-      xyz(3) = pl % origin(3) + pl % width(2) / TWO
-    else if (pl % basis == PLOT_BASIS_YZ) then
-      in_i  = 2
-      out_i = 3
-      xyz(1) = pl % origin(1)
-      xyz(2) = pl % origin(2) - pl % width(1) / TWO
-      xyz(3) = pl % origin(3) + pl % width(2) / TWO
-    end if
-
-    ! allocate and initialize particle
-    call particle_initialize(p)
-    p % coord(1) % xyz = xyz
-    p % coord(1) % uvw = [ HALF, HALF, HALF ]
-    p % coord(1) % universe = root_universe
-
-!$omp parallel do firstprivate(p) private(x, rgb, id) reduction(+ : data)
-    do y = 1, height
-      ! Set y coordinate
-      p % coord(1) % xyz(out_i) = xyz(out_i) - out_pixel*(y - 1)
-      do x = 1, width
-        ! Set x coordinate
-        p % coord(1) % xyz(in_i) = xyz(in_i) + in_pixel*(x - 1)
-
-        ! get pixel color
-        call position_rgb(p, pl, rgb, id)
-
-        ! Create a pixel at (x,y) with color (r,g,b)
-        data(:, x, y) = rgb
-      end do
-    end do
-!$omp end parallel do
-
-    ! Draw tally mesh boundaries on the image if requested
-    if (pl % index_meshlines_mesh >= 0) call draw_mesh_lines(pl, data)
-
-    ! Write out the ppm to a file
-    call output_ppm(pl, data)
-
-  end subroutine create_ppm
-
-!===============================================================================
-! DRAW_MESH_LINES draws mesh line boundaries on an image
-!===============================================================================
-
-  subroutine draw_mesh_lines(pl, data)
-    type(ObjectPlot), intent(in)    :: pl
-    integer,          intent(inout) :: data(:,:,:)
-
-    logical :: in_mesh
-    integer :: out_, in_  ! pixel location
-    integer :: rgb(3)     ! RGB color for meshlines pixels
-    integer :: outrange(2), inrange(2) ! range of pixel locations
-    integer :: i, j       ! loop indices
-    integer :: plus
-    integer :: ijk_ll(3)  ! mesh bin ijk indicies of plot lower left
-    integer :: ijk_ur(3)  ! mesh bin ijk indicies of plot upper right
-    integer :: outer, inner
-    real(8) :: frac
-    real(8) :: width(3)   ! real widths of the plot
-    real(8) :: xyz_ll_plot(3)  ! lower left xyz of plot image
-    real(8) :: xyz_ur_plot(3)  ! upper right xyz of plot image
-    real(8) :: xyz_ll(3)  ! lower left xyz
-    real(8) :: xyz_ur(3)  ! upper right xyz
-    type(RegularMesh) :: m
-
-    rgb(:) = pl % meshlines_color % rgb
-
-    select case (pl % basis)
-      case(PLOT_BASIS_XY)
-        outer = 1
-        inner = 2
-      case(PLOT_BASIS_XZ)
-        outer = 1
-        inner = 3
-      case(PLOT_BASIS_YZ)
-        outer = 2
-        inner = 3
-    end select
-
-    xyz_ll_plot = pl % origin
-    xyz_ur_plot = pl % origin
-
-    xyz_ll_plot(outer) = pl % origin(outer) - pl % width(1) / TWO
-    xyz_ll_plot(inner) = pl % origin(inner) - pl % width(2) / TWO
-    xyz_ur_plot(outer) = pl % origin(outer) + pl % width(1) / TWO
-    xyz_ur_plot(inner) = pl % origin(inner) + pl % width(2) / TWO
-
-    width = xyz_ur_plot - xyz_ll_plot
-
-    m = meshes(pl % index_meshlines_mesh)
-    call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension()), in_mesh)
-    call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension()), in_mesh)
-
-    ! sweep through all meshbins on this plane and draw borders
-    do i = ijk_ll(outer), ijk_ur(outer)
-      do j = ijk_ll(inner), ijk_ur(inner)
-        ! check if we're in the mesh for this ijk
-        if (i > 0 .and. i <= m % dimension(outer) .and. &
-              j > 0 .and. j <= m % dimension(inner)) then
-
-          ! get xyz's of lower left and upper right of this mesh cell
-          xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1)
-          xyz_ll(inner) = m % lower_left(inner) + m % width(inner) * (j - 1)
-          xyz_ur(outer) = m % lower_left(outer) + m % width(outer) * i
-          xyz_ur(inner) = m % lower_left(inner) + m % width(inner) * j
-
-          ! map the xyz ranges to pixel ranges
-
-          frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer)
-          outrange(1) = int(frac * real(pl % pixels(1), 8))
-          frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer)
-          outrange(2) = int(frac * real(pl % pixels(1), 8))
-
-          frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner)
-          inrange(1) = int((ONE - frac) * real(pl % pixels(2), 8))
-          frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner)
-          inrange(2) = int((ONE - frac) * real(pl % pixels(2), 8))
-
-          ! draw lines
-          do out_ = outrange(1), outrange(2)
-            do plus = 0, pl % meshlines_width
-              data(:, out_ + 1, inrange(1) + plus + 1) = rgb
-              data(:, out_ + 1, inrange(2) + plus + 1) = rgb
-              data(:, out_ + 1, inrange(1) - plus + 1) = rgb
-              data(:, out_ + 1, inrange(2) - plus + 1) = rgb
-            end do
-          end do
-          do in_ = inrange(1), inrange(2)
-            do plus = 0, pl % meshlines_width
-              data(:, outrange(1) + plus + 1, in_ + 1) = rgb
-              data(:, outrange(2) + plus + 1, in_ + 1) = rgb
-              data(:, outrange(1) - plus + 1, in_ + 1) = rgb
-              data(:, outrange(2) - plus + 1, in_ + 1) = rgb
-            end do
-          end do
-
-        end if
-      end do
-    end do
-
-  end subroutine draw_mesh_lines
-
-!===============================================================================
-! OUTPUT_PPM writes out a previously generated image to a PPM file
-!===============================================================================
-
-  subroutine output_ppm(pl, data)
-    type(ObjectPlot), intent(in) :: pl
-    integer,          intent(in) :: data(:,:,:)
-
-    integer :: y ! loop index for height
-    integer :: x ! loop index for width
-    integer :: unit_plot
-
-    ! Open PPM file for writing
-    open(NEWUNIT=unit_plot, FILE=pl % path_plot)
-
-    ! Write header
-    write(unit_plot, '(A2)') 'P6'
-    write(unit_plot, '(I0,'' '',I0)') pl % pixels(1), pl % pixels(2)
-    write(unit_plot, '(A)') '255'
-
-    ! Write color for each pixel
-    do y = 1, pl % pixels(2)
-      do x = 1, pl % pixels(1)
-        write(unit_plot, '(3A1)', advance='no') achar(data(RED, x, y)), &
-             achar(data(GREEN, x, y)), achar(data(BLUE, x, y))
-      end do
-    end do
-
-    ! Close plot file
-    close(UNIT=unit_plot)
-  end subroutine output_ppm
-
-!===============================================================================
-! CREATE_VOXEL outputs a binary file that can be input into silomesh for 3D
-! geometry visualization.  It works the same way as create_ppm by dragging a
-! particle across the geometry for the specified number of voxels. The first
-! 3 int(4)'s in the binary are the number of x, y, and z voxels.  The next 3
-! real(8)'s are the widths of the voxels in the x, y, and z directions. The next
-! 3 real(8)'s are the x, y, and z coordinates of the lower left point. Finally
-! the binary is filled with entries of four int(4)'s each. Each 'row' in the
-! binary contains four int(4)'s: 3 for x,y,z position and 1 for cell or material
-! id.  For 1 million voxels this produces a file of approximately 15MB.
-!===============================================================================
-
-  subroutine create_voxel(pl)
-    type(ObjectPlot), intent(in) :: pl
-
-    integer(C_INT) :: x, y, z      ! voxel location indices
-    integer :: rgb(3)       ! colors (red, green, blue) from 0-255
-    integer :: id           ! id of cell or material
-    integer(C_INT), target :: data(pl%pixels(3),pl%pixels(2))
-    integer(HID_T) :: file_id
-    integer(HID_T) :: dspace
-    integer(HID_T) :: memspace
-    integer(HID_T) :: dset
-    integer(HSIZE_T) :: dims(3)
-    real(8) :: vox(3)       ! x, y, and z voxel widths
-    real(8) :: ll(3)        ! lower left starting point for each sweep direction
-    type(Particle)    :: p
-    type(ProgressBar) :: progress
-
-    interface
-      subroutine voxel_init(file_id, dims, dspace, dset, memspace) bind(C)
-        import HID_T, HSIZE_T
-        integer(HID_T), value :: file_id
-        integer(HSIZE_T), intent(in) :: dims(*)
-        integer(HID_T), intent(out) :: dspace
-        integer(HID_T), intent(out) :: dset
-        integer(HID_T), intent(out) :: memspace
-      end subroutine voxel_init
-
-      subroutine voxel_write_slice(x, dspace, dset, memspace, buf) bind(C)
-        import C_INT, HID_T, C_PTR
-        integer(C_INT), value :: x
-        integer(HID_T), value :: dspace
-        integer(HID_T), value :: dset
-        integer(HID_T), value :: memspace
-        type(C_PTR), value :: buf
-      end subroutine voxel_write_slice
-
-      subroutine voxel_finalize(dspace, dset, memspace) bind(C)
-        import HID_T
-        integer(HID_T), value :: dspace
-        integer(HID_T), value :: dset
-        integer(HID_T), value :: memspace
-      end subroutine voxel_finalize
-    end interface
-
-    ! compute voxel widths in each direction
-    vox = pl % width/dble(pl % pixels)
-
-    ! initial particle position
-    ll = pl % origin - pl % width / TWO
-
-    ! allocate and initialize particle
-    call particle_initialize(p)
-    p % coord(1) % xyz = ll
-    p % coord(1) % uvw = [ HALF, HALF, HALF ]
-    p % coord(1) % universe = root_universe
-
-    ! Open binary plot file for writing
-    file_id = file_open(pl%path_plot, 'w')
-
-    ! write header info
-    call write_attribute(file_id, "filetype", 'voxel')
-    call write_attribute(file_id, "version", VERSION_VOXEL)
-    call write_attribute(file_id, "openmc_version", VERSION)
-#ifdef GIT_SHA1
-    call write_attribute(file_id, "git_sha1", GIT_SHA1)
-#endif
-
-    ! Write current date and time
-    call write_attribute(file_id, "date_and_time", time_stamp())
-
-    call write_attribute(file_id, "num_voxels", pl%pixels)
-    call write_attribute(file_id, "voxel_width", vox)
-    call write_attribute(file_id, "lower_left", ll)
-
-    ! Create dataset for voxel data -- note that the dimensions are reversed
-    ! since we want the order in the file to be z, y, x
-    dims(:) = pl % pixels
-    call voxel_init(file_id, dims, dspace, dset, memspace)
-
-    ! move to center of voxels
-    ll = ll + vox / TWO
-
-    do x = 1, pl % pixels(1)
-      call progress % set_value(dble(x)/dble(pl % pixels(1))*100)
-      do y = 1, pl % pixels(2)
-        do z = 1, pl % pixels(3)
-          ! get voxel color
-          call position_rgb(p, pl, rgb, id)
-
-          ! write to plot file
-          data(z,y) = id
-
-          ! advance particle in z direction
-          p % coord(1) % xyz(3) = p % coord(1) % xyz(3) + vox(3)
-        end do
-
-        ! advance particle in y direction
-        p % coord(1) % xyz(2) = p % coord(1) % xyz(2) + vox(2)
-        p % coord(1) % xyz(3) = ll(3)
-      end do
-
-      ! advance particle in y direction
-      p % coord(1) % xyz(1) = p % coord(1) % xyz(1) + vox(1)
-      p % coord(1) % xyz(2) = ll(2)
-      p % coord(1) % xyz(3) = ll(3)
-
-      ! Write to HDF5 dataset
-      call voxel_write_slice(x, dspace, dset, memspace, c_loc(data))
-    end do
-
-    call voxel_finalize(dspace, dset, memspace)
-    call file_close(file_id)
-
-  end subroutine create_voxel
-
-end module plot
diff --git a/src/plot.cpp b/src/plot.cpp
index 79f04c47f..7b0cddb63 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -1,10 +1,931 @@
+#include <fstream>
+#include <sstream>
+
 #include "openmc/plot.h"
+#include "openmc/constants.h"
+#include "openmc/settings.h"
+#include "openmc/error.h"
+#include "openmc/particle.h"
+#include "openmc/geometry.h"
+#include "openmc/cell.h"
+#include "openmc/material.h"
+#include "openmc/string_functions.h"
+#include "openmc/mesh.h"
+#include "openmc/output.h"
+#include "openmc/hdf5_interface.h"
+#include "openmc/random_lcg.h"
+#include "openmc/output.h"
+#include "openmc/progress_bar.h"
 
 namespace openmc {
 
+//==============================================================================
+// Global variables
+//==============================================================================
+
+int PLOT_LEVEL_LOWEST = -1;
+
+std::unordered_map<int, int> plot_map;
+
+int n_plots;
+
+std::vector<Plot> plots;
+
+const RGBColor WHITE = {255, 255, 255};
+
+//==============================================================================
+// RUN_PLOT controls the logic for making one or many plots
+//==============================================================================
+
+extern "C"
+int openmc_plot_geometry()
+{
+  int err;
+
+  for (auto pl : plots) {
+    std::stringstream ss;
+    ss << "Processing plot " << pl.id_ << ": "
+       << pl.path_plot_ << "...";
+    write_message(ss.str(), 5);
+
+    if (PlotType::slice == pl.type_) {
+      // create 2D image
+      create_ppm(pl);
+    } else if (PlotType::voxel == pl.type_) {
+      // create voxel file for 3D viewing
+      create_voxel(pl);
+    }
+  }
+  return 0;
+}
+
+
 void
-voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, hid_t* dset,
-           hid_t* memspace)
+read_plots(pugi::xml_node* plots_node)
+{
+  n_plots = 0;
+  for (auto node : plots_node->children("plot")) {
+    Plot pl(node);
+    plots.push_back(pl);
+    plot_map[pl.id_] = n_plots++;
+  }
+}
+
+//==============================================================================
+// CREATE_PPM creates an image based on user input from a plots.xml <plot>
+// specification in the portable pixmap format (PPM)
+//==============================================================================
+
+void create_ppm(Plot pl)
+{
+
+  size_t width = pl.pixels_[0];
+  size_t height = pl.pixels_[1];
+
+  double in_pixel = (pl.width_[0])/static_cast<double>(width);
+  double out_pixel = (pl.width_[1])/static_cast<double>(height);
+
+  ImageData data;
+  data.resize({width, height});
+
+  int in_i, out_i;
+  double xyz[3];
+  switch(pl.basis_) {
+  case PlotBasis::xy :
+    in_i = 0;
+    out_i = 1;
+    xyz[0] = pl.origin_[0] - pl.width_[0] / 2.;
+    xyz[1] = pl.origin_[1] + pl.width_[1] / 2.;
+    xyz[2] = pl.origin_[2];
+    break;
+  case PlotBasis::xz :
+    in_i = 0;
+    out_i = 2;
+    xyz[0] = pl.origin_[0] - pl.width_[0] / 2.;
+    xyz[1] = pl.origin_[1];
+    xyz[2] = pl.origin_[2] + pl.width_[1] / 2.;
+    break;
+  case PlotBasis::yz :
+    in_i = 1;
+    out_i = 2;
+    xyz[0] = pl.origin_[0];
+    xyz[1] = pl.origin_[1] - pl.width_[0] / 2.;
+    xyz[2] = pl.origin_[2] + pl.width_[1] / 2.;
+    break;
+  }
+
+  double dir[3] = {0.5, 0.5, 0.5};
+
+#pragma omp parallel
+{
+  Particle p;
+  p.initialize();
+  std::copy(xyz, xyz+3, p.coord[0].xyz);
+  std::copy(dir, dir+3, p.coord[0].uvw);
+  p.coord[0].universe = openmc_root_universe;
+
+#pragma omp for
+  for (int y = 0; y < height; y++) {
+    p.coord[0].xyz[out_i] = xyz[out_i] - out_pixel * y;
+    for (int x = 0; x < width; x++) {
+      // local variables
+      RGBColor rgb;
+      int id;
+      p.coord[0].xyz[in_i] = xyz[in_i] + in_pixel * x;
+      position_rgb(p, pl, rgb, id);
+      data(x,y) = rgb;
+    }
+  }
+}
+  // draw mesh lines if present
+  if (pl.index_meshlines_mesh_ >= 0) {draw_mesh_lines(pl, data);}
+
+  // write ppm data to file
+  output_ppm(pl, data);
+}
+
+void
+Plot::set_id(pugi::xml_node plot_node)
+{
+  // Copy data into plots
+  if (check_for_node(plot_node, "id")) {
+    id_ = std::stoi(get_node_value(plot_node, "id"));
+  } else {
+    fatal_error("Must specify plot id in plots XML file.");
+  }
+
+  // Check to make sure 'id' hasn't been used
+  if (plot_map.find(id_) != plot_map.end()) {
+    std::stringstream err_msg;
+    err_msg << "Two or more plots use the same unique ID: " << id_;
+    fatal_error(err_msg.str());
+  }
+}
+
+void
+Plot::set_type(pugi::xml_node plot_node)
+{
+  // Copy plot type
+  // Default is slice
+  type_ = PlotType::slice;
+  // check type specified on plot node
+  if (check_for_node(plot_node, "type")) {
+    std::string type_str = get_node_value(plot_node, "type", true);
+    // set type using node value
+    if (type_str == "slice") {
+      type_ = PlotType::slice;
+    }
+    else if (type_str == "voxel") {
+      type_ = PlotType::voxel;
+    } else {
+    // if we're here, something is wrong
+    std::stringstream err_msg;
+    err_msg << "Unsupported plot type '" << type_str
+            << "' in plot " << id_;
+    fatal_error(err_msg.str());
+    }
+  }
+}
+
+void
+Plot::set_output_path(pugi::xml_node plot_node)
+{
+  // Set output file path
+  std::stringstream filename;
+
+  if (check_for_node(plot_node, "filename")) {
+    filename << get_node_value(plot_node, "filename");
+  } else {
+    filename << "plot_" << id_;
+  }
+  // add appropriate file extension to name
+  switch(type_) {
+  case PlotType::slice:
+    filename << ".ppm";
+    break;
+  case PlotType::voxel:
+    filename << ".h5";
+    break;
+  }
+
+  path_plot_ = filename.str();
+
+  // Copy plot pixel size
+  std::vector<int> pxls = get_node_array<int>(plot_node, "pixels");
+  if (PlotType::slice == type_) {
+    if (pxls.size() == 2) {
+      pixels_[0] = pxls[0];
+      pixels_[1] = pxls[1];
+    } else {
+      std::stringstream err_msg;
+      err_msg << "<pixels> must be length 2 in slice plot "
+              << id_;
+      fatal_error(err_msg.str());
+    }
+  } else if (PlotType::voxel == type_) {
+    if (pxls.size() == 3) {
+      pixels_[0] = pxls[0];
+      pixels_[1] = pxls[1];
+      pixels_[2] = pxls[2];
+    } else {
+      std::stringstream err_msg;
+      err_msg << "<pixels> must be length 3 in voxel plot "
+              << id_;
+      fatal_error(err_msg.str());
+    }
+  }
+}
+
+void
+Plot::set_bg_color(pugi::xml_node plot_node)
+{
+  // Copy plot background color
+  if (check_for_node(plot_node, "background")) {
+    std::vector<int> bg_rgb = get_node_array<int>(plot_node, "background");
+    if (PlotType::voxel == type_) {
+      if (openmc_master) {
+        std::stringstream err_msg;
+        err_msg << "Background color ignored in voxel plot "
+                << id_;
+        warning(err_msg.str());
+      }
+    }
+    if (bg_rgb.size() == 3) {
+      not_found_ = bg_rgb;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Bad background RGB in plot "
+              << id_;
+      fatal_error(err_msg);
+    }
+  } else {
+    // default to a white background
+    not_found_ = WHITE;
+  }
+}
+
+void
+Plot::set_basis(pugi::xml_node plot_node)
+{
+  // Copy plot basis
+  if (PlotType::slice == type_) {
+    std::string pl_basis = "xy";
+    if (check_for_node(plot_node, "basis")) {
+      pl_basis = get_node_value(plot_node, "basis", true);
+    }
+    if ("xy" == pl_basis) {
+      basis_ = PlotBasis::xy;
+    } else if ("xz" == pl_basis) {
+      basis_ = PlotBasis::xz;
+    } else if ("yz" == pl_basis) {
+      basis_ = PlotBasis::yz;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Unsupported plot basis '" << pl_basis
+              << "' in plot " << id_;
+      fatal_error(err_msg);
+    }
+  }
+}
+
+void
+Plot::set_origin(pugi::xml_node plot_node)
+{
+  // Copy plotting origin
+  std::vector<double> pl_origin = get_node_array<double>(plot_node, "origin");
+  if (pl_origin.size() == 3) {
+    origin_[0] = pl_origin[0];
+    origin_[1] = pl_origin[1];
+    origin_[2] = pl_origin[2];
+  } else {
+    std::stringstream err_msg;
+    err_msg << "Origin must be length 3 in plot "
+            << id_;
+    fatal_error(err_msg);
+  }
+}
+
+void
+Plot::set_width(pugi::xml_node plot_node)
+{
+  // Copy plotting width
+  std::vector<double> pl_width = get_node_array<double>(plot_node, "width");
+  if (PlotType::slice == type_) {
+    if (pl_width.size() == 2) {
+      width_[0] = pl_width[0];
+      width_[1] = pl_width[1];
+    } else {
+      std::stringstream err_msg;
+      err_msg << "<width> must be length 2 in slice plot "
+              << id_;
+      fatal_error(err_msg);
+    }
+  } else if (PlotType::voxel == type_) {
+    if (pl_width.size() == 3) {
+      pl_width = get_node_array<double>(plot_node, "width");
+      width_[0] = pl_width[0];
+      width_[1] = pl_width[1];
+      width_[2] = pl_width[2];
+    } else {
+      std::stringstream err_msg;
+      err_msg << "<width> must be length 3 in voxel plot "
+              << id_;
+      fatal_error(err_msg);
+    }
+  }
+}
+
+void
+Plot::set_universe(pugi::xml_node plot_node)
+{
+  // Copy plot universe level
+  if (check_for_node(plot_node, "level")) {
+    level_ = std::stoi(get_node_value(plot_node, "level"));
+    if (level_ < 0) {
+      std::stringstream err_msg;
+      err_msg << "Bad universe level in plot " << id_;
+      fatal_error(err_msg);
+    }
+  } else {
+    level_ = PLOT_LEVEL_LOWEST;
+  }
+}
+
+void
+Plot::set_default_colors(pugi::xml_node plot_node)
+{
+  // Copy plot color type and initialize all colors randomly
+  std::string pl_color_by = "cell";
+  if (check_for_node(plot_node, "color_by")) {
+    pl_color_by = get_node_value(plot_node, "color_by", true);
+  }
+  if ("cell" == pl_color_by) {
+    color_by_ = PlotColorBy::cells;
+    colors_.resize(n_cells);
+    for (int i = 0; i < n_cells; i++) {
+      colors_[i] = random_color();
+    }
+
+  } else if("material" == pl_color_by) {
+    color_by_ = PlotColorBy::mats;
+    colors_.resize(n_materials);
+    for (int i = 0; i < materials.size(); i++) {
+      colors_[i] = random_color();
+    }
+  } else {
+    std::stringstream err_msg;
+    err_msg << "Unsupported plot color type '" << pl_color_by
+            << "' in plot " << id_;
+    fatal_error(err_msg);
+  }
+}
+
+void
+Plot::set_user_colors(pugi::xml_node plot_node)
+{
+  if (!plot_node.select_nodes("color").empty() && PlotType::voxel == type_) {
+    if (openmc_master) {
+      std::stringstream err_msg;
+      err_msg << "Color specifications ignored in voxel plot "
+              << id_;
+      warning(err_msg);
+    }
+  }
+
+  for (auto cn : plot_node.children("color")) {
+    // Make sure 3 values are specified for RGB
+    std::vector<int> user_rgb = get_node_array<int>(cn, "rgb");
+    if (user_rgb.size() != 3) {
+      std::stringstream err_msg;
+      err_msg << "Bad RGB in plot " << id_;
+      fatal_error(err_msg);
+    }
+    // Ensure that there is an id for this color specification
+    int col_id;
+    if (check_for_node(cn, "id")) {
+      col_id = std::stoi(get_node_value(cn, "id"));
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Must specify id for color specification in plot "
+              << id_;
+      fatal_error(err_msg);
+    }
+    // Add RGB
+    if (PlotColorBy::cells == color_by_) {
+      if (cell_map.find(col_id) != cell_map.end()) {
+        col_id = cell_map[col_id];
+        colors_[col_id] = user_rgb;
+      } else {
+        std::stringstream err_msg;
+        err_msg << "Could not find cell " << col_id
+                << " specified in plot " << id_;
+        fatal_error(err_msg);
+      }
+    } else if (PlotColorBy::mats == color_by_) {
+      if (material_map.find(col_id) != material_map.end()) {
+        col_id = material_map[col_id];
+        colors_[col_id] = user_rgb;
+      } else {
+        std::stringstream err_msg;
+        err_msg << "Could not find material " << col_id
+                << " specified in plot " << id_;
+        fatal_error(err_msg);
+      }
+    }
+  } // color node loop
+}
+
+void
+Plot::set_meshlines(pugi::xml_node plot_node)
+{
+  // Deal with meshlines
+  pugi::xpath_node_set mesh_line_nodes = plot_node.select_nodes("meshlines");
+
+  if (!mesh_line_nodes.empty()) {
+    if (PlotType::voxel == type_) {
+      std::stringstream msg;
+      msg << "Meshlines ignored in voxel plot " << id_;
+      warning(msg);
+    }
+
+    if (mesh_line_nodes.size() == 1) {
+      // Get first meshline node
+      pugi::xml_node meshlines_node = mesh_line_nodes[0].node();
+
+      // Check mesh type
+      std::string meshtype;
+      if (check_for_node(meshlines_node, "meshtype")) {
+        meshtype = get_node_value(meshlines_node, "meshtype");
+      } else {
+        std::stringstream err_msg;
+        err_msg << "Must specify a meshtype for meshlines specification in plot " << id_;
+        fatal_error(err_msg);
+      }
+
+      // Ensure that there is a linewidth for this meshlines specification
+      std::string meshline_width;
+      if (check_for_node(meshlines_node, "linewidth")) {
+        meshline_width = get_node_value(meshlines_node, "linewidth");
+        meshlines_width_ = std::stoi(meshline_width);
+      } else {
+        std::stringstream err_msg;
+        err_msg << "Must specify a linewidth for meshlines specification in plot " << id_;
+        fatal_error(err_msg);
+      }
+
+      // Check for color
+      if (check_for_node(meshlines_node, "color")) {
+        // Check and make sure 3 values are specified for RGB
+        std::vector<int> ml_rgb = get_node_array<int>(meshlines_node, "color");
+        if (ml_rgb.size() != 3) {
+          std::stringstream err_msg;
+          err_msg << "Bad RGB for meshlines color in plot " << id_;
+          fatal_error(err_msg);
+        }
+        meshlines_color_ = ml_rgb;
+      }
+
+      // Set mesh based on type
+      if ("ufs" == meshtype) {
+        if (settings::index_ufs_mesh < 0) {
+          std::stringstream err_msg;
+          err_msg << "No UFS mesh for meshlines on plot " << id_;
+          fatal_error(err_msg);
+        } else {
+          index_meshlines_mesh_ = settings::index_ufs_mesh;
+        }
+      } else if ("cmfd" == meshtype) {
+        if (!settings::cmfd_run) {
+          std::stringstream err_msg;
+          err_msg << "Need CMFD run to plot CMFD mesh for meshlines on plot " << id_;
+          fatal_error(err_msg);
+        } else {
+          index_meshlines_mesh_ = settings::index_cmfd_mesh;
+        }
+      } else if ("entropy" == meshtype) {
+        if (settings::index_entropy_mesh < 0) {
+          std::stringstream err_msg;
+          err_msg <<"No entropy mesh for meshlines on plot " << id_;
+          fatal_error(err_msg);
+        } else {
+          index_meshlines_mesh_ = settings::index_entropy_mesh;
+        }
+      } else if ("tally" == meshtype) {
+        // Ensure that there is a mesh id if the type is tally
+        int tally_mesh_id;
+        if (check_for_node(meshlines_node, "id")) {
+          tally_mesh_id = std::stoi(get_node_value(meshlines_node, "id"));
+        } else {
+          std::stringstream err_msg;
+          err_msg << "Must specify a mesh id for meshlines tally "
+                  << "mesh specification in plot " << id_;
+          fatal_error(err_msg);
+        }
+        // find the tally index
+        int idx;
+        int err = openmc_get_mesh_index(tally_mesh_id, &idx);
+        if (err != 0) {
+          std::stringstream err_msg;
+          err_msg << "Could not find mesh " << tally_mesh_id
+                  << " specified in meshlines for plot " << id_;
+          fatal_error(err_msg);
+        }
+        index_meshlines_mesh_ = idx;
+      } else {
+        std::stringstream err_msg;
+        err_msg << "Invalid type for meshlines on plot " << id_ ;
+        fatal_error(err_msg);
+      }
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Mutliple meshlines specified in plot " << id_;
+      fatal_error(err_msg);
+    }
+  }
+}
+
+void
+Plot::set_mask(pugi::xml_node plot_node)
+{
+  // Deal with masks
+  pugi::xpath_node_set mask_nodes = plot_node.select_nodes("mask");
+
+  if (!mask_nodes.empty()) {
+    if (PlotType::voxel == type_) {
+      if (openmc_master) {
+        std::stringstream wrn_msg;
+        wrn_msg << "Mask ignored in voxel plot " << id_;
+        warning(wrn_msg);
+      }
+    }
+
+    if (mask_nodes.size() == 1) {
+      // Get pointer to mask
+      pugi::xml_node mask_node = mask_nodes[0].node();
+
+      // Determine how many components there are and allocate
+      std::vector<int> iarray = get_node_array<int>(mask_node, "components");
+      if (iarray.size() == 0) {
+        std::stringstream err_msg;
+        err_msg << "Missing <components> in mask of plot " << id_;
+        fatal_error(err_msg);
+      }
+
+      // First we need to change the user-specified identifiers to indices
+      // in the cell and material arrays
+      for (auto& col_id : iarray) {
+        if (PlotColorBy::cells == color_by_) {
+          if (cell_map.find(col_id) != cell_map.end()) {
+            col_id  = cell_map[col_id];
+          }
+          else {
+            std::stringstream err_msg;
+            err_msg << "Could not find cell " << col_id
+                    << " specified in the mask in plot " << id_;
+            fatal_error(err_msg);
+          }
+        } else if (PlotColorBy::mats == color_by_) {
+          if (material_map.find(col_id) != material_map.end()) {
+            col_id = material_map[col_id];
+          }
+          else {
+            std::stringstream err_msg;
+            err_msg << "Could not find material " << col_id
+                    << " specified in the mask in plot " << id_;
+            fatal_error(err_msg);
+          }
+        }
+      }
+
+      // Alter colors based on mask information
+      for (int j = 0; j < colors_.size(); j++) {
+        if (std::find(iarray.begin(), iarray.end(), j) == iarray.end()) {
+          if (check_for_node(mask_node, "background")) {
+            std::vector<int> bg_rgb = get_node_array<int>(mask_node, "background");
+            colors_[j] = bg_rgb;
+          } else {
+            colors_[j] = WHITE;
+          }
+        }
+      }
+
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Mutliple masks specified in plot " << id_;
+      fatal_error(err_msg);
+    }
+  }
+}
+
+Plot::Plot(pugi::xml_node plot_node):
+index_meshlines_mesh_(-1)
+{
+  set_id(plot_node);
+  set_type(plot_node);
+  set_output_path(plot_node);
+  set_bg_color(plot_node);
+  set_basis(plot_node);
+  set_origin(plot_node);
+  set_width(plot_node);
+  set_universe(plot_node);
+  set_default_colors(plot_node);
+  set_user_colors(plot_node);
+  set_meshlines(plot_node);
+  set_mask(plot_node);
+} // End Plot constructor
+
+//==============================================================================
+// POSITION_RGB computes the red/green/blue values for a given plot with the
+// current particle's position
+//==============================================================================
+
+
+void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id)
+{
+  p.n_coord = 1;
+
+  bool found_cell = find_cell(&p, 0);
+
+  int j = p.n_coord - 1;
+
+  if (settings::check_overlaps) {check_cell_overlap(&p);}
+
+  // Set coordinate level if specified
+  if (pl.level_ >= 0) {j = pl.level_ + 1;}
+
+  if (!found_cell) {
+    // If no cell, revert to default color
+    rgb = pl.not_found_;
+    id = -1;
+  } else {
+    if (PlotColorBy::mats == pl.color_by_) {
+      // Assign color based on material
+      Cell* c = cells[p.coord[j].cell];
+      if (c->type_ == FILL_UNIVERSE) {
+        // If we stopped on a middle universe level, treat as if not found
+        rgb = pl.not_found_;
+        id = -1;
+      } else if (p.material == MATERIAL_VOID) {
+        // By default, color void cells white
+        rgb = WHITE;
+        id = -1;
+      } else {
+        rgb = pl.colors_[p.material - 1];
+        id = materials[p.material - 1]->id_;
+      }
+    } else if (PlotColorBy::cells == pl.color_by_) {
+      // Assign color based on cell
+      rgb = pl.colors_[p.coord[j].cell];
+      id = cells[p.coord[j].cell]->id_;
+    }
+  } // endif found_cell
+}
+
+//==============================================================================
+// OUTPUT_PPM writes out a previously generated image to a PPM file
+//==============================================================================
+
+void output_ppm(Plot pl, const ImageData& data)
+{
+  // Open PPM file for writing
+  std::string fname = pl.path_plot_;
+  fname = strtrim(fname);
+  std::ofstream of;
+
+  of.open(fname);
+
+  // Write header
+  of << "P6" << "\n";
+  of << pl.pixels_[0] << " " << pl.pixels_[1] << "\n";
+  of << "255" << "\n";
+  of.close();
+
+  of.open(fname, std::ios::binary | std::ios::app);
+  // Write color for each pixel
+  for (int y = 0; y < pl.pixels_[1]; y++) {
+    for (int x = 0; x < pl.pixels_[0]; x++) {
+      RGBColor rgb = data(x,y);
+      of << rgb.red << rgb.green << rgb.blue;
+    }
+  }
+
+  // Close file
+  // THIS IS HERE TO MATCH FORTRAN VERSION, NOT TECHNICALLY NECESSARY
+  of << "\n";
+  of.close();
+}
+
+//==============================================================================
+// DRAW_MESH_LINES draws mesh line boundaries on an image
+//==============================================================================
+
+void draw_mesh_lines(Plot pl, ImageData& data)
+{
+  RGBColor rgb;
+  rgb = pl.meshlines_color_;
+
+  int outer, inner;
+  switch(pl.basis_) {
+  case PlotBasis::xy :
+    outer = 0;
+    inner = 1;
+    break;
+  case PlotBasis::xz :
+    outer = 0;
+    inner = 2;
+    break;
+  case PlotBasis::yz :
+    outer = 1;
+    inner = 2;
+    break;
+  }
+
+  double xyz_ll_plot[3], xyz_ur_plot[3];
+  xyz_ll_plot[0] = pl.origin_[0];
+  xyz_ll_plot[1] = pl.origin_[1];
+  xyz_ll_plot[2] = pl.origin_[2];
+
+  xyz_ur_plot[0] = pl.origin_[0];
+  xyz_ur_plot[1] = pl.origin_[1];
+  xyz_ur_plot[2] = pl.origin_[2];
+
+  xyz_ll_plot[outer] = pl.origin_[outer] - pl.width_[0] / 2.;
+  xyz_ll_plot[inner] = pl.origin_[inner] - pl.width_[1] / 2.;
+  xyz_ur_plot[outer] = pl.origin_[outer] + pl.width_[0] / 2.;
+  xyz_ur_plot[inner] = pl.origin_[inner] + pl.width_[1] / 2.;
+
+  int width[3];
+  width[0] = xyz_ur_plot[0] - xyz_ll_plot[0];
+  width[1] = xyz_ur_plot[1] - xyz_ll_plot[1];
+  width[2] = xyz_ur_plot[2] - xyz_ll_plot[2];
+
+  auto& m = meshes[pl.index_meshlines_mesh_];
+
+  int ijk_ll[3], ijk_ur[3];
+  bool in_mesh;
+  m->get_indices(Position(xyz_ll_plot), &(ijk_ll[0]), &in_mesh);
+  m->get_indices(Position(xyz_ur_plot), &(ijk_ur[0]), &in_mesh);
+
+  // Fortran/C++ index correction
+  ijk_ur[0]++; ijk_ur[1]++; ijk_ur[2]++;
+
+  double xyz_ll[3], xyz_ur[3];
+  // sweep through all meshbins on this plane and draw borders
+  for (int i = ijk_ll[outer]; i <= ijk_ur[outer]; i++) {
+    for (int j = ijk_ll[inner]; j <= ijk_ur[inner]; j++) {
+      // check if we're in the mesh for this ijk
+      if (i > 0 && i <= m->shape_[outer] && j >0 && j <= m->shape_[inner] ) {
+        int outrange[3], inrange[3];
+        // get xyz's of lower left and upper right of this mesh cell
+        xyz_ll[outer] = m->lower_left_[outer] + m->width_[outer] * (i - 1);
+        xyz_ll[inner] = m->lower_left_[inner] + m->width_[inner] * (j - 1);
+        xyz_ur[outer] = m->lower_left_[outer] + m->width_[outer] * i;
+        xyz_ur[inner] = m->lower_left_[inner] + m->width_[inner] * j;
+
+        // map the xyz ranges to pixel ranges
+        double frac = (xyz_ll[outer] - xyz_ll_plot[outer]) / width[outer];
+        outrange[0] = int(frac * double(pl.pixels_[0]));
+        frac = (xyz_ur[outer] - xyz_ll_plot[outer]) / width[outer];
+        outrange[1] = int(frac * double(pl.pixels_[0]));
+
+        frac = (xyz_ur[inner] - xyz_ll_plot[inner]) / width[inner];
+        inrange[0] = int((1. - frac) * (double)pl.pixels_[1]);
+        frac = (xyz_ll[inner] - xyz_ll_plot[inner]) / width[inner];
+        inrange[1] = int((1. - frac) * (double)pl.pixels_[1]);
+
+        // draw lines
+        for (int out_ = outrange[0]; out_ <= outrange[1]; out_++) {
+          for (int plus = 0; plus <= pl.meshlines_width_; plus++) {
+            data(out_, inrange[0] + plus) = rgb;
+            data(out_, inrange[1] + plus) = rgb;
+            data(out_, inrange[0] - plus) = rgb;
+            data(out_, inrange[1] - plus) = rgb;
+          }
+        }
+
+        for (int in_ = inrange[0]; in_ <= inrange[1]; in_++) {
+          for (int plus = 0; plus <= pl.meshlines_width_; plus++) {
+            data(outrange[0] + plus, in_) = rgb;
+            data(outrange[1] + plus, in_) = rgb;
+            data(outrange[0] - plus, in_) = rgb;
+            data(outrange[1] - plus, in_) = rgb;
+          }
+        }
+
+      } // end if(in mesh)
+    }
+  } // end outer loops
+} // end draw_mesh_lines
+
+//==============================================================================
+// CREATE_VOXEL outputs a binary file that can be input into silomesh for 3D
+// geometry visualization.  It works the same way as create_ppm by dragging a
+// particle across the geometry for the specified number of voxels. The first 3
+// int(4)'s in the binary are the number of x, y, and z voxels.  The next 3
+// real(8)'s are the widths of the voxels in the x, y, and z directions. The
+// next 3 real(8)'s are the x, y, and z coordinates of the lower left
+// point. Finally the binary is filled with entries of four int(4)'s each. Each
+// 'row' in the binary contains four int(4)'s: 3 for x,y,z position and 1 for
+// cell or material id.  For 1 million voxels this produces a file of
+// approximately 15MB.
+// =============================================================================
+
+void create_voxel(Plot pl)
+{
+
+  // compute voxel widths in each direction
+  std::array<double, 3> vox;
+  vox[0] = pl.width_[0]/(double)pl.pixels_[0];
+  vox[1] = pl.width_[1]/(double)pl.pixels_[1];
+  vox[2] = pl.width_[2]/(double)pl.pixels_[2];
+
+  // initial particle position
+  std::array<double, 3> ll;
+  ll[0] = pl.origin_[0] - pl.width_[0] / 2.;
+  ll[1] = pl.origin_[1] - pl.width_[1] / 2.;
+  ll[2] = pl.origin_[2] - pl.width_[2] / 2.;
+
+  // allocate and initialize particle
+  double dir[3] = {0.5, 0.5, 0.5};
+  Particle p;
+  p.initialize();
+  std::copy(ll.begin(), ll.begin()+ll.size(), p.coord[0].xyz);
+  std::copy(dir, dir+3, p.coord[0].uvw);
+  p.coord[0].universe = openmc_root_universe;
+
+  // Open binary plot file for writing
+  std::ofstream of;
+  std::string fname = std::string(pl.path_plot_);
+  fname = strtrim(fname);
+  hid_t file_id = file_open(fname, 'w');
+
+  // write header info
+  write_attribute(file_id, "filetype", "voxel");
+  write_attribute(file_id, "version", VERSION_VOXEL);
+  write_attribute(file_id, "openmc_version", VERSION);
+
+#ifdef GIT_SHA1
+  write_attribute(file_id, "git_sha1", GIT_SHA1);
+#endif
+
+  // Write current date and time
+  write_attribute(file_id, "date_and_time", time_stamp().c_str());
+  hsize_t three = 3;
+  write_attribute(file_id, "num_voxels", pl.pixels_);
+  write_attribute(file_id, "voxel_width", vox);
+  write_attribute(file_id, "lower_left", ll);
+
+  // Create dataset for voxel data -- note that the dimensions are reversed
+  // since we want the order in the file to be z, y, x
+  hsize_t dims[3];
+  dims[0] = pl.pixels_[0];
+  dims[1] = pl.pixels_[1];
+  dims[2] = pl.pixels_[2];
+  hid_t dspace, dset, memspace;
+  voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace);
+
+  // move to center of voxels
+  ll[0] = ll[0] + vox[0] / 2.;
+  ll[1] = ll[1] + vox[1] / 2.;
+  ll[2] = ll[2] + vox[2] / 2.;
+
+  int data[pl.pixels_[1]][pl.pixels_[2]];
+
+  ProgressBar pb;
+
+  RGBColor rgb;
+  int id;
+  for (int x = 0; x < pl.pixels_[0]; x++) {
+    pb.set_value(100.*(double)x/(double)(pl.pixels_[0]-1));
+    for (int y = 0; y < pl.pixels_[1]; y++) {
+      for (int z = 0; z < pl.pixels_[2]; z++) {
+        // get voxel color
+        position_rgb(p, pl, rgb, id);
+        // write to plot data
+        data[y][z] = id;
+        // advance particle in z direction
+        p.coord[0].xyz[2] = p.coord[0].xyz[2] + vox[2];
+      }
+      // advance particle in y direction
+      p.coord[0].xyz[1] = p.coord[0].xyz[1] + vox[1];
+      p.coord[0].xyz[2] = ll[2];
+    }
+    // advance particle in x direction
+    p.coord[0].xyz[0] = p.coord[0].xyz[0] + vox[0];
+    p.coord[0].xyz[1] = ll[1];
+    p.coord[0].xyz[2] = ll[2];
+    // Write to HDF5 dataset
+    voxel_write_slice(x, dspace, dset, memspace, &(data[0]));
+  }
+
+  voxel_finalize(dspace, dset, memspace);
+  file_close(file_id);
+
+}
+
+void
+voxel_init(hid_t file_id, const hsize_t* dims,
+           hid_t* dspace, hid_t* dset, hid_t* memspace)
 {
   // Create dataspace/dataset for voxel data
   *dspace = H5Screate_simple(3, dims, nullptr);
@@ -25,7 +946,7 @@ voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, hid_t* dset,
 void
 voxel_write_slice(int x, hid_t dspace, hid_t dset, hid_t memspace, void* buf)
 {
-  hssize_t offset[3] {x - 1, 0, 0};
+  hssize_t offset[3] {x, 0, 0};
   H5Soffset_simple(dspace, offset);
   H5Dwrite(dset, H5T_NATIVE_INT, memspace, dspace, H5P_DEFAULT, buf);
 }
@@ -39,4 +960,8 @@ voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace)
   H5Sclose(memspace);
 }
 
+RGBColor random_color() {
+  return {int(prn()*255), int(prn()*255), int(prn()*255)};
+}
+
 } // namespace openmc
diff --git a/src/plot_header.F90 b/src/plot_header.F90
deleted file mode 100644
index 761481b33..000000000
--- a/src/plot_header.F90
+++ /dev/null
@@ -1,74 +0,0 @@
-module plot_header
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants
-  use dict_header, only: DictIntInt
-
-  implicit none
-
-!===============================================================================
-! ObjectColor holds color information for plotted objects
-!===============================================================================
-
-  type ObjectColor
-    integer :: rgb(3)
-  end type ObjectColor
-
-!===============================================================================
-! PLOTSLICE holds plot information
-!===============================================================================
-
-  type ObjectPlot
-    integer :: id                    ! Unique ID
-    character(MAX_LINE_LEN) :: path_plot ! path for plot file
-    integer :: type                  ! Type
-    integer :: color_by              ! quantity to color regions by
-    real(8) :: origin(3)             ! xyz center of plot location
-    real(8) :: width(3)              ! xyz widths of plot
-    integer :: basis                 ! direction of plot slice
-    integer :: pixels(3)             ! pixel width/height of plot slice
-    integer :: meshlines_width       ! pixel width of meshlines
-    integer :: level                 ! universe depth to plot the cells of
-    integer :: index_meshlines_mesh = -1 ! index of  mesh to plot
-    type(ObjectColor) :: meshlines_color ! Color for meshlines
-    type(ObjectColor) :: not_found   ! color for positions where no cell found
-    type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats
-  end type ObjectPlot
-
-  ! Plot type
-  integer, parameter :: PLOT_TYPE_SLICE = 1
-  integer, parameter :: PLOT_TYPE_VOXEL = 2
-
-  ! Plot level
-  integer, parameter :: PLOT_LEVEL_LOWEST = -1
-
-  ! Plot basis plane
-  integer, parameter :: PLOT_BASIS_XY = 1
-  integer, parameter :: PLOT_BASIS_XZ = 2
-  integer, parameter :: PLOT_BASIS_YZ = 3
-
-  ! Indicate whether color refers to unique cell or unique material
-  integer, parameter :: PLOT_COLOR_CELLS = 1
-  integer, parameter :: PLOT_COLOR_MATS = 2
-
-  integer(C_INT32_T), bind(C) :: n_plots     ! # of plots
-
-  type(ObjectPlot), allocatable, target :: plots(:)
-
-  ! Dictionary that maps user IDs to indices in 'plots'
-  type(DictIntInt) :: plot_dict
-
-contains
-
-!===============================================================================
-! FREE_MEMORY_PLOT deallocates global arrays defined in this module
-!===============================================================================
-
-  subroutine free_memory_plot()
-    n_plots = 0
-    if (allocated(plots)) deallocate(plots)
-    call plot_dict % clear()
-  end subroutine free_memory_plot
-
-end module plot_header
diff --git a/src/progress_bar.cpp b/src/progress_bar.cpp
new file mode 100644
index 000000000..7917f9ffd
--- /dev/null
+++ b/src/progress_bar.cpp
@@ -0,0 +1,67 @@
+
+#include "openmc/progress_bar.h"
+
+#include <sstream>
+#include <iostream>
+#include <iomanip>
+
+#ifdef UNIX
+#include <unistd.h>
+#endif
+
+#define BAR_WIDTH 72
+
+bool is_terminal() {
+#ifdef UNIX
+  return isatty(STDOUT_FILENO) != 0;
+#else
+  return false;
+#endif
+}
+
+ProgressBar::ProgressBar() {
+  // initialize bar
+  set_value(0.0);
+}
+
+void
+ProgressBar::set_value(double val) {
+
+  if (!is_terminal()) return;
+
+  // set the bar percentage
+  if (val >= 100.0) {
+    bar.append("100");
+  } else if (val <= 0.0) {
+    bar.append("  0");
+  } else {
+    std::stringstream ss;
+    ss << std::setfill(' ') << std::setw(3) << (int)val;
+    bar.append(ss.str());
+  }
+
+  bar.append("% |");
+  // remaining width of the bar
+  int remaining_width = BAR_WIDTH - bar.size() - 2;
+
+  // set the bar width
+  if (val >= 100.0) {
+    bar.append(remaining_width, '=');
+  } else if (val < 0.0) {
+    bar.append(remaining_width, ' ');
+  } else {
+    int width = (int)((double)remaining_width*val/100);
+    bar.append(width, '=');
+    bar.append(1, '>');
+    bar.append(remaining_width-width-1, ' ');
+  }
+
+  bar.append("|+");
+
+  // write the bar
+  std::cout << '\r' << bar << std::flush;
+  if (val >= 100.0) { std::cout << "\n"; }
+
+  // reset the bar value
+  bar = "";
+}
diff --git a/src/progress_header.F90 b/src/progress_header.F90
deleted file mode 100644
index 5f462e0fa..000000000
--- a/src/progress_header.F90
+++ /dev/null
@@ -1,102 +0,0 @@
-module progress_header
-
-  use, intrinsic :: ISO_FORTRAN_ENV,      only: OUTPUT_UNIT
-
-  implicit none
-
-#ifdef UNIX
-  interface
-    function check_isatty(fd) bind(C, name = 'isatty')
-      use, intrinsic :: ISO_C_BINDING, only: c_int
-      integer(c_int)        :: check_isatty
-      integer(c_int), value :: fd
-    end function check_isatty
-  end interface
-#endif
-
-!===============================================================================
-! PROGRESSBAR
-!===============================================================================
-
-  type ProgressBar
-    private
-    character(len=72) :: bar="???% |                                      " // &
-                             "                           |"
-  contains
-    procedure :: set_value => bar_set_value
-  end type ProgressBar
-
-contains
-
-!===============================================================================
-! IS_TERMINAL checks if output is currently being output to a terminal.  Relies
-! on a POSIX implementation of isatty, and defaults to false if that is not
-! available
-!===============================================================================
-
-    function is_terminal() result(istty)
-      logical :: istty
-
-      istty = .true.
-
-#ifdef UNIX
-      if (check_isatty(1) == 0) istty = .false.
-#else
-      istty = .false.
-#endif
-
-    end function is_terminal
-
-!===============================================================================
-! BAR_SET_VALUE prints the progress bar without advancing.  The value is
-! specified as percent completion, from 0 to 100.  If the value is ever set to
-! 100 or above, the bar is set to 100 and a newline is written.
-!===============================================================================
-
-  subroutine bar_set_value(self, val)
-
-    class(ProgressBar), intent(inout) :: self
-    real(8),            intent(in)    :: val
-
-    integer :: i
-
-    if (.not. is_terminal()) return
-
-    ! set the percentage
-    if (val >= 100.) then
-      write(self % bar(1:3), "(I3)") 100
-    else if (val <= 0.) then
-      write(self % bar(1:3), "(I3)") 0
-    else
-      write(self % bar(1:3), "(I3)") int(val)
-    end if
-
-    ! set the bar width
-    if (val >= 100.) then
-      do i=1,65
-        write(self % bar(i+6:i+6), '(A)') '='
-      end do
-    else
-      do i=1,int(dble(65)*val/100.)
-        write(self % bar(i+6:i+6), '(A)') '='
-      end do
-    end if
-
-    write(OUTPUT_UNIT, '(A1,A1,A72)', ADVANCE='no') '+', char(13), self % bar
-    flush(OUTPUT_UNIT)
-
-    if (val >= 100.) then
-
-      ! make new line
-      write(OUTPUT_UNIT, "(A)") ""
-      flush(OUTPUT_UNIT)
-
-      ! reset the bar in case we want to use this instance again
-      self % bar = "???% |                                      " // &
-                   "                           |"
-
-    end if
-
-  end subroutine bar_set_value
-
-end module progress_header
diff --git a/src/settings.cpp b/src/settings.cpp
index ce0925978..2394597c7 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -72,7 +72,8 @@ std::string path_statepoint;
 
 int32_t index_entropy_mesh {-1};
 int32_t index_ufs_mesh {-1};
-
+int32_t index_cmfd_mesh {-1};
+  
 int32_t n_batches;
 int32_t n_inactive {0};
 int32_t gen_per_batch {1};
diff --git a/src/string_functions.cpp b/src/string_functions.cpp
index f81b09aef..27e14c9ba 100644
--- a/src/string_functions.cpp
+++ b/src/string_functions.cpp
@@ -1,4 +1,5 @@
 #include "openmc/string_functions.h"
+#include <sstream>
 
 namespace openmc {
 
@@ -27,4 +28,13 @@ void to_lower(std::string& str)
   for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]);
 }
 
+int word_count(std::string const& str)
+{
+  std::stringstream stream(str);
+  std::string dum;
+  int count = 0;
+  while (stream >> dum) {count++;}
+  return count;
+}
+
 } // namespace openmc
diff --git a/src/surface.cpp b/src/surface.cpp
index bf1c2e58c..265a2f02c 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -8,7 +8,7 @@
 #include "openmc/error.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/xml_interface.h"
-
+#include "openmc/string_functions.h"
 
 namespace openmc {
 
@@ -35,24 +35,6 @@ std::map<int, int> surface_map;
 // Helper functions for reading the "coeffs" node of an XML surface element
 //==============================================================================
 
-int word_count(const std::string &text)
-{
-  bool in_word = false;
-  int count {0};
-  for (auto c = text.begin(); c != text.end(); c++) {
-    if (std::isspace(*c)) {
-      if (in_word) {
-        in_word = false;
-        count++;
-      }
-    } else {
-      in_word = true;
-    }
-  }
-  if (in_word) count++;
-  return count;
-}
-
 void read_coeffs(pugi::xml_node surf_node, int surf_id, double &c1)
 {
   // Check the given number of coefficients.