diff --git a/CMakeLists.txt b/CMakeLists.txt
index bb5c45c42..ed7f3b9f0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -381,6 +381,7 @@ add_library(libopenmc SHARED
   src/tallies/tally_header.F90
   src/tallies/trigger.F90
   src/tallies/trigger_header.F90
+  src/bank.cpp
   src/dagmc.cpp
   src/cell.cpp
   src/cmfd_execute.cpp
diff --git a/include/openmc/bank.h b/include/openmc/bank.h
new file mode 100644
index 000000000..b5dcb4d91
--- /dev/null
+++ b/include/openmc/bank.h
@@ -0,0 +1,47 @@
+#ifndef OPENMC_BANK_H
+#define OPENMC_BANK_H
+
+#include <cstdint>
+#include <vector>
+
+namespace openmc {
+
+struct Bank {
+  double wgt;
+  double xyz[3];
+  double uvw[3];
+  double E;
+  int delayed_group;
+  int particle;
+};
+
+} // namespace openmc
+
+// Without an explicit instantiation of vector<Bank>, the Intel compiler
+// will complain about the threadprivate directive on filter_matches. Note that
+// this has to happen *outside* of the openmc namespace
+extern template class std::vector<openmc::Bank>;
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace simulation {
+
+extern "C" int64_t n_bank;
+
+extern std::vector<Bank> source_bank;
+extern std::vector<Bank> fission_bank;
+#ifdef _OPENMP
+extern std::vector<Bank> master_fission_bank;
+#endif
+
+#pragma omp threadprivate(fission_bank, n_bank)
+
+} // namespace simulation
+
+} // namespace openmc
+
+#endif // OPENMC_BANK_H
diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index fdea48e9e..50ae72000 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -6,9 +6,11 @@
 #include <stddef.h>
 
 #ifdef __cplusplus
+#include "openmc/bank.h"
 extern "C" {
-#endif
-
+  int openmc_fission_bank(openmc::Bank** ptr, int64_t* n);
+  int openmc_source_bank(openmc::Bank** ptr, int64_t* n);
+#else
   struct Bank {
     double wgt;
     double xyz[3];
@@ -18,6 +20,10 @@ extern "C" {
     int particle;
   };
 
+  int openmc_fission_bank(struct Bank** ptr, int64_t* n);
+  int openmc_source_bank(struct Bank** ptr, int64_t* n);
+#endif
+
   int openmc_calculate_volumes();
   int openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n);
   int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n);
