diff --git a/openmc/geometry.py b/openmc/geometry.py
index 079647b79b..52f3ade3b4 100644
--- a/openmc/geometry.py
+++ b/openmc/geometry.py
@@ -42,12 +42,6 @@ class Geometry(object):
     @root_universe.setter
     def root_universe(self, root_universe):
         check_type('root universe', root_universe, openmc.Universe)
-        if root_universe.id != 0:
-            msg = 'Unable to add root Universe "{0}" to Geometry since ' \
-                  'it has ID="{1}" instead of ' \
-                  'ID=0'.format(root_universe, root_universe.id)
-            raise ValueError(msg)
-
         self._root_universe = root_universe
 
     def add_volume_information(self, volume_calc):
diff --git a/src/geometry.F90 b/src/geometry.F90
index a7fa24fbf2..df963bc2ff 100644
--- a/src/geometry.F90
+++ b/src/geometry.F90
@@ -156,7 +156,7 @@ contains
     do j = 1, n_coord
       p % n_coord = j
       univ => universes(p % coord(j) % universe)
-      n = univ % n_cells
+      n = size(univ % cells)
 
       ! loop through each cell on this level
       do i = 1, n
@@ -216,7 +216,7 @@ contains
     else
       use_search_cells = .false.
       univ => universes(p % coord(j) % universe)
-      n = univ % n_cells
+      n = size(univ % cells)
     end if
 
     CELL_LOOP: do i = 1, n
@@ -1104,7 +1104,7 @@ contains
     type(Universe), pointer :: next_univ  ! next universe to cycle through
     class(Lattice), pointer :: lat        ! pointer to current lattice
 
-    n = univ % n_cells
+    n = size(univ % cells)
     offset = 0
 
     do i = 1, n
@@ -1221,7 +1221,7 @@ contains
     end if
 
     count = 0
-    n = univ % n_cells
+    n = size(univ % cells)
 
     do i = 1, n
 
@@ -1334,7 +1334,7 @@ contains
     type(Universe), pointer :: next_univ  ! next universe to loop through
     class(Lattice), pointer :: lat        ! pointer to current lattice
 
-    n = univ % n_cells
+    n = size(univ % cells)
 
     do i = 1, n
 
@@ -1423,7 +1423,7 @@ contains
     class(Lattice), pointer :: lat        ! pointer to current lattice
 
     levels_below = 0
-    do i = 1, univ % n_cells
+    do i = 1, size(univ % cells)
       c => cells(univ % cells(i))
 
       ! ====================================================================
diff --git a/src/geometry_header.F90 b/src/geometry_header.F90
index 60499fcfa6..0d3ae26401 100644
--- a/src/geometry_header.F90
+++ b/src/geometry_header.F90
@@ -17,7 +17,6 @@ module geometry_header
   type Universe
     integer :: id                     ! Unique ID
     integer :: type                   ! Type
-    integer :: n_cells                ! # of cells within
     integer, allocatable :: cells(:)  ! List of cells within
     real(8) :: x0                     ! Translation in x-coordinate
     real(8) :: y0                     ! Translation in y-coordinate
@@ -153,8 +152,8 @@ module geometry_header
     real(8), allocatable :: rotation_matrix(:,:)
   end type Cell
 
-  ! array index of universe 0
-  integer :: BASE_UNIVERSE
+  ! array index of the root universe
+  integer :: root_universe = -1
 
 contains
 
diff --git a/src/initialize.F90 b/src/initialize.F90
index a5db7d1538..1852154e84 100644
--- a/src/initialize.F90
+++ b/src/initialize.F90
@@ -16,11 +16,11 @@ module initialize
   use geometry,        only: neighbor_lists, count_instance, calc_offsets, &
                              maximum_levels
   use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
-                             &BASE_UNIVERSE
+                             root_universe
   use global
   use hdf5_interface,  only: file_open, read_attribute, file_close, &
                              hdf5_bank_t, hdf5_integer8_t
-  use input_xml,       only: read_input_xml, cells_in_univ_dict, read_plots_xml
+  use input_xml,       only: read_input_xml, read_plots_xml
   use material_header, only: Material
   use message_passing
   use mgxs_data,       only: read_mgxs, create_macro_xs
