diff --git a/scripts/openmc-update-mgxs.py b/scripts/openmc-update-mgxs.py
index bc0c201c3e..fc9bad30f8 100755
--- a/scripts/openmc-update-mgxs.py
+++ b/scripts/openmc-update-mgxs.py
@@ -42,18 +42,20 @@ def parse_args():
     # Create argument parser
     parser = argparse.ArgumentParser(description=description,
                                      formatter_class=argparse.RawTextHelpFormatter)
-    parser.add_argument('-c', '--compression', type=int,
-                        help='input XML file')
     parser.add_argument('-i', '--input', type=argparse.FileType('r'),
                         help='input XML file')
     parser.add_argument('-o', '--output', nargs='?', default='',
                         help='output file, in HDF5 format')
+    parser.add_argument('-c', '--compression', type=int,
+                        help='HDF5 Compression Level')
     args = vars(parser.parse_args())
 
     if args['output'] == '':
-        fname = args['input'].name.split('.')
-        fname[-1] = '.h5'
-        args['output'] = ''.join(fname)
+        filename = args['input'].name
+        extension = filename[filename.rfind('.'):]
+        if extension == '.xml':
+            filename = filename[:filename.rfind('.')] + '.h5'
+        args['output'] = filename
 
     # Parse and return commandline arguments.
     return args
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 24875d1fe9..d755f21682 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -49,10 +49,18 @@ contains
 
     call read_settings_xml()
     call read_geometry_xml()
-    call read_tallies_xml()
     call read_materials()
+    call read_tallies_xml()
     if (cmfd_run) call configure_cmfd()
 
+    if (.not. run_CE) then
+      ! Create material macroscopic data for MGXS
+      call time_read_xs % start()
+      call read_mgxs()
+      call create_macro_xs()
+      call time_read_xs % stop()
+    end if
+
   end subroutine read_input_xml
 
 !===============================================================================
@@ -2092,15 +2100,9 @@ contains
     call get_temperatures(nuc_temps, sab_temps)
 
     ! Read continuous-energy cross sections
-    if (run_mode /= MODE_PLOTTING) then
+    if (run_CE .and. run_mode /= MODE_PLOTTING) then
       call time_read_xs % start()
-      if (run_CE) then
-        call read_ce_cross_sections(libraries, library_dict, nuc_temps, sab_temps)
-      else
-        ! Create material macroscopic data for MGXS
-        call read_mgxs(nuc_temps)
-        call create_macro_xs()
-      end if
+      call read_ce_cross_sections(libraries, library_dict, nuc_temps, sab_temps)
       call time_read_xs % stop()
     end if
 
diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90
index d9039a3ad7..45737ff1dc 100644
--- a/src/mgxs_data.F90
+++ b/src/mgxs_data.F90
@@ -20,9 +20,7 @@ contains
 ! nuclides and sab_tables arrays
 !===============================================================================
 
-  subroutine read_mgxs(temps)
-    type(VectorReal), allocatable :: temps(:)
-
+  subroutine read_mgxs()
     integer :: i            ! index in materials array
     integer :: j            ! index over nuclides in material
     integer :: i_xsdata     ! index in <xsdata> list
@@ -38,6 +36,7 @@ contains
     logical :: get_kfiss, get_fiss
     integer :: l
     type(DictCharInt) :: xsdata_dict
+    type(VectorReal), allocatable :: temps(:)
 
     ! Check if cross_sections.xml exists
     inquire(FILE=path_cross_sections, EXIST=file_exists)
@@ -49,6 +48,9 @@ contains
 
     call write_message("Loading Cross Section Data...", 5)
 
+    ! Get temperatures
+    call get_temperatures(temps)
+
     ! Open file for reading
     file_id = file_open(path_cross_sections, 'r', parallel=.true.)
 
@@ -234,5 +236,51 @@ contains
 
   end subroutine get_mat_kTs
 