@@ -38,7 +44,6 @@ extern "C" {
   int openmc_filter_set_type(int32_t index, const char* type);
   int openmc_finalize();
   int openmc_find_cell(double* xyz, int32_t* index, int32_t* instance);
-  int openmc_fission_bank(struct Bank** ptr, int64_t* n);
   int openmc_get_cell_index(int32_t id, int32_t* index);
   int openmc_get_filter_index(int32_t id, int32_t* index);
   void openmc_get_filter_next_id(int32_t* id);
@@ -85,7 +90,6 @@ extern "C" {
   void openmc_set_seed(int64_t new_seed);
   int openmc_simulation_finalize();
   int openmc_simulation_init();
-  int openmc_source_bank(struct Bank** ptr, int64_t* n);
   int openmc_spatial_legendre_filter_get_order(int32_t index, int* order);
   int openmc_spatial_legendre_filter_get_params(int32_t index, int* axis, double* min, double* max);
   int openmc_spatial_legendre_filter_set_order(int32_t index, int order);
diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h
index 6d918eb40..f7e93ec2b 100644
--- a/include/openmc/eigenvalue.h
+++ b/include/openmc/eigenvalue.h
@@ -25,9 +25,6 @@ extern std::array<double, 2> k_sum; //!< Used to reduce sum and sum_sq
 extern std::vector<double> entropy; //!< Shannon entropy at each generation
 extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
 
-extern "C" int64_t n_bank;
-#pragma omp threadprivate(n_bank)
-
 } // namespace simulation
 
 //==============================================================================
@@ -35,15 +32,22 @@ extern "C" int64_t n_bank;
 //==============================================================================
 
 //! Collect/normalize the tracklength keff from each process
-extern "C" void calculate_generation_keff();
+void calculate_generation_keff();
 
 //! Calculate mean/standard deviation of keff during active generations
 //!
 //! This function sets the global variables keff and keff_std which represent
 //! the mean and standard deviation of the mean of k-effective over active
 //! generations. It also broadcasts the value from the master process.
-extern "C" void calculate_average_keff();
+void calculate_average_keff();
 
+#ifdef _OPENMP
+//! Join threadprivate fission banks into a single fission bank
+//!
+//! Note that this operation is necessarily sequential to preserve the order of
+//! the bank when using varying numbers of threads.
+void join_bank_from_threads();
+#endif
 
 //! Calculates a minimum variance estimate of k-effective
 //!
@@ -60,16 +64,16 @@ extern "C" void calculate_average_keff();
 extern "C" int openmc_get_keff(double* k_combined);
 
 //! Sample/redistribute source sites from accumulated fission sites
-extern "C" void synchronize_bank();
+void synchronize_bank();
 
 //! Calculates the Shannon entropy of the fission source distribution to assess
 //! source convergence
-extern "C" void shannon_entropy();
+void shannon_entropy();
 
 //! Determines the source fraction in each UFS mesh cell and reweights the
 //! source bank so that the sum of the weights is equal to n_particles. The
 //! 'source_frac' variable is used later to bias the production of fission sites
-extern "C" void ufs_count_sites();
+void ufs_count_sites();
 
 //! Get UFS weight corresponding to particle's location
 extern "C" double ufs_get_weight(const Particle* p);
diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h
index 2002e95fe..cba0a2b71 100644
--- a/include/openmc/simulation.h
+++ b/include/openmc/simulation.h
@@ -56,6 +56,9 @@ extern "C" int thread_id;  //!< ID of a given thread
 // Functions
 //==============================================================================
 
+//! Allocate space for source and fission banks
+void allocate_banks();
+
 //! Determine number of particles to transport per process
 void calculate_work();
 
diff --git a/include/openmc/source.h b/include/openmc/source.h
index 45f4896b4..d9f5cbefa 100644
--- a/include/openmc/source.h
+++ b/include/openmc/source.h
@@ -9,9 +9,9 @@
 
 #include "pugixml.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/distribution_multi.h"
 #include "openmc/distribution_spatial.h"
-#include "openmc/capi.h"
 #include "openmc/particle.h"
 
 namespace openmc {
@@ -62,7 +62,7 @@ extern "C" void initialize_source();
 //! Sample a site from all external source distributions in proportion to their
 //! source strength
 //! \return Sampled source site
-extern "C" Bank sample_external_source();
+Bank sample_external_source();
 
 //! Fill source bank at end of generation for fixed source simulations
 void fill_source_bank_fixedsource();
diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h
index 6b019e7b8..0c1db86f3 100644
--- a/include/openmc/state_point.h
+++ b/include/openmc/state_point.h
@@ -10,8 +10,8 @@
 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_source_bank(hid_t group_id);
+extern "C" void read_source_bank(hid_t group_id);
 extern "C" void write_tally_results_nr(hid_t file_id);
 extern "C" void restart_set_keff();
 
diff --git a/src/api.F90 b/src/api.F90
index efcfa16a6..d7fedcbe3 100644
--- a/src/api.F90
+++ b/src/api.F90
@@ -82,6 +82,9 @@ contains
 
       subroutine free_memory_settings() bind(C)
       end subroutine free_memory_settings
+
+      subroutine free_memory_bank() bind(C)
+      end subroutine free_memory_bank
     end interface
 
     call free_memory_geometry()
diff --git a/src/bank.cpp b/src/bank.cpp
new file mode 100644
index 000000000..89c493bab
--- /dev/null
+++ b/src/bank.cpp
@@ -0,0 +1,107 @@
+#include "openmc/bank.h"
+
+#include "openmc/capi.h"
+#include "openmc/error.h"
+
+#include <cstdint>
+
+// Explicit template instantiation definition
+template class std::vector<openmc::Bank>;
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace simulation {
+
+int64_t n_bank;
+
+std::vector<Bank> source_bank;
+std::vector<Bank> fission_bank;
+#ifdef _OPENMP
+std::vector<Bank> master_fission_bank;
+#endif
+
+} // namespace simulation
+
+//==============================================================================
+// C API
+//==============================================================================
+
+extern "C" int openmc_source_bank(Bank** ptr, int64_t* n)
+{
+  if (simulation::source_bank.size() == 0) {
+    set_errmsg("Source bank has not been allocated.");
+    return OPENMC_E_ALLOCATE;
+  } else {
+    *ptr = simulation::source_bank.data();
+    *n = simulation::source_bank.size();
+    return 0;
+  }
+}
+
+extern "C" int openmc_fission_bank(Bank** ptr, int64_t* n)
+{
+  if (simulation::fission_bank.size() == 0) {
+    set_errmsg("Fission bank has not been allocated.");
+    return OPENMC_E_ALLOCATE;
+  } else {
+    *ptr = simulation::fission_bank.data();
+    *n = simulation::fission_bank.size();
+    return 0;
+  }
+}
+
+//==============================================================================
+// Fortran compatibility
+//==============================================================================
+
+extern "C" void free_memory_bank()
+{
+  simulation::source_bank.clear();
+#pragma omp parallel
+  {
+    simulation::fission_bank.clear();
+  }
+#ifdef _OPENMP
+  simulation::master_fission_bank.clear();
+#endif
+}
+
+extern "C" int fission_bank_delayed_group(int64_t i) {
+  return simulation::fission_bank[i-1].delayed_group;
+}
+
+extern "C" double fission_bank_E(int64_t i) {
+  return simulation::fission_bank[i-1].E;
+}
+
+extern "C" double fission_bank_wgt(int64_t i) {
+  return simulation::fission_bank[i-1].wgt;
+}
+
+extern "C" void source_bank_xyz(int64_t i, double* xyz)
+{
+  xyz[0] = simulation::source_bank[i-1].xyz[0];
+  xyz[1] = simulation::source_bank[i-1].xyz[1];
+  xyz[2] = simulation::source_bank[i-1].xyz[2];
+}
+
+extern "C" double source_bank_E(int64_t i)
+{
+  return simulation::source_bank[i-1].E;
+}
+
+extern "C" double source_bank_wgt(int64_t i)
+{
+  return simulation::source_bank[i-1].wgt;
+}
+
+extern "C" void source_bank_set_wgt(int64_t i, double wgt)
+{
+  simulation::source_bank[i-1].wgt = wgt;
+}
+
+} // namespace openmc
diff --git a/src/bank_header.F90 b/src/bank_header.F90
index c39bf8b78..992d357fb 100644
--- a/src/bank_header.F90
+++ b/src/bank_header.F90
@@ -2,8 +2,6 @@ module bank_header
 
   use, intrinsic :: ISO_C_BINDING
 
-  use error, only: E_ALLOCATE, set_errmsg
-
   implicit none
 
 !===============================================================================
@@ -21,70 +19,58 @@ module bank_header
     integer(C_INT) :: particle      ! particle type (neutron, photon, etc.)
   end type Bank
 
-  ! Source and fission bank
-  type(Bank), allocatable, target :: source_bank(:)
-  type(Bank), allocatable, target :: fission_bank(:)
-#ifdef _OPENMP
-  type(Bank), allocatable, target :: master_fission_bank(:)
-#endif
-
   integer(C_INT64_T), bind(C) :: n_bank       ! # of sites in fission bank
+!$omp threadprivate(n_bank)
 
-!$omp threadprivate(fission_bank, n_bank)
+  interface
+    function openmc_fission_bank(ptr, n) result(err) bind(C)
+      import C_PTR, C_INT64_T, C_INT
+      type(C_PTR), intent(out) :: ptr
+      integer(C_INT64_T), intent(out) :: n
+      integer(C_INT) :: err
+    end function
 
-contains
+    function fission_bank_delayed_group(i) result(g) bind(C)
+      import C_INT64_T, C_INT
+      integer(C_INT64_T), value :: i
+      integer(C_INT) :: g
+    end function
 
-!===============================================================================
-! FREE_MEMORY_BANK deallocates global arrays defined in this module
-!===============================================================================
+    function fission_bank_E(i) result(E) bind(C, name='fission_bank_E')
+      import C_INT64_T, C_DOUBLE
+      integer(C_INT64_T), value :: i
+      real(C_DOUBLE) :: E
+    end function
 
-  subroutine free_memory_bank()
+    function fission_bank_wgt(i) result(wgt) bind(C)
+      import C_INT64_T, C_DOUBLE
+      integer(C_INT64_T), value :: i
+      real(C_DOUBLE) :: wgt
+    end function
 
-    ! Deallocate fission and source bank and entropy
-!$omp parallel
-    if (allocated(fission_bank)) deallocate(fission_bank)
-!$omp end parallel
-#ifdef _OPENMP
-    if (allocated(master_fission_bank)) deallocate(master_fission_bank)
-#endif
-    if (allocated(source_bank)) deallocate(source_bank)
+    subroutine source_bank_xyz(i, xyz) bind(C)
+      import C_INT64_T, C_DOUBLE
+      integer(C_INT64_T), value :: i
+      real(C_DOUBLE), intent(in) :: xyz(*)
+    end subroutine
 
-  end subroutine free_memory_bank
+    function source_bank_E(i) result(E) bind(C, name='source_bank_E')
+      import C_INT64_T, C_DOUBLE
+      integer(C_INT64_T), value :: i
+      real(C_DOUBLE) :: E
+    end function
 
-!===============================================================================
-!                               C API FUNCTIONS
-!===============================================================================
+    function source_bank_wgt(i) result(wgt) bind(C)
+      import C_INT64_T, C_DOUBLE
+      integer(C_INT64_T), value :: i
+      real(C_DOUBLE) :: wgt
+    end function
 
-  function openmc_source_bank(ptr, n) result(err) bind(C)
-    ! Return a pointer to the source bank
-    type(C_PTR), intent(out) :: ptr
-    integer(C_INT64_T), intent(out) :: n
-    integer(C_INT) :: err
-
-    if (.not. allocated(source_bank)) then
-      err = E_ALLOCATE
-      call set_errmsg("Source bank has not been allocated.")
-    else
-      err = 0
-      ptr = C_LOC(source_bank)
-      n = size(source_bank)
-    end if
-  end function openmc_source_bank
-
-  function openmc_fission_bank(ptr, n) result(err) bind(C)
-    ! Return a pointer to the source bank
-    type(C_PTR), intent(out) :: ptr
-    integer(C_INT64_T), intent(out) :: n
-    integer(C_INT) :: err
-
-    if (.not. allocated(fission_bank)) then
-      err = E_ALLOCATE
-      call set_errmsg("Fission bank has not been allocated.")
-    else
-      err = 0
-      ptr = C_LOC(fission_bank)
-      n = size(fission_bank)
-    end if
-  end function openmc_fission_bank
+    subroutine source_bank_set_wgt(i, wgt) bind(C)
+      import C_INT64_T, C_DOUBLE
+      integer(C_INT64_T), value :: i
+      real(C_DOUBLE), value :: wgt
+    end subroutine
+  end interface
 
 end module bank_header
diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90
index 681d08c5c..4b97b5e32 100644
--- a/src/cmfd_execute.F90
+++ b/src/cmfd_execute.F90
@@ -218,7 +218,7 @@ contains
   subroutine cmfd_reweight(new_weights)
 
     use algorithm,   only: binary_search
-    use bank_header, only: source_bank
+    use bank_header
     use constants,   only: ZERO, ONE
     use error,       only: warning, fatal_error
     use message_passing
@@ -230,13 +230,14 @@ contains
     integer :: ny       ! maximum number of cells in y direction
     integer :: nz       ! maximum number of cells in z direction
     integer(C_INT) :: ng       ! maximum number of energy groups
-    integer :: i        ! iteration counter
+    integer(8) :: i        ! iteration counter
     integer :: g        ! index for group
     integer :: ijk(3)   ! spatial bin location
     integer :: e_bin    ! energy bin of source particle
     integer :: mesh_bin ! mesh bin of soruce particle
     integer :: n_groups ! number of energy groups
     real(8) :: norm     ! normalization factor
+    real(C_DOUBLE) :: xyz(3)
     logical(C_BOOL) :: outside  ! any source sites outside mesh
     logical :: in_mesh  ! source site is inside mesh
 
@@ -309,21 +310,22 @@ contains
     end if
 
     ! begin loop over source bank
-    do i = 1, int(work,4)
+    do i = 1, work
 
       ! Determine spatial bin
-      call cmfd_mesh % get_indices(source_bank(i) % xyz, ijk, in_mesh)
+      call source_bank_xyz(i, xyz)
+      call cmfd_mesh % get_indices(xyz, ijk, in_mesh)
 
       ! Determine energy bin
       n_groups = size(cmfd % egrid) - 1
-      if (source_bank(i) % E < cmfd % egrid(1)) then
+      if (source_bank_E(i) < cmfd % egrid(1)) then
         e_bin = 1
         if (master) call warning('Source pt below energy grid')
-      elseif (source_bank(i) % E > cmfd % egrid(n_groups + 1)) then
+      elseif (source_bank_E(i) > cmfd % egrid(n_groups + 1)) then
         e_bin = n_groups
         if (master) call warning('Source pt above energy grid')
       else
-        e_bin = binary_search(cmfd % egrid, n_groups + 1, source_bank(i) % E)
+        e_bin = binary_search(cmfd % egrid, n_groups + 1, source_bank_E(i))
       end if
 
       ! Reverese energy bin (lowest grp is highest energy bin)
@@ -335,8 +337,8 @@ contains
       end if
 
       ! Reweight particle
-      source_bank(i) % wgt = source_bank(i) % wgt * &
-           cmfd % weightfactors(e_bin, ijk(1), ijk(2), ijk(3))
+      call source_bank_set_wgt(i, source_bank_wgt(i) * &
+           cmfd % weightfactors(e_bin, ijk(1), ijk(2), ijk(3)))
     end do
 
   end subroutine cmfd_reweight
diff --git a/src/cmfd_execute.cpp b/src/cmfd_execute.cpp
index 0af584c1b..28561b00a 100644
--- a/src/cmfd_execute.cpp
+++ b/src/cmfd_execute.cpp
@@ -5,6 +5,7 @@
 #include "xtensor/xarray.hpp"
 #include "xtensor/xio.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/mesh.h"
 #include "openmc/message_passing.h"
@@ -18,14 +19,10 @@ extern "C" void
 cmfd_populate_sourcecounts(int n_energy, const double* energies,
   double* source_counts, bool* outside)
 {
-  // Get pointer to source bank
-  Bank* source_bank;
-  int64_t n;
-  openmc_source_bank(&source_bank, &n);
-
   // Get source counts in each mesh bin / energy bin
   auto& m = model::meshes.at(settings::index_cmfd_mesh);
-  xt::xarray<double> counts = m->count_sites(simulation::work, source_bank, n_energy, energies, outside);
+  xt::xarray<double> counts = m->count_sites(simulation::work,
+    simulation::source_bank.data(), n_energy, energies, outside);
 
   // Copy data from the xarray into the source counts array
   std::copy(counts.begin(), counts.end(), source_counts);
diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90
index bdbbd1c2b..4de3ae8a6 100644
--- a/src/eigenvalue.F90
+++ b/src/eigenvalue.F90
@@ -2,9 +2,6 @@ module eigenvalue
 
   use, intrinsic :: ISO_C_BINDING
 
-  use bank_header
-  use simulation_header
-
   implicit none
 
   interface
@@ -15,49 +12,4 @@ module eigenvalue
     end function
   end interface
 
-contains
-
-#ifdef _OPENMP
-!===============================================================================
-! JOIN_BANK_FROM_THREADS joins threadprivate fission banks into a single fission
-! bank that can be sampled. Note that this operation is necessarily sequential
-! to preserve the order of the bank when using varying numbers of threads.
-!===============================================================================
-
-  subroutine join_bank_from_threads() bind(C)
-
-    integer(8) :: total ! total number of fission bank sites
-    integer    :: i     ! loop index for threads
-
-    ! Initialize the total number of fission bank sites
-    total = 0
-
-!$omp parallel
-
-    ! Copy thread fission bank sites to one shared copy
-!$omp do ordered schedule(static)
-    do i = 1, n_threads
-!$omp ordered
-      master_fission_bank(total+1:total+n_bank) = fission_bank(1:n_bank)
-      total = total + n_bank
-!$omp end ordered
-    end do
-!$omp end do
-
-    ! Make sure all threads have made it to this point
-!$omp barrier
-
-    ! Now copy the shared fission bank sites back to the master thread's copy.
-    if (thread_id == 0) then
-      n_bank = total
-      fission_bank(1:n_bank) = master_fission_bank(1:n_bank)
-    else
-      n_bank = 0
-    end if
-
-!$omp end parallel
-
-  end subroutine join_bank_from_threads
-#endif
-
 end module eigenvalue
diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp
index 64e0c4e42..3aa140a3c 100644
--- a/src/eigenvalue.cpp
+++ b/src/eigenvalue.cpp
@@ -5,6 +5,7 @@
 #include "xtensor/xtensor.hpp"
 #include "xtensor/xview.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/constants.h"
 #include "openmc/error.h"
@@ -69,13 +70,6 @@ void synchronize_bank()
 {
   simulation::time_bank.start();
 
-  // Get pointers to source/fission bank
-  Bank* source_bank;
-  Bank* fission_bank;
-  int64_t n;
-  openmc_source_bank(&source_bank, &n);
-  openmc_fission_bank(&fission_bank, &n);
-
   // In order to properly understand the fission bank algorithm, you need to
   // think of the fission and source bank as being one global array divided
   // over multiple processors. At the start, each processor has a random amount
@@ -147,14 +141,14 @@ void synchronize_bank()
     // and the remaining 100 would be randomly sampled.
     if (total < settings::n_particles) {
       for (int64_t j = 1; j <= settings::n_particles / total; ++j) {
-        temp_sites[index_temp] = fission_bank[i];
+        temp_sites[index_temp] = simulation::fission_bank[i];
         ++index_temp;
       }
     }
 
     // Randomly sample sites needed
     if (prn() < p_sample) {
-      temp_sites[index_temp] = fission_bank[i];
+      temp_sites[index_temp] = simulation::fission_bank[i];
       ++index_temp;
     }
   }
@@ -195,7 +189,8 @@ void synchronize_bank()
       // fission bank
       sites_needed = settings::n_particles - finish;
       for (int i = 0; i < sites_needed; ++i) {
-        temp_sites[index_temp] = fission_bank[simulation::n_bank - sites_needed + i];
+        int i_bank = simulation::n_bank - sites_needed + i;
+        temp_sites[index_temp] = simulation::fission_bank[i_bank];
         ++index_temp;
       }
     }
@@ -300,7 +295,8 @@ void synchronize_bank()
   MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE);
 
 #else
