diff --git a/CMakeLists.txt b/CMakeLists.txt
index b7517a130..99d4e4a5e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -445,7 +445,8 @@ add_library(libopenmc SHARED
   src/surface.cpp
   src/thermal.cpp
   src/xml_interface.cpp
-  src/xsdata.cpp)
+  src/xsdata.cpp
+  src/tallies/tally.cpp)
 
 set_target_properties(libopenmc PROPERTIES
   OUTPUT_NAME openmc
diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index 9c96258fb..6ead514d7 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -314,9 +314,9 @@ constexpr int MG_GET_XS_CHI_DELAYED        {14};
 // TALLY-RELATED CONSTANTS
 
 // Tally result entries
-constexpr int RESULT_VALUE  {1};
-constexpr int RESULT_SUM    {2};
-constexpr int RESULT_SUM_SQ {3};
+constexpr int RESULT_VALUE  {0};
+constexpr int RESULT_SUM    {1};
+constexpr int RESULT_SUM_SQ {2};
 
 // Tally type
 // TODO: Convert to enum
@@ -407,6 +407,7 @@ constexpr int RELATIVE_ERROR     {2};
 constexpr int STANDARD_DEVIATION {3};
 
 // Global tally parameters
+constexpr int N_GLOBAL_TALLIES {4};
 constexpr int K_COLLISION   {1};
 constexpr int K_ABSORPTION  {2};
 constexpr int K_TRACKLENGTH {3};
diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
new file mode 100644
index 000000000..a68278bb1
--- /dev/null
+++ b/include/openmc/tallies/tally.h
@@ -0,0 +1,23 @@
+#ifndef OPENMC_TALLIES_TALLY_H
+#define OPENMC_TALLIES_TALLY_H
+
+#include "openmc/constants.h"
+
+#include "xtensor/xtensor.hpp"
+
+namespace openmc {
+
+extern "C" double total_weight;
+
+xt::xtensor<double, 3> tally_results(int idx);
+
+xt::xtensor<double, 2> global_tallies();
+
+#ifdef OPENMC_MPI
+//! Collect all tally results onto master process
+extern "C" void reduce_tally_results();
+#endif
+
+} // namespace openmc
+
+#endif // OPENMC_TALLIES_TALLY_H
diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index a618f351e..ddba2dc07 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -3778,6 +3778,11 @@ contains
     real(C_DOUBLE) :: val
 
 #ifdef OPENMC_MPI
+    interface
+      subroutine reduce_tally_results() bind(C)
+      end subroutine
+    end interface
+
     ! Combine tally results onto master process
     if (reduce_tallies) call reduce_tally_results()
 #endif
@@ -3821,67 +3826,6 @@ contains
 
   end subroutine accumulate_tallies
 
-!===============================================================================
-! REDUCE_TALLY_RESULTS collects all the results from tallies onto one processor
-!===============================================================================
-
-#ifdef OPENMC_MPI
-  subroutine reduce_tally_results()
-
-    integer :: i
-    integer :: n       ! number of filter bins
-    integer :: m       ! number of score bins
-    integer :: n_bins  ! total number of bins
-    integer :: mpi_err ! MPI error code
-    real(C_DOUBLE), allocatable :: tally_temp(:,:)  ! contiguous array of results
-    real(C_DOUBLE), allocatable :: tally_temp2(:,:) ! reduced contiguous results
-    real(C_DOUBLE) :: temp(N_GLOBAL_TALLIES), temp2(N_GLOBAL_TALLIES)
-
-    do i = 1, active_tallies % size()
-      associate (t => tallies(active_tallies % data(i)) % obj)
-
-        m = size(t % results, 2)
-        n = size(t % results, 3)
-        n_bins = m*n
-
-        allocate(tally_temp(m,n), tally_temp2(m,n))
-
-        ! Reduce contiguous set of tally results
-        tally_temp = t % results(RESULT_VALUE,:,:)
-        call MPI_REDUCE(tally_temp, tally_temp2, n_bins, MPI_DOUBLE, &
-             MPI_SUM, 0, mpi_intracomm, mpi_err)
-
-        if (master) then
-          ! Transfer values to value on master
-          t % results(RESULT_VALUE,:,:) = tally_temp2
-        else
-          ! Reset value on other processors
-          t % results(RESULT_VALUE,:,:) = ZERO
-        end if
-
-        deallocate(tally_temp, tally_temp2)
-      end associate
-    end do
-
-    ! Reduce global tallies onto master
-    temp = global_tallies(RESULT_VALUE, :)
-    call MPI_REDUCE(temp, temp2, N_GLOBAL_TALLIES, MPI_DOUBLE, MPI_SUM, &
-         0, mpi_intracomm, mpi_err)
-    if (master) then
-      global_tallies(RESULT_VALUE, :) = temp2
-    else
-      global_tallies(RESULT_VALUE, :) = ZERO
-    end if
-
-    ! We also need to determine the total starting weight of particles from the
-    ! last realization
-    temp(1) = total_weight
-    call MPI_REDUCE(temp, total_weight, 1, MPI_REAL8, MPI_SUM, &
-         0, mpi_intracomm, mpi_err)
-
-  end subroutine reduce_tally_results
-#endif
-
 !===============================================================================
 ! SETUP_ACTIVE_TALLIES
 !===============================================================================
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
new file mode 100644
index 000000000..c98442b2c
--- /dev/null
+++ b/src/tallies/tally.cpp
@@ -0,0 +1,96 @@
+#include "openmc/tallies/tally.h"
+
+#include "openmc/capi.h"
+#include "openmc/constants.h"
+#include "openmc/message_passing.h"
+
+#include "xtensor/xadapt.hpp"
+#include "xtensor/xbuilder.hpp" // for empty_like
+#include "xtensor/xview.hpp"
+
+#include <cstddef>
+
+namespace openmc {
+
+xt::xtensor<double, 2> global_tallies()
+{
+  // Get pointer to global tallies
+  double* buffer;
+  openmc_global_tallies(&buffer);
+
+  // Adapt into xtensor
+  std::vector<int> shape {N_GLOBAL_TALLIES, 3};
+  int size {3*N_GLOBAL_TALLIES};
+  return xt::adapt(buffer, size, xt::no_ownership(), shape);
+}
+
+xt::xtensor<double, 3> tally_results(int idx)
+{
+  // Get pointer to tally results
+  double* results;
+  std::vector<int> shape(3);
+  openmc_tally_results(idx, &results, shape.data());
+
+  // Adapt array into xtensor with no ownership
+  int size {shape[0] * shape[1] * shape[2]};
+  return xt::adapt(results, size, xt::no_ownership(), shape);
+}
+
+#ifdef OPENMC_MPI
+void reduce_tally_results()
+{
+  for (int i = 1; i <= n_tallies; ++i) {
+    // Skip any tallies that are not active
+    bool active;
+    openmc_tally_get_active(i, &active);
+    if (!active) continue;
+
+    // Get view of accumulated tally values
+    auto results = tally_results(i);
+    auto values_view = xt::view(results, xt::all(), xt::all(), RESULT_VALUE);
+
+    // Make copy of tally values in contiguous array
+    xt::xtensor<double, 2> values = values_view;
+    xt::xtensor<double, 2> values_reduced = xt::empty_like(values);
+
+    // Reduce contiguous set of tally results
+    MPI_Reduce(values.data(), values_reduced.data(), values.size(),
+      MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm);
+
+    // Transfer values on master and reset on other ranks
+    if (mpi::master) {
+      values_view = mpi::master;
+    } else {
+      values_view = 0.0;
+    }
+  }
+
+  // Get view of global tally values
+  auto gt = global_tallies();
+  auto gt_values_view = xt::view(gt, xt::all(), RESULT_VALUE);
+
+  // Make copy of values in contiguous array
+  xt::xtensor<double, 1> gt_values = gt_values_view;
+  xt::xtensor<double, 1> gt_values_reduced = xt::empty_like(gt_values);
+
+  // Reduce contiguous data
+  MPI_Reduce(gt_values.data(), gt_values_reduced.data(), N_GLOBAL_TALLIES,
+    MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm);
+
+  // Transfer values on master and reset on other ranks
+  if (mpi::master) {
+    gt_values_view = gt_values_reduced;
+  } else {
+    gt_values_view = 0.0;
+  }
+
+  // We also need to determine the total starting weight of particles from the
+  // last realization
+  double weight_reduced;
+  MPI_Reduce(&total_weight, &weight_reduced, 1, MPI_DOUBLE, MPI_SUM,
+    0, mpi::intracomm);
+  if (mpi::master) total_weight = weight_reduced;
+}
+#endif
+
+} // namespace openmc