diff --git a/CMakeLists.txt b/CMakeLists.txt
index ed9021cf4..5c6636cbc 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -436,6 +436,7 @@ set(LIBOPENMC_CXX_SRC
   src/initialize.cpp
   src/hdf5_interface.cpp
   src/message_passing.cpp
+  src/plot.cpp
   src/random_lcg.cpp
   src/simulation.cpp
   src/state_point.cpp
diff --git a/src/plot.F90 b/src/plot.F90
index 92b8f350d..c2ffa2d39 100644
--- a/src/plot.F90
+++ b/src/plot.F90
@@ -2,8 +2,6 @@ module plot
 
   use, intrinsic :: ISO_C_BINDING
 
-  use hdf5
-
   use constants
   use error,           only: fatal_error, write_message
   use geometry,        only: find_cell, check_cell_overlap
@@ -340,23 +338,46 @@ contains
   subroutine create_voxel(pl)
     type(ObjectPlot), intent(in) :: pl
 
-    integer :: x, y, z      ! voxel location indices
+    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 :: hdf5_err
-    integer, target :: data(pl%pixels(3),pl%pixels(2))
+    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)
-    integer(HSIZE_T) :: dims_slab(3)
-    integer(HSIZE_T) :: offset(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
-    type(c_ptr)       :: f_ptr
+
+    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)
@@ -390,20 +411,8 @@ contains
 
     ! 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(3), pl%pixels(2), pl%pixels(1)]
-    call h5screate_simple_f(3, dims, dspace, hdf5_err)
-    call h5dcreate_f(file_id, "data", H5T_NATIVE_INTEGER, dspace, dset, hdf5_err)
-
-    ! Create another dataspace for 2D array in memory
-    dims_slab(1) = pl%pixels(3)
-    dims_slab(2) = pl%pixels(2)
-    dims_slab(3) = 1
-    call h5screate_simple_f(2, dims_slab(1:2), memspace, hdf5_err)
-
-    ! Initialize offset and get pointer to data
-    offset(:) = 0
-    call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims_slab, hdf5_err)
-    f_ptr = c_loc(data)
+    dims(:) = pl % pixels
+    call voxel_init(file_id, dims, dspace, dset, memspace)
 
     ! move to center of voxels
     ll = ll + vox / TWO
@@ -433,15 +442,10 @@ contains
       p % coord(1) % xyz(3) = ll(3)
 
       ! Write to HDF5 dataset
-      offset(3) = x - 1
-      call h5soffset_simple_f(dspace, offset, hdf5_err)
-      call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, &
-           mem_space_id=memspace, file_space_id=dspace)
+      call voxel_write_slice(x, dspace, dset, memspace, c_loc(data))
     end do
 
-    call h5dclose_f(dset, hdf5_err)
-    call h5sclose_f(dspace, hdf5_err)
-    call h5sclose_f(memspace, hdf5_err)
+    call voxel_finalize(dspace, dset, memspace)
     call file_close(file_id)
 
   end subroutine create_voxel
diff --git a/src/plot.cpp b/src/plot.cpp
new file mode 100644
index 000000000..0dca6d08c
--- /dev/null
+++ b/src/plot.cpp
@@ -0,0 +1,42 @@
+#include "plot.h"
+
+namespace openmc {
+
+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);
+  *dset = H5Dcreate(file_id, "data", H5T_NATIVE_INT, *dspace, H5P_DEFAULT,
+                    H5P_DEFAULT, H5P_DEFAULT);
+
+  // Create dataspace for a slice of the voxel
+  hsize_t dims_slice[2] {dims[1], dims[2]};
+  *memspace = H5Screate_simple(2, dims_slice, nullptr);
+
+  // Select hyperslab in dataspace
+  hsize_t start[3] {0, 0, 0};
+  hsize_t count[3] {1, dims[1], dims[2]};
+  H5Sselect_hyperslab(*dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
+}
+
+
+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};
+  H5Soffset_simple(dspace, offset);
+  H5Dwrite(dset, H5T_NATIVE_INT, memspace, dspace, H5P_DEFAULT, buf);
+}
+
+
+void
+voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace)
+{
+  H5Dclose(dset);
+  H5Sclose(dspace);
+  H5Sclose(memspace);
+}
+
+} // namespace openmc
diff --git a/src/plot.h b/src/plot.h
new file mode 100644
index 000000000..75406b9ac
--- /dev/null
+++ b/src/plot.h
@@ -0,0 +1,15 @@
+#ifndef PLOT_H
+#define PLOT_H
+
+#include "hdf5.h"
+
+namespace openmc {
+
+extern "C" 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,
+                                  hid_t memspace, void* buf);
+extern "C" void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace);
+
+} // namespace openmc
+#endif // PLOT_H