-  std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles, source_bank);
+  std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles,
+    simulation::source_bank.begin());
 #endif
 
   simulation::time_bank_sendrecv.stop();
@@ -347,6 +343,46 @@ void calculate_average_keff()
   }
 }
 
+#ifdef _OPENMP
+void join_bank_from_threads()
+{
+  // Initialize the total number of fission bank sites
+  int64_t total = 0;
+
+#pragma omp parallel
+  {
+    // Copy thread fission bank sites to one shared copy
+#pragma omp for ordered schedule(static)
+    for (int i = 0; i < simulation::n_threads; ++i) {
+#pragma omp ordered
+      {
+        std::copy(
+          &simulation::fission_bank[0],
+          &simulation::fission_bank[0] + simulation::n_bank,
+          &simulation::master_fission_bank[total]
+        );
+        total += simulation::n_bank;
+      }
+    }
+
+    // Make sure all threads have made it to this point
+#pragma omp barrier
+
+    // Now copy the shared fission bank sites back to the master thread's copy.
+    if (simulation::thread_id == 0) {
+      simulation::n_bank = total;
+      std::copy(
+        &simulation::master_fission_bank[0],
+        &simulation::master_fission_bank[0] + simulation::n_bank,
+        &simulation::fission_bank[0]
+      );
+    } else {
+      simulation::n_bank = 0;
+    }
+  }
+}
+#endif
+
 int openmc_get_keff(double* k_combined)
 {
   k_combined[0] = 0.0;
@@ -499,15 +535,10 @@ void shannon_entropy()
   // Get pointer to entropy mesh
   auto& m = model::meshes[settings::index_entropy_mesh];
 
-  // Get pointer to fission bank
-  Bank* fission_bank;
-  int64_t n;
-  openmc_fission_bank(&fission_bank, &n);
-
   // Get source weight in each mesh bin
   bool sites_outside;
-  xt::xtensor<double, 1> p = m->count_sites(
-    simulation::n_bank, fission_bank, 0, nullptr, &sites_outside);
+  xt::xtensor<double, 1> p = m->count_sites(simulation::n_bank,
+    simulation::fission_bank.data(), 0, nullptr, &sites_outside);
 
   // display warning message if there were sites outside entropy box
   if (sites_outside) {
@@ -544,15 +575,10 @@ void ufs_count_sites()
     s = m->volume_frac_;
 
   } else {
-    // Get pointer to source bank
-    Bank* source_bank;
-    int64_t n;
-    openmc_source_bank(&source_bank, &n);
-
     // count number of source sites in each ufs mesh cell
     bool sites_outside;
-    simulation::source_frac = m->count_sites(simulation::work, source_bank, 0, nullptr,
-      &sites_outside);
+    simulation::source_frac = m->count_sites(simulation::work,
+      simulation::source_bank.data(), 0, nullptr, &sites_outside);
 
     // Check for sites outside of the mesh
     if (mpi::master && sites_outside) {
@@ -572,7 +598,7 @@ void ufs_count_sites()
     // Since the total starting weight is not equal to n_particles, we need to
     // renormalize the weight of the source sites
     for (int i = 0; i < simulation::work; ++i) {
-      source_bank[i].wgt *= settings::n_particles / total;
+      simulation::source_bank[i].wgt *= settings::n_particles / total;
     }
   }
 }
diff --git a/src/particle.cpp b/src/particle.cpp
index f87ff3066..c8a753801 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -3,6 +3,7 @@
 #include <algorithm>
 #include <sstream>
 
+#include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/constants.h"
 #include "openmc/error.h"
@@ -184,17 +185,12 @@ Particle::write_restart() const
     write_dataset(file_id, "id", id);
     write_dataset(file_id, "type", type);
 
-    // Get pointer to source bank
-    Bank* src;
-    int64_t n;
-    openmc_source_bank(&src, &n);
-
     int64_t i = simulation::current_work;
-    write_dataset(file_id, "weight", src[i-1].wgt);
-    write_dataset(file_id, "energy", src[i-1].E);
+    write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt);
+    write_dataset(file_id, "energy", simulation::source_bank[i-1].E);
     hsize_t dims[] {3};
