diff --git a/include/openmc/output.h b/include/openmc/output.h
index 5c7ebcd59..c80ca1167 100644
--- a/include/openmc/output.h
+++ b/include/openmc/output.h
@@ -61,5 +61,9 @@ void print_generation();
 
 void print_batch_keff();
 
+//! Display time elapsed for various stages of a run
+
+void print_runtime();
+
 } // namespace openmc
 #endif // OPENMC_OUTPUT_H
diff --git a/src/output.F90 b/src/output.F90
index 3f864536f..78e784220 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -31,12 +31,6 @@ module output
   integer :: eu = ERROR_UNIT
 
   interface
-    function entropy(i) result(h) bind(C, name='entropy_c')
-      import C_INT, C_DOUBLE
-      integer(C_INT), value :: i
-      real(C_DOUBLE) :: h
-    end function
-
     subroutine print_particle(p) bind(C)
       import Particle
       type(Particle), intent(in) :: p
@@ -99,77 +93,6 @@ contains
 
   end subroutine header
 
-!===============================================================================
-! PRINT_RUNTIME displays the total time elapsed for the entire run, for
-! initialization, for computation, and for intergeneration synchronization.
-!===============================================================================
-
-  subroutine print_runtime() bind(C)
-
-    integer       :: n_active
-    real(8)       :: speed_inactive  ! # of neutrons/second in inactive batches
-    real(8)       :: speed_active    ! # of neutrons/second in active batches
-    character(15) :: string
-
-    ! display header block
-    call header("Timing Statistics", 6)
-
-    ! display time elapsed for various sections
-    write(ou,100) "Total time for initialization", time_initialize_elapsed()
-    write(ou,100) "  Reading cross sections", time_read_xs_elapsed()
-    write(ou,100) "Total time in simulation", time_inactive_elapsed() + &
-         time_active_elapsed()
-    write(ou,100) "  Time in transport only", time_transport_elapsed()
-    if (run_mode == MODE_EIGENVALUE) then
-      write(ou,100) "  Time in inactive batches", time_inactive_elapsed()
-    end if
-    write(ou,100) "  Time in active batches", time_active_elapsed()
-    if (run_mode == MODE_EIGENVALUE) then
-      write(ou,100) "  Time synchronizing fission bank", time_bank_elapsed()
-      write(ou,100) "    Sampling source sites", time_bank_sample_elapsed()
-      write(ou,100) "    SEND/RECV source sites", time_bank_sendrecv_elapsed()
-    end if
-    write(ou,100) "  Time accumulating tallies", time_tallies_elapsed()
-    write(ou,100) "Total time for finalization", time_finalize_elapsed()
-    write(ou,100) "Total time elapsed", time_total_elapsed()
-
-    ! Calculate particle rate in active/inactive batches
-    n_active = current_batch - n_inactive
-    if (restart_run) then
-      if (restart_batch < n_inactive) then
-        speed_inactive = real(n_particles * (n_inactive - restart_batch) * &
-             gen_per_batch) / time_inactive_elapsed()
-        speed_active = real(n_particles * n_active * gen_per_batch) / &
-             time_active_elapsed()
-      else
-        speed_inactive = ZERO
-        speed_active = real(n_particles * (n_batches - restart_batch) * &
-             gen_per_batch) / time_active_elapsed()
-      end if
-    else
-      if (n_inactive > 0) then
-        speed_inactive = real(n_particles * n_inactive * gen_per_batch) / &
-             time_inactive_elapsed()
-      end if
-      speed_active = real(n_particles * n_active * gen_per_batch) / &
-           time_active_elapsed()
-    end if
-
-    ! display calculation rate
-    if (.not. (restart_run .and. (restart_batch >= n_inactive)) &
-         .and. n_inactive > 0) then
-      string = to_str(speed_inactive)
-      write(ou,101) "Calculation Rate (inactive)", trim(string)
-    end if
-    string = to_str(speed_active)
-    write(ou,101) "Calculation Rate (active)", trim(string)
-
-    ! format for write statements
-100 format (1X,A,T36,"= ",ES11.4," seconds")
-101 format (1X,A,T36,"=  ",A," particles/second")
-
-  end subroutine print_runtime
-
 !===============================================================================
 ! PRINT_RESULTS displays various estimates of k-effective as well as the global
 ! leakage rate.