+!===============================================================================
+! GET_TEMPERATURES returns a list of temperatures that each MGXS table
+! appears at in the model. Later, this list is used to determine the actual
+! temperatures to read (which may be different if interpolation is used)
+!===============================================================================
+
+  subroutine get_temperatures(temps)
+    type(VectorReal), allocatable, intent(out) :: temps(:)
+
+    integer :: i, j, k
+    integer :: i_nuclide    ! index in nuclides array
+    integer :: i_material
+    real(8) :: temperature  ! temperature in Kelvin
+
+    allocate(temps(n_nuclides_total))
+
+    do i = 1, size(cells)
+      do j = 1, size(cells(i) % material)
+        ! Skip any non-material cells and void materials
+        if (cells(i) % material(j) == NONE .or. &
+             cells(i) % material(j) == MATERIAL_VOID) cycle
+
+        ! Get temperature of cell (rounding to nearest integer)
+        if (size(cells(i) % sqrtkT) > 1) then
+          temperature = cells(i) % sqrtkT(j)**2 / K_BOLTZMANN
+        else
+          temperature = cells(i) % sqrtkT(1)**2 / K_BOLTZMANN
+        end if
+
+        i_material = material_dict % get_key(cells(i) % material(j))
+        associate (mat => materials(i_material))
+          NUC_NAMES_LOOP: do k = 1, size(mat % names)
+            ! Get index in temps array
+            i_nuclide = nuclide_dict % get_key(to_lower(mat % names(k)))
+
+            ! Add temperature if it hasn't already been added
+            if (find(temps(i_nuclide), temperature) == -1) then
+              call temps(i_nuclide) % push_back(temperature)
+            end if
+          end do NUC_NAMES_LOOP
+        end associate
+      end do
+    end do
+
+  end subroutine get_temperatures
+
 
 end module mgxs_data
diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90
index cef4685172..73a11476c5 100644
--- a/src/mgxs_header.F90
+++ b/src/mgxs_header.F90
@@ -691,7 +691,7 @@ module mgxs_header
               if (check_dataset(xsdata_grp, "nu-fission")) then
                 allocate(temp_arr(groups * groups * this % n_azi * this % n_pol))
                 call read_dataset(temp_arr, xsdata_grp, "nu-fission")
-                allocate(temp_4d(groups, groups, this % n_azi,this % n_pol))
+                allocate(temp_4d(groups, groups, this % n_azi, this % n_pol))
                 temp_4d(:, :, :, :) = reshape(temp_arr, (/groups, groups, &
                      this % n_azi, this % n_pol/))
                 deallocate(temp_arr)
@@ -1241,25 +1241,25 @@ module mgxs_header
       ! Create the Xs Data for each temperature
       TEMP_LOOP: do t = 1, temps % size()
         ! Allocate and initialize the data needed for macro_xs(i_mat) object
-        allocate(this % xs(t) % total(n_azi, n_pol, groups))
+        allocate(this % xs(t) % total(groups, n_azi, n_pol))
         this % xs(t) % total = ZERO
-        allocate(this % xs(t) % absorption(n_azi, n_pol, groups))
+        allocate(this % xs(t) % absorption(groups, n_azi, n_pol))
         this % xs(t) % absorption = ZERO
-        allocate(this % xs(t) % fission(n_azi, n_pol, groups))
+        allocate(this % xs(t) % fission(groups, n_azi, n_pol))
         this % xs(t) % fission = ZERO
-        allocate(this % xs(t) % k_fission(n_azi, n_pol, groups))
+        allocate(this % xs(t) % k_fission(groups, n_azi, n_pol))
         this % xs(t) % k_fission = ZERO
-        allocate(this % xs(t) % nu_fission(n_azi, n_pol, groups))
+        allocate(this % xs(t) % nu_fission(groups, n_azi, n_pol))
         this % xs(t) % nu_fission = ZERO
-        allocate(this % xs(t) % chi(n_azi, n_pol, groups, groups))
+        allocate(this % xs(t) % chi(groups, groups, n_azi, n_pol))
         this % xs(t) % chi = ZERO
-        allocate(temp_mult(n_azi, n_pol, groups, groups))
+        allocate(temp_mult(groups, groups, n_azi, n_pol))
         temp_mult = ZERO
-        allocate(mult_num(n_azi, n_pol, groups, groups))
+        allocate(mult_num(groups, groups, n_azi, n_pol))
         mult_num = ZERO
-        allocate(mult_denom(n_azi, n_pol, groups, groups))
+        allocate(mult_denom(groups, groups, n_azi, n_pol))
         mult_denom = ZERO
-        allocate(scatt_coeffs(n_azi, n_pol, order_dim, groups, groups))
+        allocate(scatt_coeffs(order_dim, groups, groups, n_azi, n_pol))
         scatt_coeffs = ZERO
 
         allocate(this % xs(t) % scatter(n_azi, n_pol))
@@ -1342,18 +1342,18 @@ module mgxs_header
                            nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % gmax(gin)
                       nuscatt = nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % scattxs(gin) * &
                            nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % energy(gin) % data(gout)
