@@ -95,7 +95,7 @@ extern "C" {
   int openmc_sphharm_filter_get_cosine(int32_t index, char cosine[]);
   int openmc_sphharm_filter_set_order(int32_t index, int order);
   int openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]);
-  int openmc_statepoint_write(const char filename[]);
+  int openmc_statepoint_write(const char filename[], bool* write_source);
   int openmc_tally_allocate(int32_t index, const char* type);
   int openmc_tally_get_active(int32_t index, bool* active);
   int openmc_tally_get_estimator(int32_t index, int32_t* estimator);
 
@@ -0,0 +1,17 @@
+#ifndef OPENMC_CONTAINER_UTIL_H
+#define OPENMC_CONTAINER_UTIL_H
+
+#include <algorithm> // for find
+#include <iterator> // for begin, end
+
+namespace openmc {
+
+template<class C, class T>
+inline bool contains(const C& v, const T& x)
+{
+  return std::end(v) != std::find(std::begin(v), std::end(v), x);
+}
+
+}
+
+#endif // OPENMC_CONTAINER_UTIL_H
@@ -7,6 +7,7 @@
 #include <array>
 #include <cstdint>
 #include <string>
+#include <unordered_set>
 #include <vector>
 
 #include "pugixml.hpp"
@@ -76,6 +77,8 @@ extern "C" int res_scat_method;          //!< resonance upscattering method
 extern "C" double res_scat_energy_min;   //!< Min energy in [eV] for res. upscattering
 extern "C" double res_scat_energy_max;   //!< Max energy in [eV] for res. upscattering
 extern "C" int run_mode;                 //!< Run mode (eigenvalue, fixed src, etc.)
+extern std::unordered_set<int> sourcepoint_batch; //!< Batches when source should be written
+extern std::unordered_set<int> statepoint_batch; //!< Batches when state should be written
 extern "C" int temperature_method;       //!< method for choosing temperatures
 extern "C" double temperature_tolerance; //!< Tolerance in [K] on choosing temperatures
 extern "C" double temperature_default;   //!< Default T in [K]
 
@@ -66,6 +66,13 @@ extern "C" void initialize_generation();
 
 extern "C" void initialize_history(Particle* p, int64_t index_source);
 
+//! Finalize a batch
+//!
+//! Handles synchronization and accumulation of tallies, calculation of Shannon
+//! entropy, getting single-batch estimate of keff, and turning on tallies when
+//! appropriate
+extern "C" void finalize_batch();
+
 //! Finalize a fission generation
 extern "C" void finalize_generation();
 
@@ -74,7 +81,6 @@ extern "C" int overall_generation();
 
 #ifdef OPENMC_MPI
 extern "C" void broadcast_results();
-extern "C" void broadcast_triggers();
 #endif
 
 } // namespace openmc
 
@@ -9,6 +9,7 @@
 
 namespace openmc {
 
+void write_source_point(const char* filename);
 extern "C" void write_source_bank(hid_t group_id, Bank* source_bank);
 extern "C" void read_source_bank(hid_t group_id, Bank* source_bank);
 extern "C" void write_tally_results_nr(hid_t file_id);
 
@@ -45,6 +45,7 @@ extern Timer time_bank_sendrecv;
 extern Timer time_finalize;
 extern Timer time_inactive;
 extern Timer time_initialize;
+extern Timer time_tallies;
 extern Timer time_total;
 
 } // namespace openmc
 
@@ -132,8 +132,6 @@ end subroutine openmc_free_bank
     n_batch_interval = 1
     n_lost_particles = 0
     n_particles = -1
-    n_source_points = 0
-    n_state_points = 0
     n_tallies = 0
     output_summary = .true.
     output_tallies = .true.
 
@@ -29,7 +29,7 @@ subroutine cmfd_broadcast(n, buffer) bind(C)
 ! EXECUTE_CMFD runs the CMFD calculation
 !==============================================================================
 