-    write_double(file_id, 1, dims, "xyz", src[i-1].xyz, false);
-    write_double(file_id, 1, dims, "uvw", src[i-1].uvw, false);
+    write_double(file_id, 1, dims, "xyz", simulation::source_bank[i-1].xyz, false);
+    write_double(file_id, 1, dims, "uvw", simulation::source_bank[i-1].uvw, false);
 
     // Close file
     file_close(file_id);
diff --git a/src/physics.F90 b/src/physics.F90
index 29a806495..99125915e 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -1,6 +1,7 @@
 module physics
 
   use algorithm,              only: binary_search
+  use bank_header
   use constants
   use endf,                   only: reaction_name
   use error,                  only: fatal_error, warning, write_message
@@ -86,7 +87,11 @@ contains
     integer :: i_nuclide    ! index in nuclides array
     integer :: i_nuc_mat    ! index in material's nuclides array
     integer :: i_reaction   ! index in nuc % reactions array
+    integer(C_INT) :: err
+    integer(C_INT64_T) :: n
     type(Nuclide), pointer :: nuc
+    type(C_PTR) :: ptr
+    type(Bank), pointer :: fission_bank(:)
 
     call sample_nuclide(p, 'total  ', i_nuclide, i_nuc_mat)
 
@@ -103,6 +108,10 @@ contains
 
     if (nuc % fissionable) then
       if (run_mode == MODE_EIGENVALUE) then
