diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h
index f9f7a9488..a6279f285 100644
--- a/include/openmc/simulation.h
+++ b/include/openmc/simulation.h
@@ -18,6 +18,7 @@ namespace simulation {
 extern "C" int current_batch;    //!< current batch
 extern "C" int current_gen;      //!< current fission generation
 extern "C" int64_t current_work; //!< index in source back of current particle
+extern "C" bool initialized;     //!< has simulation been initialized?
 extern "C" double keff;          //!< average k over batches
 extern "C" double keff_std;      //!< standard deviation of average k
 extern "C" double k_col_abs;     //!< sum over batches of k_collision * k_absorption
@@ -28,7 +29,6 @@ extern "C" int n_lost_particles; //!< cumulative number of lost particles
 extern "C" bool need_depletion_rx; //!< need to calculate depletion rx?
 extern "C" int restart_batch;   //!< batch at which a restart job resumed
 extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
-extern "C" bool simulation_initialized; //!< has simulation been initialized?
 extern "C" int total_gen;        //!< total number of generations simulated
 extern "C" int64_t work;         //!< number of particles per process
 
diff --git a/src/api.F90 b/src/api.F90
index a85a769b2..869660ae8 100644
--- a/src/api.F90
+++ b/src/api.F90
@@ -78,7 +78,6 @@ module openmc_api
   public :: openmc_plot_geometry
   public :: openmc_reset
   public :: openmc_set_seed
-  public :: openmc_simulation_init
   public :: openmc_source_bank
   public :: openmc_tally_allocate
   public :: openmc_tally_get_estimator
diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90
index 400522f74..418c2f2e1 100644
--- a/src/cmfd_execute.F90
+++ b/src/cmfd_execute.F90
@@ -375,7 +375,7 @@ contains
 ! CMFD_TALLY_INIT
 !===============================================================================
 
-  subroutine cmfd_tally_init()
+  subroutine cmfd_tally_init() bind(C)
     integer :: i
     if (cmfd_run) then
       do i = 1, size(cmfd_tallies)
diff --git a/src/output.F90 b/src/output.F90
index 58965456b..3b3bd3c4d 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -285,7 +285,7 @@ contains
 ! below them
 !===============================================================================
 
-  subroutine print_columns()
+  subroutine print_columns() bind(C)
 
     write(UNIT=ou, FMT='(2X,A9,3X)', ADVANCE='NO') "Bat./Gen."
     write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "   k    "
diff --git a/src/simulation.F90 b/src/simulation.F90
index d8a6e2d98..2bcdf68ab 100644
--- a/src/simulation.F90
+++ b/src/simulation.F90
@@ -20,8 +20,7 @@ module simulation
   use message_passing
   use mgxs_interface,  only: energy_bins, energy_bin_avg
   use nuclide_header,  only: micro_xs, n_nuclides
-  use output,          only: header, print_columns, &
-                             print_batch_keff, print_generation
+  use output,          only: print_batch_keff
   use particle_header
   use photon_header,   only: micro_photon_xs, n_elements
   use random_lcg,      only: set_particle_seed
@@ -41,16 +40,12 @@ module simulation
   implicit none
   private
   public :: openmc_next_batch
-  public :: openmc_simulation_init
 
   integer(C_INT), parameter :: STATUS_EXIT_NORMAL = 0
   integer(C_INT), parameter :: STATUS_EXIT_MAX_BATCH = 1
   integer(C_INT), parameter :: STATUS_EXIT_ON_TRIGGER = 2
 
   interface
-    subroutine openmc_simulation_init_c() bind(C)
-    end subroutine
-
     subroutine initialize_source() bind(C)
     end subroutine
 
@@ -256,32 +251,10 @@ contains
 ! INITIALIZE_SIMULATION
 !===============================================================================
 
-  function openmc_simulation_init() result(err) bind(C)
-    integer(C_INT) :: err
+  subroutine simulation_init_f() bind(C)
 
     integer :: i
 
-    err = 0
-
-    ! Skip if simulation has already been initialized
-    if (simulation_initialized) return
-
-    ! Call initialization on C++ side
-    call openmc_simulation_init_c()
-
-    ! Set up tally procedure pointers
-    call init_tally_routines()
-
-    ! Allocate source bank, and for eigenvalue simulations also allocate the
-    ! fission bank
-    call allocate_banks()
-
-    ! Allocate tally results arrays if they're not allocated yet
-    call configure_tallies()
-
-    ! Activate the CMFD tallies
-    call cmfd_tally_init()
-
     ! Set up material nuclide index mapping
     do i = 1, n_materials
       call materials(i) % init_nuclide_index()
@@ -300,36 +273,7 @@ contains
     end do
 !$omp end parallel
 
-    ! Reset global variables -- this is done before loading state point (as that
-    ! will potentially populate k_generation and entropy)
-    current_batch = 0
-    call k_generation_clear()
-    call entropy_clear()
-    need_depletion_rx = .false.
-
-    ! If this is a restart run, load the state point data and binary source
-    ! file
-    if (restart_run) then
-      call load_state_point()
-      call write_message("Resuming simulation...", 6)
-    else
-      call initialize_source()
-    end if
-
-    ! Display header
-    if (master) then
-      if (run_mode == MODE_FIXEDSOURCE) then
-        call header("FIXED SOURCE TRANSPORT SIMULATION", 3)
-      elseif (run_mode == MODE_EIGENVALUE) then
-        call header("K EIGENVALUE SIMULATION", 3)
-        if (verbosity >= 7) call print_columns()
-      end if
-    end if
-
-    ! Set flag indicating initialization is done
-    simulation_initialized = .true.
-
-  end function openmc_simulation_init
+  end subroutine
 
 !===============================================================================
 ! FINALIZE_SIMULATION calculates tally statistics, writes tallies, and displays
@@ -358,7 +302,7 @@ contains
 ! ALLOCATE_BANKS allocates memory for the fission and source banks
 !===============================================================================
 
-  subroutine allocate_banks()
+  subroutine allocate_banks() bind(C)
 
     integer :: alloc_err  ! allocation error code
 
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 39cb31361..60c573740 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -14,15 +14,26 @@
 #include <string>
 
 // data/functions from Fortran side
+extern "C" void allocate_banks();
 extern "C" void cmfd_init_batch();
+extern "C" void cmfd_tally_init();
+extern "C" void configure_tallies();
+extern "C" void init_tally_routines();
 extern "C" void join_bank_from_threads();
+extern "C" void load_state_point();
+extern "C" void print_columns();
 extern "C" void print_generation();
 extern "C" void print_results();
 extern "C" void print_runtime();
 extern "C" void setup_active_tallies();
+extern "C" void simulation_init_f();
 extern "C" void simulation_finalize_f();
 extern "C" void write_tallies();
 
+//==============================================================================
+// C API functions
+//==============================================================================
+
 // OPENMC_RUN encompasses all the main logic where iterations are performed
 // over the batches, generations, and histories in a fixed source or k-eigenvalue
 // calculation.
@@ -40,12 +51,69 @@ int openmc_run() {
   return err;
 }
 
+int openmc_simulation_init()
+{
+  using namespace openmc;
+
+  // Skip if simulation has already been initialized
+  if (simulation::initialized) return 0;
+
+  // Determine how much work each process should do
+  calculate_work();
+
+  // Set up tally procedure pointers
+  init_tally_routines();
+
+  // Allocate source bank, and for eigenvalue simulations also allocate the
+  // fission bank
+  allocate_banks();
+
+  // Allocate tally results arrays if they're not allocated yet
+  configure_tallies();
+
+  // Activate the CMFD tallies
+  cmfd_tally_init();
+
+  // Call Fortran initialization
+  simulation_init_f();
+
+  // Reset global variables -- this is done before loading state point (as that
+  // will potentially populate k_generation and entropy)
+  simulation::current_batch = 0;
+  simulation::k_generation.clear();
+  entropy.clear();
+  simulation::need_depletion_rx = false;
+
+  // If this is a restart run, load the state point data and binary source
+  // file
+  if (settings::restart_run) {
+    load_state_point();
+    write_message("Resuming simulation...", 6);
+  } else {
+    initialize_source();
+  }
+
+  // Display header
+  if (mpi::master) {
+    if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
+      header("FIXED SOURCE TRANSPORT SIMULATION", 3);
+    } else if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+      header("K EIGENVALUE SIMULATION", 3);
+      if (settings::verbosity >= 7) print_columns();
+    }
+  }
+
+  // Set flag indicating initialization is done
+  simulation::initialized = true;
+  return 0;
+}
+
 int openmc_simulation_finalize()
 {
   using namespace openmc;
 
   // Skip if simulation was never run
-  if (!simulation::simulation_initialized) return 0;
+  if (!simulation::initialized) return 0;
 
   // Stop active batch timer and start finalization timer
   time_active.stop();
@@ -80,7 +148,7 @@ int openmc_simulation_finalize()
 
   // Reset flags
   simulation::need_depletion_rx = false;
-  simulation::simulation_initialized = false;
+  simulation::initialized = false;
   return 0;
 }
 