-  subroutine execute_cmfd()
+  subroutine execute_cmfd() bind(C)
 
     use cmfd_data,              only: set_up_cmfd
     use cmfd_solver,            only: cmfd_solver_execute
 
@@ -143,7 +143,7 @@ module cmfd_header
   logical, public :: cmfd_run_adjoint = .false.
 
   ! CMFD run logicals
-  logical, public :: cmfd_on = .false.
+  logical(C_BOOL), public, bind(C) :: cmfd_on = .false.
 
   ! CMFD display info
   character(len=25), public :: cmfd_display = 'balance'
 
@@ -210,99 +210,14 @@ subroutine read_settings_xml_f(root_ptr) bind(C)
 
     integer :: i
     integer :: n
-    integer, allocatable :: temp_int_array(:)
     type(XMLNode) :: root
-    type(XMLNode) :: node_sp
     type(XMLNode) :: node_res_scat
     type(XMLNode) :: node_vol
     type(XMLNode), allocatable :: node_vol_list(:)
 
     ! Get proper XMLNode type given pointer
     root % ptr = root_ptr
 
-    ! Check if the user has specified to write state points
-    if (check_for_node(root, "state_point")) then
-
-      ! Get pointer to state_point node
-      node_sp = root % child("state_point")
-
-      ! Determine number of batches at which to store state points
-      if (check_for_node(node_sp, "batches")) then
-        n_state_points = node_word_count(node_sp, "batches")
-      else
-        n_state_points = 0
-      end if
-
-      if (n_state_points > 0) then
-        ! User gave specific batches to write state points
-        allocate(temp_int_array(n_state_points))
-        call get_node_array(node_sp, "batches", temp_int_array)
-        do i = 1, n_state_points
-          call statepoint_batch % add(temp_int_array(i))
-        end do
-        deallocate(temp_int_array)
-      else
-        ! If neither were specified, write state point at last batch
-        n_state_points = 1
-        call statepoint_batch % add(n_batches)
-      end if
-    else
-      ! If no <state_point> tag was present, by default write state point at
-      ! last batch only
-      n_state_points = 1
-      call statepoint_batch % add(n_batches)
-    end if
-
-    ! Check if the user has specified to write source points
-    if (check_for_node(root, "source_point")) then
-
-      ! Get pointer to source_point node
-      node_sp = root % child("source_point")
-
-      ! Determine number of batches at which to store source points
-      if (check_for_node(node_sp, "batches")) then
-        n_source_points = node_word_count(node_sp, "batches")
-      else
-        n_source_points = 0
-      end if
-
-      if (n_source_points > 0) then
-        ! User gave specific batches to write source points
-        allocate(temp_int_array(n_source_points))
-        call get_node_array(node_sp, "batches", temp_int_array)
-        do i = 1, n_source_points
-          call sourcepoint_batch % add(temp_int_array(i))
-        end do
-        deallocate(temp_int_array)
-      else
-        ! If neither were specified, write source points with state points
-        n_source_points = n_state_points
-        do i = 1, n_state_points
-          call sourcepoint_batch % add(statepoint_batch % get_item(i))
-        end do
-      end if
-    else
-      ! If no <source_point> tag was present, by default we keep source bank in
-      ! statepoint file and write it out at statepoints intervals
-      n_source_points = n_state_points
-      do i = 1, n_state_points
-        call sourcepoint_batch % add(statepoint_batch % get_item(i))
-      end do
-    end if
-
-    ! If source is not seperate and is to be written out in the statepoint file,
-    ! make sure that the sourcepoint batch numbers are contained in the
-    ! statepoint list
-    if (.not. source_separate) then
-      do i = 1, n_source_points
-        if (.not. statepoint_batch % contains(sourcepoint_batch % &
-             get_item(i))) then
-          call fatal_error('Sourcepoint batches are not a subset&
-               & of statepoint batches.')
-        end if
-      end do
-    end if
-
     ! Resonance scattering parameters
     if (check_for_node(root, "resonance_scattering")) then
       node_res_scat = root % child("resonance_scattering")