diff --git a/include/openmc/settings.h b/include/openmc/settings.h
index 4b76c3dc59..fc6fea6e7f 100644
--- a/include/openmc/settings.h
+++ b/include/openmc/settings.h
@@ -36,7 +36,7 @@ extern bool particle_restart_run;     //!< particle restart run?
 extern "C" bool photon_transport;     //!< photon transport turned on?
 extern "C" bool reduce_tallies;       //!< reduce tallies at end of batch?
 extern bool res_scat_on;              //!< use resonance upscattering method?
-extern bool restart_run;              //!< restart run?
+extern "C" bool restart_run;          //!< restart run?
 extern "C" bool run_CE;               //!< run with continuous-energy data?
 extern bool source_latest;            //!< write latest source at each batch?
 extern bool source_separate;          //!< write source to separate file?
@@ -57,7 +57,7 @@ extern std::string path_output;           //!< directory where output files are
 extern std::string path_particle_restart; //!< path to a particle restart file
 extern std::string path_source;
 extern std::string path_sourcepoint;      //!< path to a source file
-extern std::string path_statepoint;       //!< path to a statepoint file
+extern "C" std::string path_statepoint;   //!< path to a statepoint file
 
 extern int32_t index_entropy_mesh;  //!< Index of entropy mesh in global mesh array
 extern int32_t index_ufs_mesh;      //!< Index of UFS mesh in global mesh array
diff --git a/openmc/capi/settings.py b/openmc/capi/settings.py
index 115962395d..3dd672354f 100644
--- a/openmc/capi/settings.py
+++ b/openmc/capi/settings.py
@@ -22,6 +22,7 @@ class _Settings(object):
     generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch')
     inactive = _DLLGlobal(c_int32, 'n_inactive')
     particles = _DLLGlobal(c_int64, 'n_particles')
+    restart_run = _DLLGlobal(c_bool, 'restart_run')
     run_CE = _DLLGlobal(c_bool, 'run_CE')
     verbosity = _DLLGlobal(c_int, 'verbosity')
 
@@ -43,6 +44,11 @@ class _Settings(object):
         else:
             raise ValueError('Invalid run mode: {}'.format(mode))
 
+    @property
+    def path_statepoint(self):
+        path = c_char_p.in_dll(_dll, 'path_statepoint').value
+        return path.decode()
+
     @property
     def seed(self):
         return _dll.openmc_get_seed()
diff --git a/openmc/cmfd.py b/openmc/cmfd.py
index 18a5caa1ce..a7acda23a4 100644
--- a/openmc/cmfd.py
+++ b/openmc/cmfd.py
@@ -337,11 +337,11 @@ class CMFDRun(object):
         self._openmc_src = None
         self._sourcecounts = None
         self._weightfactors = None
-        self._entropy = None
-        self._balance = None
-        self._src_cmp = None
-        self._dom = None
-        self._k_cmfd = None
+        self._entropy = []
+        self._balance = []
+        self._src_cmp = []
+        self._dom = []
+        self._k_cmfd = []
         self._resnb = None
         self._time_cmfd = None
         self._time_cmfdbuild = None
@@ -636,8 +636,8 @@ class CMFDRun(object):
             # Configure CMFD parameters and tallies
             self._configure_cmfd()
 
-            # Set up CMFD coremap
-            self._set_coremap()
+            # Initialize all arrays used for CMFD solver
+            self._allocate_cmfd()
 
             # Compute and store array indices used to build cross section
             # arrays
@@ -731,16 +731,24 @@ class CMFDRun(object):
 
     def _configure_cmfd(self):
         """Initialize CMFD parameters and set CMFD input variables"""
-        # Read in cmfd input defined in Python
-        self._read_cmfd_input()
+        # Check if restarting simulation from statepoint file
+        if not openmc.capi.settings.restart_run:
+            # Read in cmfd input defined in Python
+            self._read_cmfd_input()
 
-        # Initialize timers
-        self._time_cmfd = 0.0
-        self._time_cmfdbuild = 0.0
-        self._time_cmfdsolve = 0.0
+            # Set up CMFD coremap
+            self._set_coremap()
 
-        # Initialize all numpy arrays used for CMFD solver
-        self._allocate_cmfd()
+            # Initialize timers
+            self._time_cmfd = 0.0
+            self._time_cmfdbuild = 0.0
+            self._time_cmfdsolve = 0.0
+
+        else:
+            # Reset CMFD parameters from statepoint file
+            # TODO implement reset_cmfd
+            path_statepoint = openmc.capi.settings.path_statepoint
+            sys.exit()
 
     def _read_cmfd_input(self):
         """Sets values of additional instance variables based on user input"""