diff --git a/src/output.cpp b/src/output.cpp
index 0c9d35cef..39ebf0dce 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -2,7 +2,8 @@
 
 #include <algorithm>  // for std::transform
 #include <cstring>  // for strlen
-#include <iomanip>  // for setw
+#include <iomanip>  // for setw, setprecision
+#include <ios> // for fixed, scientific, left
 #include <iostream>
 #include <sstream>
 #include <ctime>
@@ -20,6 +21,7 @@
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
 #include "openmc/surface.h"
+#include "openmc/timer.h"
 
 namespace openmc {
 
@@ -341,51 +343,154 @@ void print_columns()
 
 void print_generation()
 {
+  // Save state of cout
+  auto f {std::cout.flags()};
+
   // Determine overall generation and number of active generations
   int i = overall_generation() - 1;
   int n = simulation::current_batch > settings::n_inactive ?
     settings::gen_per_batch*n_realizations + simulation::current_gen : 0;
 
+  // Set format for values
+  std::cout << std::fixed << std::setprecision(5);
+
   // write out information batch and option independent output
   std::cout << "  "  << std::setw(9) <<  std::to_string(simulation::current_batch)
-    + "/" + std::to_string(simulation::current_gen) << "   " <<
-    std::fixed << std::setw(8) << std::setprecision(5) << simulation::k_generation[i];
+    + "/" + std::to_string(simulation::current_gen) << "   " << std::setw(8)
+    << simulation::k_generation[i];
 
   // write out entropy info
-  if (settings::entropy_on) std::cout << "   " <<
-    std::fixed << std::setw(8) << std::setprecision(5) << simulation::entropy[i];
+  if (settings::entropy_on) {
+    std::cout << "   " << std::setw(8) << simulation::entropy[i];
+  }
 
   if (n > 1) {
-    std::cout << "   " <<
-      std::fixed << std::setw(8) << std::setprecision(5) << simulation::keff << " +/-" <<
-      std::fixed << std::setw(8) << std::setprecision(5) << simulation::keff_std;
+    std::cout << "   " << std::setw(8) << simulation::keff << " +/-"
+      << std::setw(8) << simulation::keff_std;
   }
   std::cout << '\n';
+
+  // Restore state of cout
+  std::cout.flags(f);
 }
 
 //==============================================================================
 
 void print_batch_keff()
 {
+  // Save state of cout
+  auto f {std::cout.flags()};
+
   // Determine overall generation and number of active generations
   int i = simulation::current_batch*settings::gen_per_batch - 1;
   int n = n_realizations*settings::gen_per_batch;
 
+  // Set format for values
+  std::cout << std::fixed << std::setprecision(5);
+
   // write out information batch and option independent output
   std::cout << "  "  << std::setw(9) <<  std::to_string(simulation::current_batch)
-    + "/" + std::to_string(settings::gen_per_batch) << "   " <<
-    std::fixed << std::setw(8) << std::setprecision(5) << simulation::k_generation[i];
+    + "/" + std::to_string(settings::gen_per_batch) << "   " << std::setw(8)
+    << simulation::k_generation[i];
 
   // write out entropy info
-  if (settings::entropy_on) std::cout << "   " <<
-    std::fixed << std::setw(8) << std::setprecision(5) << simulation::entropy[i];
+  if (settings::entropy_on) {
+    std::cout << "   " << std::setw(8) << simulation::entropy[i];
+  }
 
   if (n > 1) {
-    std::cout << "   " <<
-      std::fixed << std::setw(8) << std::setprecision(5) << simulation::keff << " +/-" <<
-      std::fixed << std::setw(8) << std::setprecision(5) << simulation::keff_std;
+    std::cout << "   " << std::setw(8) << simulation::keff << " +/-"
+      << std::setw(8) << simulation::keff_std;
   }
   std::cout << '\n';
+
+  // Restore state of cout
+  std::cout.flags(f);
+}
+
+//==============================================================================
+
+void show_time(const char* label, double value)
+{
+  std::cout << " " << std::setw(33) << std::left << label << " = "
+    << std::setw(10) << std::right << value << " seconds\n";
+}
+
+void show_rate(const char* label, double value)
+{
+  std::cout << " " << std::setw(33) << std::left << label << " = " <<
+    value << " particles/second\n";
+}
+
+void print_runtime()
+{
+  using namespace simulation;
+
+  // display header block
+  header("Timing Statistics", 6);
+
+  // Save state of cout
+  auto f {std::cout.flags()};
+
+  // display time elapsed for various sections
+  std::cout << std::scientific << std::setprecision(4);
+  show_time("Total time for initialization", time_initialize.elapsed());
+  show_time("  Reading cross sections", time_read_xs.elapsed());
+  show_time("Total time in simulation", time_inactive.elapsed() +
+    time_active.elapsed());
+  show_time("  Time in transport only", time_transport.elapsed());
+  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+    show_time("  Time in inactive batches", time_inactive.elapsed());
+  }
+  show_time("  Time in active batches", time_active.elapsed());
+  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+    show_time("  Time synchronizing fission bank", time_bank.elapsed());
+    show_time("    Sampling source sites", time_bank_sample.elapsed());
+    show_time("    SEND/RECV source sites", time_bank_sendrecv.elapsed());
+  }
+  show_time("  Time accumulating tallies", time_tallies.elapsed());
+  show_time("Total time for finalization", time_finalize.elapsed());
+  show_time("Total time elapsed", time_total.elapsed());
+
+  // Restore state of cout
+  std::cout.flags(f);
+
+  // Calculate particle rate in active/inactive batches
+  int n_active = simulation::current_batch - settings::n_inactive;
+  double speed_inactive;
+  double speed_active;
+  if (settings::restart_run) {
+    if (simulation::restart_batch < settings::n_inactive) {
+      speed_inactive = (settings::n_particles * (settings::n_inactive
+        - simulation::restart_batch) * settings::gen_per_batch)
+        / time_inactive.elapsed();
+      speed_active = (settings::n_particles * n_active
+        * settings::gen_per_batch) / time_active.elapsed();
+    } else {
+      speed_inactive = 0.0;
+      speed_active = (settings::n_particles * (settings::n_batches
+        - simulation::restart_batch) * settings::gen_per_batch)
+        / time_active.elapsed();
+    }
+  } else {
+    if (settings::n_inactive > 0) {
+      speed_inactive = (settings::n_particles * settings::n_inactive
+        * settings::gen_per_batch) / time_inactive.elapsed();
+    }
+    speed_active = (settings::n_particles * n_active * settings::gen_per_batch)
+      / time_active.elapsed();
+  }
+
+  // display calculation rate
+  std::cout << std::setprecision(6) << std::showpoint;
+  if (!(settings::restart_run && (simulation::restart_batch >= settings::n_inactive))
+      && settings::n_inactive > 0) {
+    show_rate("Calculation Rate (inactive)", speed_inactive);
+  }
+  show_rate("Calculation Rate (active)", speed_active);
+
+  // Restore state of cout
+  std::cout.flags(f);
 }
 
 } // namespace openmc
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 7f4f0dc03..3d362890c 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -32,7 +32,6 @@ extern "C" void check_triggers();
 extern "C" void init_tally_routines();
 extern "C" void load_state_point();
 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();