-                      mult_num(iazi, ipol, gout, gin) = &
-                           mult_num(iazi, ipol, gout, gin) + atom_density * &
+                      mult_num(gout, gin, iazi, ipol) = &
+                           mult_num(gout, gin, iazi, ipol) + atom_density * &
                            nuscatt
                       if (nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % mult(gin) % data(gout) > ZERO) then
-                        mult_denom(iazi, ipol, gout, gin) = &
-                             mult_denom(iazi, ipol, gout, gin) + atom_density * &
+                        mult_denom(gout, gin, iazi, ipol) = &
+                             mult_denom(gout, gin, iazi, ipol) + atom_density * &
                              nuscatt / &
                              nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % mult(gin) % data(gout)
                       else
                         ! Avoid division by zero
-                        mult_denom(iazi, ipol, gout, gin) = &
-                             mult_denom(iazi, ipol, gout, gin) + atom_density
+                        mult_denom(gout, gin, iazi, ipol) = &
+                             mult_denom(gout, gin, iazi, ipol) + atom_density
                       end if
                     end do
                   end do
@@ -1361,8 +1361,8 @@ module mgxs_header
                   ! Get the complete scattering matrix
                   nuc_order_dim = &
                        size(nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % dist(1) % data, dim=1)
-                  scatt_coeffs(iazi, ipol, 1:min(nuc_order_dim, order_dim), :, :) = &
-                       scatt_coeffs(iazi, ipol, 1:min(nuc_order_dim, order_dim), :, :) + &
+                  scatt_coeffs(1:min(nuc_order_dim, order_dim), :, :, iazi, ipol) = &
+                       scatt_coeffs(1:min(nuc_order_dim, order_dim), :, :, iazi, ipol) + &
                        atom_density * &
                        nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % get_matrix(min(nuc_order_dim, order_dim))
                 end do
@@ -1376,12 +1376,12 @@ module mgxs_header
               do iazi = 1, n_azi
                 do gin = 1, groups
                   do gout = 1, groups
-                    if (mult_denom(iazi, ipol, gout, gin) > ZERO) then
-                      temp_mult(iazi, ipol, gout, gin) = &
-                           mult_num(iazi, ipol, gout, gin) / &
-                           mult_denom(iazi, ipol, gout, gin)
+                    if (mult_denom(gout, gin, iazi, ipol) > ZERO) then
+                      temp_mult(gout, gin, iazi, ipol) = &
+                           mult_num(gout, gin, iazi, ipol) / &
+                           mult_denom(gout, gin, iazi, ipol)
                     else
-                      temp_mult(iazi, ipol, gout, gin) = ONE
+                      temp_mult(gout, gin, iazi, ipol) = ONE
                     end if
                   end do
                 end do
@@ -1392,8 +1392,8 @@ module mgxs_header
             do ipol = 1, n_pol
               do iazi = 1, n_azi
                 call this % xs(t) % scatter(iazi, ipol) % obj % init(&
-                     temp_mult(iazi, ipol, :, :), &
-                     scatt_coeffs(iazi, ipol, :, :, :))
+                     temp_mult(:, :, iazi, ipol), &
+                     scatt_coeffs(:, :, :, iazi, ipol))
               end do
             end do
 
@@ -1402,10 +1402,10 @@ module mgxs_header
               do ipol = 1, n_pol
                 do iazi = 1, n_azi
                   do gin = 1, groups
-                    norm =  sum(this % xs(t) % chi(iazi, ipol, :, gin))
+                    norm =  sum(this % xs(t) % chi(:, gin, iazi, ipol))
                     if (norm > ZERO) then
-                      this % xs(t) % chi(iazi, ipol, :, gin) = &
-                           this % xs(t) % chi(iazi, ipol, :, gin) / norm
+                      this % xs(t) % chi(:, gin, iazi, ipol) = &
+                           this % xs(t) % chi(:, gin, iazi, ipol) / norm
                     end if
                   end do
                 end do
diff --git a/tests/input_set.py b/tests/input_set.py
index e409fba58f..2dcef059a5 100644
--- a/tests/input_set.py
+++ b/tests/input_set.py
@@ -893,7 +893,7 @@ class MGInputSet(InputSet):
         self.settings.source = Source(space=Box([0.0, 0.0, 0.0],
                                                 [10.0, 10.0, 2.0]))
         self.settings.energy_mode = "multi-group"
-        self.settings.cross_sections = "../1d_mgxs.xml"
+        self.settings.cross_sections = "../1d_mgxs.h5"
 
     def build_defualt_plots(self):
         plot = openmc.Plot()