@@ -95,6 +163,7 @@ namespace simulation {
 int current_batch;
 int current_gen;
 int64_t current_work;
+bool initialized {false};
 double keff {1.0};
 double keff_std;
 double k_col_abs {0.0};
@@ -105,7 +174,6 @@ int n_lost_particles {0};
 bool need_depletion_rx {false};
 int restart_batch;
 bool satisfy_triggers {false};
-bool simulation_initialized {false};
 int total_gen {0};
 int64_t work;
 
@@ -125,12 +193,6 @@ int thread_id;  //!< ID of a given thread
 // Non-member functions
 //==============================================================================
 
-void openmc_simulation_init_c()
-{
-  // Determine how much work each process should do
-  calculate_work();
-}
-
 void initialize_batch()
 {
   // Increment current batch
diff --git a/src/simulation_header.F90 b/src/simulation_header.F90
index 1844c8da6..288a8422c 100644
--- a/src/simulation_header.F90
+++ b/src/simulation_header.F90
@@ -23,7 +23,7 @@ module simulation_header
   integer(C_INT), bind(C) :: current_batch     ! current batch
   integer(C_INT), bind(C) :: current_gen       ! current generation within a batch
   integer(C_INT), bind(C) :: total_gen         ! total number of generations simulated
-  logical(C_BOOL), bind(C) :: simulation_initialized
+  logical(C_BOOL), bind(C, name='initialized') :: simulation_initialized
   logical(C_BOOL), bind(C) :: need_depletion_rx ! need to calculate depletion reaction rx?
 
   ! ============================================================================
diff --git a/src/state_point.F90 b/src/state_point.F90
index 97e2a8e8a..63a4d1dd4 100644
--- a/src/state_point.F90
+++ b/src/state_point.F90
@@ -485,7 +485,7 @@ contains
 ! LOAD_STATE_POINT
 !===============================================================================
 
-  subroutine load_state_point()
+  subroutine load_state_point() bind(C)
 
     integer :: i
     integer :: int_array(3)
diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index a6aa54265..8dfcf5934 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -52,7 +52,7 @@ contains
 ! with the CE and MG modes.
 !===============================================================================
 
-  subroutine init_tally_routines()
+  subroutine init_tally_routines() bind(C)
     if (run_CE) then
       score_general      => score_general_ce
       score_analog_tally => score_analog_tally_ce
diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90
index 8a5f9b550..d468ac7b7 100644
--- a/src/tallies/tally_header.F90
+++ b/src/tallies/tally_header.F90
@@ -392,7 +392,7 @@ contains
 ! tallies.xml file.
 !===============================================================================
 
-  subroutine configure_tallies()
+  subroutine configure_tallies() bind(C)
 
     integer :: i