+        ! Get fission bank pointer
+        err = openmc_fission_bank(ptr, n)
+        call c_f_pointer(ptr, fission_bank, [n])
+
         call sample_fission(i_nuclide, p % E, i_reaction)
         call create_fission_sites(p, i_nuclide, i_reaction, fission_bank, n_bank)
       elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index a4795f82a..202e5abe3 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -5,6 +5,7 @@
 
 #include "xtensor/xarray.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/constants.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
@@ -47,12 +48,8 @@ sample_reaction(Particle* p, const MaterialMacroXS* material_xs)
 
   if (model::materials[p->material - 1]->fissionable) {
     if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-      Bank* result_bank;
-      int64_t result_bank_size;
-      // Get pointer to fission bank from Fortran side
-      openmc_fission_bank(&result_bank, &result_bank_size);
-      create_fission_sites(p, result_bank, &simulation::n_bank, result_bank_size,
-                           material_xs);
+      create_fission_sites(p, simulation::fission_bank.data(), &simulation::n_bank,
+                           simulation::fission_bank.size(), material_xs);
     } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
                (settings::create_fission_neutrons)) {
       create_fission_sites(p, p->secondary_bank, &(p->n_secondary),
diff --git a/src/simulation.F90 b/src/simulation.F90
index 3eabbfecb..f00b58ed1 100644
--- a/src/simulation.F90
+++ b/src/simulation.F90
@@ -2,20 +2,9 @@ module simulation
 
   use, intrinsic :: ISO_C_BINDING
 
-#ifdef _OPENMP
-  use omp_lib
-#endif
-
-  use bank_header,     only: source_bank
-  use constants,       only: ZERO
-  use error,           only: fatal_error
   use material_header, only: n_materials, materials
-  use message_passing
   use nuclide_header,  only: micro_xs, n_nuclides
   use photon_header,   only: micro_photon_xs, n_elements
-  use settings
-  use simulation_header
-  use tally_header
   use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer
 
   implicit none
@@ -69,56 +58,4 @@ contains
 
   end subroutine
 
-!===============================================================================
-! ALLOCATE_BANKS allocates memory for the fission and source banks
-!===============================================================================
-
-  subroutine allocate_banks() bind(C)
-
-    integer :: alloc_err  ! allocation error code
-
-    ! Allocate source bank
-    if (allocated(source_bank)) deallocate(source_bank)
-    allocate(source_bank(work), STAT=alloc_err)
-
-    ! Check for allocation errors
-    if (alloc_err /= 0) then
-      call fatal_error("Failed to allocate source bank.")
-    end if
-
-    if (run_mode == MODE_EIGENVALUE) then
-
-#ifdef _OPENMP
-      ! If OpenMP is being used, each thread needs its own private fission
-      ! bank. Since the private fission banks need to be combined at the end of
-      ! a generation, there is also a 'master_fission_bank' that is used to
-      ! collect the sites from each thread.
-
-      n_threads = omp_get_max_threads()
-
-!$omp parallel
-      thread_id = omp_get_thread_num()
-
-      if (allocated(fission_bank)) deallocate(fission_bank)
-      if (thread_id == 0) then
-        allocate(fission_bank(3*work))
-      else
-        allocate(fission_bank(3*work/n_threads))
-      end if
-!$omp end parallel
-      if (allocated(master_fission_bank)) deallocate(master_fission_bank)
-      allocate(master_fission_bank(3*work), STAT=alloc_err)
-#else
-      if (allocated(fission_bank)) deallocate(fission_bank)
-      allocate(fission_bank(3*work), STAT=alloc_err)
-#endif
-
-      ! Check for allocation errors
-      if (alloc_err /= 0) then
-        call fatal_error("Failed to allocate fission bank.")
-      end if
-    end if
-
-  end subroutine allocate_banks
-
 end module simulation
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 50bf12521..0cfb99a4a 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -1,5 +1,6 @@
 #include "openmc/simulation.h"
 
+#include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/container_util.h"
 #include "openmc/eigenvalue.h"
@@ -15,6 +16,8 @@
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
 
+#include <omp.h>
+
 #include <algorithm>
 #include <string>
 
@@ -24,14 +27,12 @@ namespace openmc {
 extern "C" bool cmfd_on;
 
 extern "C" void accumulate_tallies();
-extern "C" void allocate_banks();
 extern "C" void allocate_tally_results();
 extern "C" void check_triggers();
 extern "C" void cmfd_init_batch();
 extern "C" void cmfd_tally_init();
 extern "C" void execute_cmfd();
 extern "C" void init_tally_routines();
-extern "C" void join_bank_from_threads();
 extern "C" void load_state_point();
 extern "C" void print_batch_keff();
 extern "C" void print_columns();
@@ -280,6 +281,36 @@ int thread_id;  //!< ID of a given thread
 // Non-member functions
 //==============================================================================
 
+void allocate_banks()
+{
+  // Allocate source bank
+  simulation::source_bank.resize(simulation::work);
+
+  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+#ifdef _OPENMP
+    // If OpenMP is being used, each thread needs its own private fission
+    // bank. Since the private fission banks need to be combined at the end of
+    // a generation, there is also a 'master_fission_bank' that is used to
+    // collect the sites from each thread.
+
+    simulation::n_threads = omp_get_max_threads();
+
+#pragma omp parallel
+    {
+      simulation::thread_id = omp_get_thread_num();
+      if (simulation::thread_id == 0) {
+        simulation::fission_bank.resize(3*simulation::work);
+      } else {
+        simulation::fission_bank.resize(3*simulation::work / simulation::n_threads);
+      }
+    }
+    simulation::master_fission_bank.resize(3*simulation::work);
+#else
+    simulation::fission_bank.resize(3*simulation::work);
+#endif
+  }
+}
+
 void initialize_batch()
 {
   // Increment current batch
@@ -453,13 +484,8 @@ void finalize_generation()
 
 void initialize_history(Particle* p, int64_t index_source)
 {
-  // Get pointer to source bank
-  Bank* source_bank;
-  int64_t n;
-  openmc_source_bank(&source_bank, &n);
-
   // set defaults
-  p->from_source(&source_bank[index_source - 1]);
+  p->from_source(&simulation::source_bank[index_source - 1]);
 
   // set identifier for particle
   p->id = simulation::work_index[mpi::rank] + index_source;
diff --git a/src/source.cpp b/src/source.cpp
index 6d071c0c7..bfa119823 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -5,6 +5,7 @@
 
 #include "xtensor/xadapt.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/cell.h"
 #include "openmc/error.h"
 #include "openmc/file_utils.h"
@@ -242,11 +243,6 @@ void initialize_source()
 {
   write_message("Initializing source particles...", 5);
 
-  // Get pointer to source bank
-  Bank* source_bank;
-  int64_t n;
-  openmc_source_bank(&source_bank, &n);
-
   if (settings::path_source != "") {
     // Read the source from a binary file instead of sampling from some
     // assumed source distribution
@@ -268,7 +264,7 @@ void initialize_source()
     }
 
     // Read in the source bank
-    read_source_bank(file_id, source_bank);
+    read_source_bank(file_id);
 
     // Close file
     file_close(file_id);
@@ -282,7 +278,7 @@ void initialize_source()
       set_particle_seed(id);
 
       // sample external source distribution
-      source_bank[i] = sample_external_source();
+      simulation::source_bank[i] = sample_external_source();
     }
   }
 
@@ -291,7 +287,7 @@ void initialize_source()
     write_message("Writing out initial source...", 5);
     std::string filename = settings::path_output + "initial_source.h5";
     hid_t file_id = file_open(filename, 'w', true);
-    write_source_bank(file_id, source_bank);
+    write_source_bank(file_id);
     file_close(file_id);
   }
 }
@@ -354,11 +350,6 @@ extern "C" double total_source_strength()
 void fill_source_bank_fixedsource()
 {
   if (settings::path_source.empty()) {
-    // Get pointer to source bank
-    Bank* source_bank;
-    int64_t n;
-    openmc_source_bank(&source_bank, &n);
-
     for (int64_t i = 0; i < simulation::work; ++i) {
       // initialize random number seed
       int64_t id = (simulation::total_gen + overall_generation()) *
@@ -366,7 +357,7 @@ void fill_source_bank_fixedsource()
       set_particle_seed(id);
 
       // sample external source distribution
-      source_bank[i] = sample_external_source();
+      simulation::source_bank[i] = sample_external_source();
     }
   }
 }
diff --git a/src/state_point.F90 b/src/state_point.F90
index 2692856dc..9172f1fba 100644
--- a/src/state_point.F90
+++ b/src/state_point.F90
@@ -35,16 +35,14 @@ module state_point
   implicit none
 
   interface
-    subroutine write_source_bank(group_id, bank_) bind(C)
-      import HID_T, C_INT64_T, Bank
+    subroutine write_source_bank(group_id) bind(C)
+      import HID_T
       integer(HID_T), value :: group_id
-      type(Bank), intent(in) :: bank_(*)
     end subroutine write_source_bank
 
-    subroutine read_source_bank(group_id, bank_) bind(C)
-      import HID_T, C_INT64_T, Bank
+    subroutine read_source_bank(group_id) bind(C)
+      import HID_T
       integer(HID_T), value :: group_id
-      type(Bank), intent(out) :: bank_(*)
     end subroutine read_source_bank
   end interface
 
@@ -437,7 +435,7 @@ contains
       if (master .or. parallel) then
         file_id = file_open(filename_, 'a', parallel=.true.)
       end if
-      call write_source_bank(file_id, source_bank)
+      call write_source_bank(file_id)
       if (master .or. parallel) call file_close(file_id)
     end if
   end function openmc_statepoint_write
@@ -634,8 +632,8 @@ contains
         file_id = file_open(path_source_point, 'r', parallel=.true.)
       end if
 
-      ! Write out source
-      call read_source_bank(file_id, source_bank)
+      ! Read source
+      call read_source_bank(file_id)
 
     end if
 
diff --git a/src/state_point.cpp b/src/state_point.cpp
index fdf700c44..ce1c23836 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -8,6 +8,7 @@
 #include "xtensor/xbuilder.hpp" // for empty_like
 #include "xtensor/xview.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/constants.h"
 #include "openmc/eigenvalue.h"
@@ -70,16 +71,13 @@ write_source_point(const char* filename)
   }
 
   // Get pointer to source bank and write to file