@@ -803,38 +811,13 @@ class CMFDRun(object):
         nz = self._indices[2]
         ng = self._indices[3]
 
-        # Allocate flux, cross sections and diffusion coefficient
-        self._flux = np.zeros((nx, ny, nz, ng))
-        self._totalxs = np.zeros((nx, ny, nz, ng))
-        self._p1scattxs = np.zeros((nx, ny, nz, ng))
-        self._scattxs = np.zeros((nx, ny, nz, ng, ng))  # Incoming, outgoing
-        self._nfissxs = np.zeros((nx, ny, nz, ng, ng))  # Incoming, outgoing
-        self._diffcof = np.zeros((nx, ny, nz, ng))
-
         # Allocate dtilde and dhat
         self._dtilde = np.zeros((nx, ny, nz, ng, 6))
         self._dhat = np.zeros((nx, ny, nz, ng, 6))
 
         # Allocate dimensions for each mesh cell
         self._hxyz = np.zeros((nx, ny, nz, 3))
-
-        # Allocate surface currents
-        self._current = np.zeros((nx, ny, nz, 12, ng))
-
-        # Allocate source distributions
-        self._cmfd_src = np.zeros((nx, ny, nz, ng))
-        self._openmc_src = np.zeros((nx, ny, nz, ng))
-
-        # Allocate source weight modification variables
-        self._sourcecounts = np.zeros((nx*ny*nz, ng))
-        self._weightfactors = np.ones((nx, ny, nz, ng))
-
-        # Allocate batchwise parameters
-        self._entropy = []
-        self._balance = []
-        self._src_cmp = []
-        self._dom = []
-        self._k_cmfd = []
+        self._hxyz[:,:,:,:] = openmc.capi.meshes[self._mesh_id].width
 
     def _cmfd_init_batch(self):
         """Handles CMFD options at the beginning of each batch"""
@@ -1048,9 +1031,6 @@ class CMFDRun(object):
         ng = self._indices[3]
         n = self._mat_dim
 
-        # Reset CMFD source to 0
-        self._cmfd_src.fill(0.)
-
         # Compute cmfd_src in a vecotorized manner by phi to the spatial
         # indices of the actual problem so that cmfd_flux can be multiplied by
         # nfissxs
@@ -1117,7 +1097,7 @@ class CMFDRun(object):
             ng = self._indices[3]
 
             # Set weight factors to default 1.0
-            self._weightfactors.fill(1.0)
+            self._weightfactors = np.ones((nx, ny, nz, ng))
 
             # Count bank site in mesh and reverse due to egrid structured
             outside = self._count_bank_sites()
@@ -1150,7 +1130,7 @@ class CMFDRun(object):
                 with np.errstate(divide='ignore', invalid='ignore'):
                     self._weightfactors = (np.divide(self._cmfd_src * norm,
                                            sourcecounts, where=div_condition,
-                                           out=np.ones_like(self._cmfd_src), dtype='float16'))
+                                           out=np.ones_like(self._cmfd_src), dtype='float32'))
 
             if not self._feedback:
                 return
@@ -1208,9 +1188,11 @@ class CMFDRun(object):
         bank = openmc.capi.source_bank()
         energy = self._egrid
         sites_outside = np.zeros(1, dtype=bool)
+        nxnynz = np.prod(self._indices[0:3])
         ng = self._indices[3]
 
         outside = np.zeros(1, dtype=bool)
+        self._sourcecounts = np.zeros((nxnynz, ng))
         count = np.zeros(self._sourcecounts.shape)
 
         # Get location and energy of each particle in source bank
@@ -1589,13 +1571,6 @@ class CMFDRun(object):
         nz = self._indices[2]
         ng = self._indices[3]
 
-        # Set flux object and source distribution all to zeros
-        self._flux.fill(0.)
-        self._openmc_src.fill(0.)
-
-        # Set mesh widths
-        self._hxyz[:,:,:,:] = openmc.capi.meshes[self._mesh_id].width
-
         # Reset keff_bal to zero
         self._keff_bal = 0.
 
@@ -1993,7 +1968,7 @@ class CMFDRun(object):
     def _precompute_matrix_indices(self):
         """Computes the indices and row/column data used to populate CMFD CSR
         matrices. These indices are used in _build_loss_matrix and
-        _build_prod_matrix
+        _build_prod_matrix.
 
         """
         # Extract energy group indices
diff --git a/src/state_point.cpp b/src/state_point.cpp
index ab14734033..a849efe10f 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -29,9 +29,6 @@
 
 namespace openmc {
 
-extern "C" void statepoint_write_f(hid_t file_id);
-extern "C" void load_state_point_f(hid_t file_id);
-
 extern "C" int
 openmc_statepoint_write(const char* filename, bool* write_source)
 {