diff --git a/src/geometry.F90 b/src/geometry.F90
index c3b647a822..7fdd03ad8b 100644
--- a/src/geometry.F90
+++ b/src/geometry.F90
@@ -534,73 +534,4 @@ contains
 
   end subroutine neighbor_lists
 
-!===============================================================================
-! MAXIMUM_LEVELS determines the maximum number of nested coordinate levels in
-! the geometry
-!===============================================================================
-
-  recursive function maximum_levels(univ) result(levels)
-
-    type(Universe), intent(in) :: univ  ! universe to search through
-    integer :: levels                   ! maximum number of levels for this universe
-
-    integer :: i                          ! index over cells
-    integer :: nx, ny, nz                 ! lattice shape
-    integer :: j, k, m                    ! indices in lattice
-    integer :: levels_below               ! max levels below this universe
-    type(Cell),     pointer :: c          ! pointer to current cell
-    type(Universe), pointer :: next_univ  ! next universe to loop through
-    class(Lattice), pointer :: lat        ! pointer to current lattice
-
-    levels_below = 0
-    do i = 1, size(univ % cells)
-      c => cells(univ % cells(i))
-
-      ! ====================================================================
-      ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
-      if (c % type() == FILL_UNIVERSE) then
-
-        next_univ => universes(c % fill)
-        levels_below = max(levels_below, maximum_levels(next_univ))
-
-      ! ====================================================================
-      ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
-      elseif (c % type() == FILL_LATTICE) then
-
-        ! Set current lattice
-        lat => lattices(c % fill) % obj
-
-        select type (lat)
-
-        type is (RectLattice)
-          nx = lat % n_cells(1)
-          ny = lat % n_cells(2)
-          nz = lat % n_cells(3)
-
-        type is (HexLattice)
-          nx = 2 * lat % n_rings - 1
-          ny = 2 * lat % n_rings - 1
-          nz = lat % n_axial
-
-        end select
-
-        ! Loop over lattice coordinates
-        do j = 1, nx
-          do k = 1, ny
-            do m = 1, nz
-              if (lat % are_valid_indices([j, k, m])) then
-                next_univ => universes(lat % get([j-1, k-1, m-1]) + 1)
-                levels_below = max(levels_below, maximum_levels(next_univ))
-              end if
-            end do
-          end do
-        end do
-
-      end if
-    end do
-
-    levels = 1 + levels_below
-
-  end function maximum_levels
-
 end module geometry
diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp
index 0685542409..6efe1cf99e 100644
--- a/src/geometry_aux.cpp
+++ b/src/geometry_aux.cpp
@@ -1,5 +1,6 @@
 #include "geometry_aux.h"
 
+#include <algorithm>  // for std::max
 #include <sstream>
 #include <unordered_set>
 
@@ -197,4 +198,29 @@ fill_offset_tables(int32_t target_univ_id, int map)
   }
 }
 
+//==============================================================================
+
+int
+maximum_levels(int32_t univ)
+{
+  int levels_below {0};
+
+  for (int32_t cell_indx : universes_c[univ]->cells) {
+    Cell &c = *cells_c[cell_indx];
+    if (c.type == FILL_UNIVERSE) {
+      int32_t next_univ = c.fill - 1; // TODO: off-by-one
+      levels_below = std::max(levels_below, maximum_levels(next_univ));
+    } else if (c.type == FILL_LATTICE) {
+      Lattice &lat = *lattices_c[c.fill - 1]; //TODO: off-by-one
+      for (auto it = lat.begin(); it != lat.end(); ++it) {
+        int32_t next_univ = *it;
+        levels_below = std::max(levels_below, maximum_levels(next_univ));
+      }
+    }
+  }
+
+  ++levels_below;
+  return levels_below;
+}
+
 } // namespace openmc
diff --git a/src/geometry_aux.h b/src/geometry_aux.h
index 0a0500794b..877f51e4e0 100644
--- a/src/geometry_aux.h
+++ b/src/geometry_aux.h
@@ -62,5 +62,14 @@ count_universe_instances(int32_t search_univ, int32_t target_univ_id);
 
 extern "C" void fill_offset_tables(int32_t target_univ_id, int map);
 
+//==============================================================================
+//! Determine the maximum number of nested coordinate levels in the geometry.
+//! @param univ The index of the universe to begin seraching from (probably the
+//!   root universe).
+//! @return The number of coordinate levels.
+//==============================================================================
+
+extern "C" int maximum_levels(int32_t univ);
+
 } // namespace openmc
 #endif // GEOMETRY_AUX_H
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index b708974b86..d325928a66 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -11,7 +11,7 @@ module input_xml
   use distribution_univariate
   use endf,             only: reaction_name
   use error,            only: fatal_error, warning, write_message, openmc_err_msg
-  use geometry,         only: maximum_levels, neighbor_lists
+  use geometry,         only: neighbor_lists
   use geometry_header
   use hdf5_interface
   use list_header,      only: ListChar, ListInt, ListReal
@@ -85,6 +85,12 @@ module input_xml
       import C_INT32_T
       integer(C_INT32_T) :: root
     end function find_root_universe
+
+    function maximum_levels(univ) bind(C) result(n)
+      import C_INT32_T, C_INT
+      integer(C_INT32_T), intent(in), value :: univ
+      integer(C_INT)                        :: n
+    end function maximum_levels
   end interface
 
 contains
@@ -169,7 +175,7 @@ contains
     ! Check to make sure there are not too many nested coordinate levels in the
     ! geometry since the coordinate list is statically allocated for performance
     ! reasons
-    if (maximum_levels(universes(root_universe)) > MAX_COORD) then
+    if (maximum_levels(root_universe - 1) > MAX_COORD) then
       call fatal_error("Too many nested coordinate levels in the geometry. &
            &Try increasing the maximum number of coordinate levels by &
            &providing the CMake -Dmaxcoord= option.")