-  Bank* source_bank;
-  int64_t n;
-  openmc_source_bank(&source_bank, &n);
-  write_source_bank(file_id, source_bank);
+  write_source_bank(file_id);
 
   if (mpi::master || parallel) file_close(file_id);
 }
 
 void
-write_source_bank(hid_t group_id, Bank* source_bank)
+write_source_bank(hid_t group_id)
 {
   hid_t banktype = h5banktype();
 
@@ -103,7 +101,7 @@ write_source_bank(hid_t group_id, Bank* source_bank)
   H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
 
   // Write data to file in parallel
-  H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank);
+  H5Dwrite(dset, banktype, memspace, dspace, plist, simulation::source_bank.data());
 
   // Free resources
   H5Sclose(dspace);
@@ -122,7 +120,8 @@ write_source_bank(hid_t group_id, Bank* source_bank)
 
     // Save source bank sites since the souce_bank array is overwritten below
 #ifdef OPENMC_MPI
-    std::vector<Bank> temp_source {source_bank, source_bank + simulation::work};
+    std::vector<Bank> temp_source {simulation::source_bank.begin(),
+      simulation::source_bank.begin() + simulation::work};
 #endif
 
     for (int i = 0; i < mpi::n_procs; ++i) {
@@ -134,7 +133,7 @@ write_source_bank(hid_t group_id, Bank* source_bank)
 #ifdef OPENMC_MPI
       // Receive source sites from other processes
       if (i > 0)
-        MPI_Recv(source_bank, count[0], mpi::bank, i, i,
+        MPI_Recv(source_bank.data(), count[0], mpi::bank, i, i,
                  mpi::intracomm, MPI_STATUS_IGNORE);
 #endif
 
@@ -144,7 +143,8 @@ write_source_bank(hid_t group_id, Bank* source_bank)
       H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
 
       // Write data to hyperslab
-      H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank);
+      H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT,
+        simulation::source_bank.data());
 
       H5Sclose(memspace);
       H5Sclose(dspace);