@@ -81,9 +81,6 @@ contains
     ! XML files because we need the PRNG to be initialized first
     if (run_mode == MODE_PLOTTING) call read_plots_xml()
 
-    ! Set up universe structures
-    call prepare_universes()
-
     ! Use dictionaries to redefine index pointers
     call adjust_indices()
 
@@ -97,7 +94,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(BASE_UNIVERSE)) > MAX_COORD) then
+    if (maximum_levels(universes(root_universe)) > 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.")
@@ -416,79 +413,6 @@ contains
 
   end subroutine read_command_line
 
-!===============================================================================
-! PREPARE_UNIVERSES allocates the universes array and determines the cells array
-! for each universe.
-!===============================================================================
-
-  subroutine prepare_universes()
-
-    integer              :: i                     ! index in cells array
-    integer              :: i_univ                ! index in universes array
-    integer              :: n_cells_in_univ       ! number of cells in a universe
-    integer, allocatable :: index_cell_in_univ(:) ! the index in the univ%cells
-                                                  ! array for each universe
-    type(ElemKeyValueII), pointer :: pair_list => null()
-    type(ElemKeyValueII), pointer :: current => null()
-    type(ElemKeyValueII), pointer :: next => null()
-    type(Universe),       pointer :: univ => null()
-    type(Cell),           pointer :: c => null()
-
-    allocate(universes(n_universes))
-
-    ! We also need to allocate the cell count lists for each universe. The logic
-    ! for this is a little more convoluted. In universe_dict, the (key,value)
-    ! pairs are the id of the universe and the index in the array. In
-    ! cells_in_univ_dict, it's the id of the universe and the number of cells.
-
-    pair_list => universe_dict%keys()
-    current => pair_list
-    do while (associated(current))
-      ! Find index of universe in universes array
-      i_univ = current%value
-      univ => universes(i_univ)
-      univ%id = current%key
-
-      ! Check for lowest level universe
-      if (univ%id == 0) BASE_UNIVERSE = i_univ
-
-      ! Find cell count for this universe
-      n_cells_in_univ = cells_in_univ_dict%get_key(univ%id)
-
-      ! Allocate cell list for universe
-      allocate(univ%cells(n_cells_in_univ))
-      univ%n_cells = n_cells_in_univ
-
-      ! Move to next universe
-      next => current%next
-      deallocate(current)
-      current => next
-    end do
-
-    ! Also allocate a list for keeping track of where cells have been assigned
-    ! in each universe
-
-    allocate(index_cell_in_univ(n_universes))
-    index_cell_in_univ = 0
-
-    do i = 1, n_cells
-      c => cells(i)
-
-      ! Get pointer to corresponding universe
-      i_univ = universe_dict%get_key(c%universe)
-      univ => universes(i_univ)
-
-      ! Increment the index for the cells array within the Universe object and
-      ! then store the index of the Cell object in that array
-      index_cell_in_univ(i_univ) = index_cell_in_univ(i_univ) + 1
-      univ%cells(index_cell_in_univ(i_univ)) = i
-    end do
-
-    ! Clear dictionary
-    call cells_in_univ_dict%clear()
-
-  end subroutine prepare_universes
-
 !===============================================================================
 ! ADJUST_INDICES changes the values for 'surfaces' for each cell and the
 ! material index assigned to each to the indices in the surfaces and material
@@ -794,7 +718,7 @@ contains
     if (.not. distribcell_active) return
 
     ! Count the number of instances of each cell.
-    call count_instance(universes(BASE_UNIVERSE))
+    call count_instance(universes(root_universe))
 
     ! Set the number of bins in all distribcell filters.
     do i = 1, n_tallies
@@ -881,7 +805,7 @@ contains
     ! Compute the number of unique universes containing these distribcells
     ! to determine the number of offset tables to allocate
     do i = 1, n_universes
-      do j = 1, universes(i) % n_cells
+      do j = 1, size(universes(i) % cells)
         if (cell_list % contains(universes(i) % cells(j))) then
           n_maps = n_maps + 1
         end if
@@ -904,7 +828,7 @@ contains
     ! unique distribcell array index.
     k = 1
     do i = 1, n_universes
-      do j = 1, universes(i) % n_cells
+      do j = 1, size(universes(i) % cells)
         if (cell_list % contains(universes(i) % cells(j))) then
           cells(universes(i) % cells(j)) % distribcell_index = k
           univ_list(k) = universes(i) % id
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 19af7eaf19..6d536828d0 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -10,7 +10,7 @@ module input_xml
   use energy_grid,      only: grid_method, n_log_bins
   use error,            only: fatal_error, warning
   use geometry_header,  only: Cell, Lattice, RectLattice, HexLattice, &
-                              get_temperatures
+                              get_temperatures, root_universe
   use global
   use hdf5_interface
   use list_header,      only: ListChar, ListInt, ListReal
@@ -35,9 +35,6 @@ module input_xml
   implicit none
   save
 
-  type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
-                                         ! universe contains
-
 contains
 
 !===============================================================================
@@ -1061,7 +1058,7 @@ contains
 
     integer :: i, j, k, m, i_x, i_a, input_index
     integer :: n, n_mats, n_x, n_y, n_z, n_rings, n_rlats, n_hlats
-    integer :: universe_num
+    integer :: univ_id
     integer :: n_cells_in_univ
     integer :: coeffs_reqd
     integer :: i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax
@@ -1089,6 +1086,11 @@ contains
     type(XMLNode), allocatable :: node_hlat_list(:)
     type(VectorInt) :: tokens
     type(VectorInt) :: rpn
+    type(VectorInt) :: fill_univ_ids ! List of fill universe IDs
+    type(VectorInt) :: univ_ids      ! List of all universe IDs
+    type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
+                                           ! universe contains
+
 
     ! Display output message
     call write_message("Reading geometry XML file...", 5)
@@ -1153,10 +1155,12 @@ contains
       if (check_for_node(node_cell, "universe")) then
         call get_node_value(node_cell, "universe", c % universe)
       else
-        c % universe = NONE
+        c % universe = 0
       end if
       if (check_for_node(node_cell, "fill")) then
         call get_node_value(node_cell, "fill", c % fill)
+        if (find(fill_univ_ids, c % fill) == -1) &
+             call fill_univ_ids % push_back(c % fill)
       else
         c % fill = NONE
       end if
@@ -1352,15 +1356,16 @@ contains
       ! For cells, we also need to check if there's a new universe --
       ! also for every cell add 1 to the count of cells for the
       ! specified universe
-      universe_num = c % universe
-      if (.not. cells_in_univ_dict % has_key(universe_num)) then
+      univ_id = c % universe
+      if (.not. cells_in_univ_dict % has_key(univ_id)) then
         n_universes = n_universes + 1
         n_cells_in_univ = 1
-        call universe_dict % add_key(universe_num, n_universes)
+        call universe_dict % add_key(univ_id, n_universes)
+        call univ_ids % push_back(univ_id)
       else
-        n_cells_in_univ = 1 + cells_in_univ_dict % get_key(universe_num)
+        n_cells_in_univ = 1 + cells_in_univ_dict % get_key(univ_id)
       end if
-      call cells_in_univ_dict % add_key(universe_num, n_cells_in_univ)
+      call cells_in_univ_dict % add_key(univ_id, n_cells_in_univ)
 
     end do
 
@@ -1769,7 +1774,9 @@ contains
         do k = 0, n_y - 1
           do j = 1, n_x
             lat % universes(j, n_y - k, m) = &
-                 &temp_int_array(j + n_x*k + n_x*n_y*(m-1))
+                 temp_int_array(j + n_x*k + n_x*n_y*(m-1))
+            if (find(fill_univ_ids, lat % universes(j, n_y - k, m)) == -1) &
+                 call fill_univ_ids % push_back(lat % universes(j, n_y - k, m))
           end do
         end do
       end do
@@ -1779,6 +1786,8 @@ contains
       lat % outer = NO_OUTER_UNIVERSE
       if (check_for_node(node_lat, "outer")) then
         call get_node_value(node_lat, "outer", lat % outer)
+        if (find(fill_univ_ids, lat % outer) == -1) &
+             call fill_univ_ids % push_back(lat % outer)
       end if
 
       ! Check for 'outside' nodes which are no longer supported.
@@ -1895,7 +1904,9 @@ contains
           do j = 1, k
             ! Place universe in array.
             lat % universes(i_x + n_rings, i_a + n_rings, m) = &
-                 &temp_int_array(input_index)
+                 temp_int_array(input_index)
+            if (find(fill_univ_ids, temp_int_array(input_index)) == -1) &
+                 call fill_univ_ids % push_back(temp_int_array(input_index))
             ! Walk index to closest non-adjacent right neighbor.
             i_x = i_x + 2
             i_a = i_a - 1
@@ -1920,7 +1931,9 @@ contains
           do j = 1, n_rings - mod(k-1, 2)
             ! Place universe in array.
             lat % universes(i_x + n_rings, i_a + n_rings, m) = &
-                 &temp_int_array(input_index)
+                 temp_int_array(input_index)
+            if (find(fill_univ_ids, temp_int_array(input_index)) == -1) &
+                 call fill_univ_ids % push_back(temp_int_array(input_index))
             ! Walk index to closest non-adjacent right neighbor.
             i_x = i_x + 2
             i_a = i_a - 1
@@ -1940,7 +1953,9 @@ contains
           do j = 1, n_rings - k
             ! Place universe in array.
             lat % universes(i_x + n_rings, i_a + n_rings, m) = &
-                 &temp_int_array(input_index)
+                 temp_int_array(input_index)
+            if (find(fill_univ_ids, temp_int_array(input_index)) == -1) &
+                 call fill_univ_ids % push_back(temp_int_array(input_index))
             ! Walk index to closest non-adjacent right neighbor.
             i_x = i_x + 2
             i_a = i_a - 1
@@ -1958,6 +1973,8 @@ contains
       lat % outer = NO_OUTER_UNIVERSE
       if (check_for_node(node_lat, "outer")) then
         call get_node_value(node_lat, "outer", lat % outer)
+        if (find(fill_univ_ids, lat % outer) == -1) &
+             call fill_univ_ids % push_back(lat % outer)
       end if
 
       ! Check for 'outside' nodes which are no longer supported.
@@ -1974,6 +1991,47 @@ contains
       end select
     end do HEX_LATTICES
 
+    ! ==========================================================================
+    ! SETUP UNIVERSES
+
+    ! Allocate universes, universe cell arrays, and assign base universe
+    allocate(universes(n_universes))
+    do i = 1, n_universes
+      associate (u => universes(i))
+        u % id = univ_ids % data(i)
+
+        ! Allocate cell list
+        n_cells_in_univ = cells_in_univ_dict % get_key(u % id)
+        allocate(u % cells(n_cells_in_univ))
+        u % cells(:) = 0
+
+        ! Check whether universe is a fill universe
+        if (find(fill_univ_ids, u % id) == -1) then
+          if (root_universe > 0) then
+            call fatal_error("Two or more universes are not used as fill &
+                 &universes, so it is not possible to distinguish which one &
+                 &is the root universe.")
+          else
+            root_universe = i
+          end if
+        end if
+      end associate
+    end do
+
+    do i = 1, n_cells
+      ! Get index in universes array
+      j = universe_dict%get_key(cells(i) % universe)
+
+      ! Set the first zero entry in the universe cells array to the index in the
+      ! global cells array
+      associate (u => universes(j))
+        u % cells(find(u % cells, 0)) = i
+      end associate
+    end do
+
+    ! Clear dictionary
+    call cells_in_univ_dict%clear()
+
     ! Close geometry XML file
     call doc % clear()
 
diff --git a/src/output.F90 b/src/output.F90
index 600d9a0f85..9003225ac3 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -6,7 +6,7 @@ module output
   use endf,            only: reaction_name
   use error,           only: fatal_error, warning
   use geometry_header, only: Cell, Universe, Lattice, RectLattice, &
-                             HexLattice, BASE_UNIVERSE
+                             HexLattice
   use global
   use math,            only: t_percentile
   use mesh_header,     only: RegularMesh
diff --git a/src/particle_header.F90 b/src/particle_header.F90
index 32ad6adfc4..a7309393b2 100644
--- a/src/particle_header.F90
+++ b/src/particle_header.F90
@@ -4,7 +4,7 @@ module particle_header
   use constants,       only: NEUTRON, ONE, NONE, ZERO, MAX_SECONDARY, &
                              MAX_DELAYED_GROUPS, ERROR_REAL
   use error,           only: fatal_error
-  use geometry_header, only: BASE_UNIVERSE
+  use geometry_header, only: root_universe
 
   implicit none
 
@@ -144,7 +144,7 @@ contains
     this % g = 1
 
     ! Set up base level coordinates
-    this % coord(1) % universe = BASE_UNIVERSE
+    this % coord(1) % universe = root_universe
     this % n_coord = 1
 
   end subroutine initialize_particle
diff --git a/src/particle_restart.F90 b/src/particle_restart.F90
index b4ce841fe2..262cef8021 100644
--- a/src/particle_restart.F90
+++ b/src/particle_restart.F90
@@ -4,7 +4,6 @@ module particle_restart
 
   use bank_header,      only: Bank
   use constants
-  use geometry_header,  only: BASE_UNIVERSE
   use global
   use hdf5_interface,   only: file_open, file_close, read_dataset
   use output,           only: write_message, print_particle
diff --git a/src/plot.F90 b/src/plot.F90
index ee8ffdb14f..6a3a191703 100644
--- a/src/plot.F90
+++ b/src/plot.F90
@@ -3,7 +3,7 @@ module plot
   use constants
   use error,           only: fatal_error
   use geometry,        only: find_cell, check_cell_overlap
-  use geometry_header, only: Cell, BASE_UNIVERSE
+  use geometry_header, only: Cell, root_universe
   use global
   use hdf5_interface
   use mesh,            only: get_mesh_indices
@@ -163,7 +163,7 @@ contains
     call p % initialize()
     p % coord(1) % xyz = xyz
     p % coord(1) % uvw = [ HALF, HALF, HALF ]
-    p % coord(1) % universe = BASE_UNIVERSE
+    p % coord(1) % universe = root_universe
 
 !$omp parallel do firstprivate(p) private(x, rgb, id) reduction(+ : data)
     do y = 1, height
@@ -367,7 +367,7 @@ contains
     call p % initialize()
     p % coord(1) % xyz = ll
     p % coord(1) % uvw = [ HALF, HALF, HALF ]
-    p % coord(1) % universe = BASE_UNIVERSE
+    p % coord(1) % universe = root_universe
 
     ! Open binary plot file for writing
     file_id = file_create(pl%path_plot)
diff --git a/src/source.F90 b/src/source.F90
index f77fe1c050..e4a01ce049 100644
--- a/src/source.F90
+++ b/src/source.F90
@@ -12,7 +12,6 @@ module source
   use distribution_multivariate, only: SpatialBox
   use error,            only: fatal_error
   use geometry,         only: find_cell
-  use geometry_header,  only: BASE_UNIVERSE
   use global
   use hdf5_interface,   only: file_create, file_open, file_close, read_dataset
   use message_passing,  only: rank
diff --git a/src/summary.F90 b/src/summary.F90
index 354f1db551..2745c4576e 100644
--- a/src/summary.F90
+++ b/src/summary.F90
@@ -4,7 +4,7 @@ module summary
 
   use constants
   use endf,            only: reaction_name
-  use geometry_header, only: BASE_UNIVERSE, Cell, Universe, Lattice, &
+  use geometry_header, only: root_universe, Cell, Universe, Lattice, &
                              RectLattice, HexLattice
   use global
   use hdf5_interface
@@ -127,7 +127,6 @@ contains
     character(MAX_LINE_LEN)              :: path
     type(Cell),     pointer :: c
     class(Surface), pointer :: s
-    type(Universe), pointer :: u
     class(Lattice), pointer :: lat
 
     ! Use H5LT interface to write number of geometry objects
@@ -236,7 +235,7 @@ contains
         do k = 1, c % instances
           path = ''
           offset = 1
-          call find_offset(i, universes(BASE_UNIVERSE), k, offset, path)
+          call find_offset(i, universes(root_universe), k, offset, path)
           paths(k) = path
         end do
         call write_dataset(cell_group, "paths", paths)
@@ -354,21 +353,22 @@ contains
 
     ! Write information on each universe
     UNIVERSE_LOOP: do i = 1, n_universes
-      u => universes(i)
-      univ_group = create_group(universes_group, "universe " // &
-           trim(to_str(u%id)))
+      associate (u => universes(i))
+        univ_group = create_group(universes_group, "universe " // &
+             trim(to_str(u%id)))
 
-      ! Write list of cells in this universe
-      if (u % n_cells > 0) then
-        allocate(cell_ids(u % n_cells))
-        do j = 1, u % n_cells
-          cell_ids(j) = cells(u % cells(j)) % id
-        end do
-        call write_dataset(univ_group, "cells", cell_ids)
-        deallocate(cell_ids)
-      end if
+        ! Write list of cells in this universe
+        if (size(u % cells) > 0) then
+          allocate(cell_ids(size(u % cells)))
+          do j = 1, size(u % cells)
+            cell_ids(j) = cells(u % cells(j)) % id
+          end do
+          call write_dataset(univ_group, "cells", cell_ids)
+          deallocate(cell_ids)
+        end if
 
-      call close_group(univ_group)
+        call close_group(univ_group)
+      end associate
     end do UNIVERSE_LOOP
 
     call close_group(universes_group)
diff --git a/src/tally_filter.F90 b/src/tally_filter.F90
index 6eadc42429..f88362be33 100644
--- a/src/tally_filter.F90
+++ b/src/tally_filter.F90
@@ -3,7 +3,7 @@ module tally_filter
   use algorithm,           only: binary_search
   use constants,           only: ONE, NO_BIN_FOUND, FP_PRECISION, ERROR_REAL
   use dict_header,         only: DictIntInt
-  use geometry_header,     only: BASE_UNIVERSE, RectLattice, HexLattice
+  use geometry_header,     only: root_universe, RectLattice, HexLattice
   use global
   use hdf5_interface
   use mesh_header,         only: RegularMesh
@@ -791,7 +791,7 @@ contains
     integer                 :: offset
     type(Universe), pointer :: univ
 
-    univ => universes(BASE_UNIVERSE)
+    univ => universes(root_universe)
     offset = 0
     label = ''
     call find_offset(this % cell, univ, bin-1, offset, label)
@@ -1387,6 +1387,7 @@ contains
     integer :: cell_index           ! Index in cells array
     integer :: lat_offset           ! Offset from lattice
     integer :: temp_offset          ! Looped sum of offsets
+    integer :: i_univ               ! index in universes array
     logical :: this_cell = .false.  ! Advance in this cell?
     logical :: later_cell = .false. ! Fill cells after this one?
     type(Cell), pointer :: c           ! Pointer to current cell
@@ -1396,10 +1397,11 @@ contains
     ! Get the distribcell index for this cell
     map = cells(goal) % distribcell_index
 
-    n = univ % n_cells
+    n = size(univ % cells)
 
     ! Write to the geometry stack
-    if (univ%id == 0) then
+    i_univ = universe_dict % get_key(univ % id)
+    if (i_univ == root_universe) then
       path = trim(path) // to_str(univ%id)
     else
       path = trim(path) // "->" // to_str(univ%id)
diff --git a/src/tracking.F90 b/src/tracking.F90
index a673b67516..7650f598ab 100644
--- a/src/tracking.F90
+++ b/src/tracking.F90
@@ -5,7 +5,6 @@ module tracking
   use error,              only: fatal_error, warning
   use geometry,           only: find_cell, distance_to_boundary, cross_surface, &
                                 cross_lattice, check_cell_overlap
-  use geometry_header,    only: Universe, BASE_UNIVERSE
   use global
   use output,             only: write_message
   use message_passing