diff --git a/tests/test_mg_basic/inputs_true.dat b/tests/test_mg_basic/inputs_true.dat
index bfb0cc09e9..1f95579c71 100644
--- a/tests/test_mg_basic/inputs_true.dat
+++ b/tests/test_mg_basic/inputs_true.dat
@@ -1 +1 @@
-d90692c6bf8db3672d70103ded682326d7f14d416e79e18ec055bdeafe1dbc1dbe642e8410cad95bf247f27c907b10b1567127bd2868d466ce2586aacff0411f
\ No newline at end of file
+ed4d8a6c1a00d4dad10bc5522bf42d0c7a13ef5f86fc93a3f249e77681cda4424ce8879919015b05272c276ffea41e3aecef882592e8dbd97085273a356fb630
\ No newline at end of file
diff --git a/tests/test_mg_max_order/inputs_true.dat b/tests/test_mg_max_order/inputs_true.dat
index d827cfaa49..9b5d32ce6c 100644
--- a/tests/test_mg_max_order/inputs_true.dat
+++ b/tests/test_mg_max_order/inputs_true.dat
@@ -1 +1 @@
-7d508b1f3a2661566b8e8cb76fee61aecb96e8b60d633b06f13c4600bd854ea3366cdebd033c0a71ffb8adb90a9aeb64fe5ac0ef3235260921f5689c93e54305
\ No newline at end of file
+8752d3167989d876c6f29b4b1ca7eed4b91793d8579c92de85633fd90235c67d63ee76d49081006047e0078f6f86eb6e3420c1f8fe44564e7a1fb62e828587d0
\ No newline at end of file
diff --git a/tests/test_mg_nuclide/inputs_true.dat b/tests/test_mg_nuclide/inputs_true.dat
index 18811f2854..93de08656d 100644
--- a/tests/test_mg_nuclide/inputs_true.dat
+++ b/tests/test_mg_nuclide/inputs_true.dat
@@ -1 +1 @@
-be296da93031694b2915e1a10e3e6fd663d612cdd29a84745e3ccb0065c088b4bcda45f6919962f6ee784efe52ebe178bf7d3ce019859637ae57c2da1e240a04
\ No newline at end of file
+10b64fa97cc4feb1b1b8ba401b25b68ab682aba55ad4303a76a9a293ae79cdfff2f65d74310eb917da0aac243c24713d91eafae7ec5f0c4ae535d49fc928a869
\ No newline at end of file
diff --git a/tests/test_mg_tallies/inputs_true.dat b/tests/test_mg_tallies/inputs_true.dat
index c74d7759f7..5707fe6de0 100644
--- a/tests/test_mg_tallies/inputs_true.dat
+++ b/tests/test_mg_tallies/inputs_true.dat
@@ -1 +1 @@
-707285b77a091904a69720e78bf87a6779f2223c0fa2a3725353d00ca30f2624c7d53af016d78c7ad30b70a24bb53eda7dea2f592499b098f541591825b67143
\ No newline at end of file
+34af3b6ee8f65c6257d85a7c9594bc9bbabdabcb055e42bf6dd3ebfc0e41503367674223093363798950524e04f4377eae80f6d535301e55d268ff0c74e4261f
\ No newline at end of file
diff --git a/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py b/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py
index baa54b018b..a86f73796d 100644
--- a/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py
+++ b/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py
@@ -61,14 +61,14 @@ class MGXSTestHarness(PyAPITestHarness):
             self._input_set.geometry = self.mgxs_lib.create_mg_mode()
 
         # Modify settings so we can run in MG mode
-        self._input_set.settings.cross_sections = './mgxs.xml'
+        self._input_set.settings.cross_sections = './mgxs.h5'
         self._input_set.settings.energy_mode = 'multi-group'
 
         # Write modified input files
         self._input_set.settings.export_to_xml()
         self._input_set.geometry.export_to_xml()
         self._input_set.materials.export_to_xml()
-        self._input_set.mgxs_file.export_to_xml()
+        self._input_set.mgxs_file.export_to_hdf5()
         # Dont need tallies.xml, so remove the file
         if os.path.exists('./tallies.xml'):
             os.remove('./tallies.xml')
@@ -85,7 +85,7 @@ class MGXSTestHarness(PyAPITestHarness):
     def _cleanup(self):
         return
         super(MGXSTestHarness, self)._cleanup()
-        f = os.path.join(os.getcwd(), 'mgxs.xml')
+        f = os.path.join(os.getcwd(), 'mgxs.h5')
         if os.path.exists(f):
             os.remove(f)