@@ -155,12 +155,12 @@ write_source_bank(hid_t group_id, Bank* source_bank)
 
 #ifdef OPENMC_MPI
     // Restore state of source bank
-    std::copy(temp_source.begin(), temp_source.end(), source_bank);
+    std::copy(temp_source.begin(), temp_source.end(), source_bank.begin());
 #endif
   } else {
 #ifdef OPENMC_MPI
-    MPI_Send(source_bank, simulation::work, mpi::bank, 0, mpi::rank,
-             mpi::intracomm);
+    MPI_Send(simulation::source_bank.data(), simulation::work, mpi::bank,
+      0, mpi::rank, mpi::intracomm);
 #endif
   }
 #endif
@@ -169,7 +169,7 @@ write_source_bank(hid_t group_id, Bank* source_bank)
 }
 
 
-void read_source_bank(hid_t group_id, Bank* source_bank)
+void read_source_bank(hid_t group_id)
 {
   hid_t banktype = h5banktype();
 
@@ -197,10 +197,10 @@ void read_source_bank(hid_t group_id, Bank* source_bank)
     // Read data in parallel
   hid_t plist = H5Pcreate(H5P_DATASET_XFER);
   H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
-  H5Dread(dset, banktype, memspace, dspace, plist, source_bank);
+  H5Dread(dset, banktype, memspace, dspace, plist, simulation::source_bank.data());
   H5Pclose(plist);
 #else
-  H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank);
+  H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, simulation::source_bank.data());
 #endif
 
   // Close all ids
diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index b82837ea1..1889d0dcb 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -3,6 +3,7 @@ module tally
   use, intrinsic :: ISO_C_BINDING
 
   use algorithm,        only: binary_search
+  use bank_header
   use constants
   use dict_header,      only: EMPTY
   use error,            only: fatal_error
@@ -687,7 +688,7 @@ contains
             do k = 1, p % n_bank
 
               ! get the delayed group
-              g = fission_bank(n_bank - p % n_bank + k) % delayed_group
+              g = fission_bank_delayed_group(n_bank - p % n_bank + k)
 
               ! Case for tallying delayed emissions
               if (g /= 0) then
@@ -697,8 +698,8 @@ contains
                      reactions(nuclides(p % event_nuclide) % index_fission(1)))
 
                   ! determine score based on bank site weight and keff.
-                  score = score + keff * fission_bank(n_bank - p % n_bank + k) &
-                       % wgt * rxn % product_decay_rate(1 + g) * flux
+                  score = score + keff * fission_bank_wgt(n_bank - p % n_bank + k) &
+                       * rxn % product_decay_rate(1 + g) * flux
                 end associate
 
                 ! if the delayed group filter is present, tally to corresponding
@@ -1843,7 +1844,7 @@ contains
             do k = 1, p % n_bank
 
               ! get the delayed group
-              g = fission_bank(n_bank - p % n_bank + k) % delayed_group
+              g = fission_bank_delayed_group(n_bank - p % n_bank + k)
 
               ! Case for tallying delayed emissions
               if (g /= 0) then
@@ -1851,13 +1852,13 @@ contains
                 ! determine score based on bank site weight and keff.
                 if (i_nuclide > 0) then
                   score = score + keff * atom_density * &
-                       fission_bank(n_bank - p % n_bank + k) % wgt * &
+                       fission_bank_wgt(n_bank - p % n_bank + k) * &
                        get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
                        get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
                        get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
                 else
                   score = score + keff * &
-                       fission_bank(n_bank - p % n_bank + k) % wgt * &
+                       fission_bank_wgt(n_bank - p % n_bank + k) * &
                        get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux
                 end if
 
@@ -2375,10 +2376,10 @@ contains
       do k = 1, p % n_bank
 
         ! get the delayed group
-        g = fission_bank(n_bank - p % n_bank + k) % delayed_group
+        g = fission_bank_delayed_group(n_bank - p % n_bank + k)
 
         ! determine score based on bank site weight and keff
-        score = keff * fission_bank(n_bank - p % n_bank + k) % wgt
+        score = keff * fission_bank_wgt(n_bank - p % n_bank + k)
 
         ! Add derivative information for differential tallies.  Note that the
         ! i_nuclide and atom_density arguments do not matter since this is an
@@ -2390,7 +2391,7 @@ contains
         if (.not. run_CE .and. eo_filt % matches_transport_groups) then
 
           ! determine outgoing energy group from fission bank
-          g_out = int(fission_bank(n_bank - p % n_bank + k) % E)
+          g_out = int(fission_bank_E(n_bank - p % n_bank + k))
 
           ! modify the value so that g_out = 1 corresponds to the highest
           ! energy bin
@@ -2403,10 +2404,9 @@ contains
 
           ! determine outgoing energy from fission bank
           if (run_CE) then
-            E_out = fission_bank(n_bank - p % n_bank + k) % E
+            E_out = fission_bank_E(n_bank - p % n_bank + k)
           else
-            E_out = energy_bin_avg(int(fission_bank(n_bank - p % n_bank + k) &
-                 % E))
+            E_out = energy_bin_avg(int(fission_bank_E(n_bank - p % n_bank + k)))
           end if
 
           ! If this outgoing energy falls within the energyout filter's range,