From ac29b7c719a7d07e2f3a761fca8564edf3cb3633 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 12 Jan 2018 12:36:53 -0600
Subject: [PATCH 1/8] Make cross section arrays in sum_xs contiguous

---
 src/cross_section.F90  | 24 ++++++++---------
 src/nuclide_header.F90 | 60 +++++++++++++++++++++---------------------
 2 files changed, 42 insertions(+), 42 deletions(-)

diff --git a/src/cross_section.F90 b/src/cross_section.F90
index 69e577966..f0d81f892 100644
--- a/src/cross_section.F90
+++ b/src/cross_section.F90
@@ -266,29 +266,29 @@ contains
           micro_xs(i_nuclide) % interp_factor = f
 
           ! Calculate microscopic nuclide total cross section
-          micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) &
-               + f * xs % total(i_grid + 1)
+          micro_xs(i_nuclide) % total = (ONE - f) * xs % value(SUM_XS_TOTAL,i_grid) &
+               + f * xs % value(SUM_XS_TOTAL,i_grid + 1)
 
           ! Calculate microscopic nuclide elastic cross section
-          micro_xs(i_nuclide) % elastic = (ONE - f) * xs % elastic(i_grid) &
-               + f * xs % elastic(i_grid + 1)
-
-          ! Calculate microscopic nuclide absorption cross section
-          micro_xs(i_nuclide) % absorption = (ONE - f) * xs % absorption( &
-               i_grid) + f * xs % absorption(i_grid + 1)
+          micro_xs(i_nuclide) % elastic = (ONE - f) * xs % value(SUM_XS_ELASTIC,i_grid) &
+               + f * xs % value(SUM_XS_ELASTIC,i_grid + 1)
 
           if (nuc % fissionable) then
             ! Calculate microscopic nuclide total cross section
-            micro_xs(i_nuclide) % fission = (ONE - f) * xs % fission(i_grid) &
-                 + f * xs % fission(i_grid + 1)
+            micro_xs(i_nuclide) % fission = (ONE - f) * xs % value(SUM_XS_FISSION,i_grid) &
+                 + f * xs % value(SUM_XS_FISSION,i_grid + 1)
 
             ! Calculate microscopic nuclide nu-fission cross section
-            micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % nu_fission( &
-                 i_grid) + f * xs % nu_fission(i_grid + 1)
+            micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % value(SUM_XS_NU_FISSION, &
+                 i_grid) + f * xs % value(SUM_XS_NU_FISSION,i_grid + 1)
           else
             micro_xs(i_nuclide) % fission         = ZERO
             micro_xs(i_nuclide) % nu_fission      = ZERO
           end if
+
+          ! Calculate microscopic nuclide absorption cross section
+          micro_xs(i_nuclide) % absorption = (ONE - f) * xs % value(SUM_XS_ABSORPTION, &
+               i_grid) + f * xs % value(SUM_XS_ABSORPTION,i_grid + 1)
         end associate
 
         ! Depletion-related reactions
diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90
index 38e89b648..e204f68de 100644
--- a/src/nuclide_header.F90
+++ b/src/nuclide_header.F90
@@ -36,13 +36,22 @@ module nuclide_header
     real(8), allocatable :: energy(:)     ! energy values corresponding to xs
   end type EnergyGrid
 
+  integer, parameter :: &
+       SUM_XS_TOTAL      = 1, &
+       SUM_XS_ELASTIC    = 2, &
+       SUM_XS_FISSION    = 3, &
+       SUM_XS_NU_FISSION = 4, &
+       SUM_XS_ABSORPTION = 5, &
+       SUM_XS_HEATING    = 6
+
   type SumXS
-    real(8), allocatable :: total(:)      ! total cross section
-    real(8), allocatable :: elastic(:)    ! elastic scattering
-    real(8), allocatable :: fission(:)    ! fission
-    real(8), allocatable :: nu_fission(:) ! neutron production
-    real(8), allocatable :: absorption(:) ! absorption (MT > 100)
-    real(8), allocatable :: heating(:)    ! heating
+    real(8), allocatable :: value(:,:)
+!!$    real(8), allocatable :: total(:)      ! total cross section
+!!$    real(8), allocatable :: elastic(:)    ! elastic scattering
+!!$    real(8), allocatable :: fission(:)    ! fission
+!!$    real(8), allocatable :: nu_fission(:) ! neutron production
+!!$    real(8), allocatable :: absorption(:) ! absorption (MT > 100)
+!!$    real(8), allocatable :: heating(:)    ! heating
   end type SumXS
 
   type :: Nuclide
@@ -571,16 +580,8 @@ contains
     do i = 1, n_temperature
       ! Allocate and initialize derived cross sections
       n_grid = size(this % grid(i) % energy)
-      allocate(this % sum_xs(i) % total(n_grid))
-      allocate(this % sum_xs(i) % elastic(n_grid))
-      allocate(this % sum_xs(i) % fission(n_grid))
-      allocate(this % sum_xs(i) % nu_fission(n_grid))
-      allocate(this % sum_xs(i) % absorption(n_grid))
-      this % sum_xs(i) % total(:) = ZERO
-      this % sum_xs(i) % elastic(:) = ZERO
-      this % sum_xs(i) % fission(:) = ZERO
-      this % sum_xs(i) % nu_fission(:) = ZERO
-      this % sum_xs(i) % absorption(:) = ZERO
+      allocate(this % sum_xs(i) % value(6,n_grid))
+      this % sum_xs(i) % value(:,:) = ZERO
     end do
 
     i_fission = 0
@@ -608,16 +609,17 @@ contains
           n = size(rx % xs(t) % value)
 
           ! Copy elastic
-          if (rx % MT == ELASTIC) this % sum_xs(t) % elastic(:) = rx % xs(t) % value
+          if (rx % MT == ELASTIC) this % sum_xs(t) % value(SUM_XS_ELASTIC,:) = &
+               rx % xs(t) % value
 
           ! Add contribution to total cross section
-          this % sum_xs(t) % total(j:j+n-1) = this % sum_xs(t) % total(j:j+n-1) + &
-               rx % xs(t) % value
+          this % sum_xs(t) % value(SUM_XS_TOTAL,j:j+n-1) = this % sum_xs(t) % &
+               value(SUM_XS_TOTAL,j:j+n-1) + rx % xs(t) % value
 
           ! Add contribution to absorption cross section
           if (is_disappearance(rx % MT)) then
-            this % sum_xs(t) % absorption(j:j+n-1) = this % sum_xs(t) % &
-                 absorption(j:j+n-1) + rx % xs(t) % value
+            this % sum_xs(t) % value(SUM_XS_ABSORPTION,j:j+n-1) = this % sum_xs(t) % &
+                 value(SUM_XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value
           end if
 
           ! Information about fission reactions
@@ -633,12 +635,12 @@ contains
           ! Add contribution to fission cross section
           if (is_fission(rx % MT)) then
             this % fissionable = .true.
-            this % sum_xs(t) % fission(j:j+n-1) = this % sum_xs(t) % &
-                 fission(j:j+n-1) + rx % xs(t) % value
+            this % sum_xs(t) % value(SUM_XS_FISSION,j:j+n-1) = this % sum_xs(t) % &
+                 value(SUM_XS_FISSION,j:j+n-1) + rx % xs(t) % value
 
             ! Also need to add fission cross sections to absorption
-            this % sum_xs(t) % absorption(j:j+n-1) = this % sum_xs(t) % &
-                 absorption(j:j+n-1) + rx % xs(t) % value
+            this % sum_xs(t) % value(SUM_XS_ABSORPTION,j:j+n-1) = this % sum_xs(t) % &
+                 value(SUM_XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value
 
             ! If total fission reaction is present, there's no need to store the
             ! reaction cross-section since it was copied to this % fission
@@ -689,12 +691,10 @@ contains
     ! Calculate nu-fission cross section
     do t = 1, n_temperature
       if (this % fissionable) then
-        do i = 1, size(this % sum_xs(t) % fission)
-          this % sum_xs(t) % nu_fission(i) = this % nu(this % grid(t) % energy(i), &
-               EMISSION_TOTAL) * this % sum_xs(t) % fission(i)
+        do i = 1, n_grid
+          this % sum_xs(t) % value(SUM_XS_NU_FISSION,i) = this % nu(this % grid(t) % energy(i), &
+               EMISSION_TOTAL) * this % sum_xs(t) % value(SUM_XS_FISSION,i)
         end do
-      else
-        this % sum_xs(t) % nu_fission(:) = ZERO
       end if
     end do
   end subroutine nuclide_create_derived

From d2d1703e7671285efb855af929d9ff57cb293131 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 12 Jan 2018 15:06:04 -0600
Subject: [PATCH 2/8] Aesthetic changes. Rename sum_xs, rename SUM_XS -> XS
 constants

---
 src/cross_section.F90  | 22 +++++++++---------
 src/nuclide_header.F90 | 53 +++++++++++++++++++++---------------------
 2 files changed, 37 insertions(+), 38 deletions(-)

diff --git a/src/cross_section.F90 b/src/cross_section.F90
index f0d81f892..1dfa7f52b 100644
--- a/src/cross_section.F90
+++ b/src/cross_section.F90
@@ -233,7 +233,7 @@ contains
           if (f > prn()) i_temp = i_temp + 1
         end select
 
-        associate (grid => nuc % grid(i_temp), xs => nuc % sum_xs(i_temp))
+        associate (grid => nuc % grid(i_temp), xs => nuc % xs(i_temp))
           ! Determine the energy grid index using a logarithmic mapping to
           ! reduce the energy range over which a binary search needs to be
           ! performed
@@ -266,29 +266,29 @@ contains
           micro_xs(i_nuclide) % interp_factor = f
 
           ! Calculate microscopic nuclide total cross section
-          micro_xs(i_nuclide) % total = (ONE - f) * xs % value(SUM_XS_TOTAL,i_grid) &
-               + f * xs % value(SUM_XS_TOTAL,i_grid + 1)
+          micro_xs(i_nuclide) % total = (ONE - f) * xs % value(XS_TOTAL,i_grid) &
+               + f * xs % value(XS_TOTAL,i_grid + 1)
 
           ! Calculate microscopic nuclide elastic cross section
-          micro_xs(i_nuclide) % elastic = (ONE - f) * xs % value(SUM_XS_ELASTIC,i_grid) &
-               + f * xs % value(SUM_XS_ELASTIC,i_grid + 1)
+          micro_xs(i_nuclide) % elastic = (ONE - f) * xs % value(XS_ELASTIC,i_grid) &
+               + f * xs % value(XS_ELASTIC,i_grid + 1)
 
           if (nuc % fissionable) then
             ! Calculate microscopic nuclide total cross section
-            micro_xs(i_nuclide) % fission = (ONE - f) * xs % value(SUM_XS_FISSION,i_grid) &
-                 + f * xs % value(SUM_XS_FISSION,i_grid + 1)
+            micro_xs(i_nuclide) % fission = (ONE - f) * xs % value(XS_FISSION,i_grid) &
+                 + f * xs % value(XS_FISSION,i_grid + 1)
 
             ! Calculate microscopic nuclide nu-fission cross section
-            micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % value(SUM_XS_NU_FISSION, &
-                 i_grid) + f * xs % value(SUM_XS_NU_FISSION,i_grid + 1)
+            micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % value(XS_NU_FISSION, &
+                 i_grid) + f * xs % value(XS_NU_FISSION,i_grid + 1)
           else
             micro_xs(i_nuclide) % fission         = ZERO
             micro_xs(i_nuclide) % nu_fission      = ZERO
           end if
 
           ! Calculate microscopic nuclide absorption cross section
-          micro_xs(i_nuclide) % absorption = (ONE - f) * xs % value(SUM_XS_ABSORPTION, &
-               i_grid) + f * xs % value(SUM_XS_ABSORPTION,i_grid + 1)
+          micro_xs(i_nuclide) % absorption = (ONE - f) * xs % value(XS_ABSORPTION, &
+               i_grid) + f * xs % value(XS_ABSORPTION,i_grid + 1)
         end associate
 
         ! Depletion-related reactions
diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90
index e204f68de..88e934856 100644
--- a/src/nuclide_header.F90
+++ b/src/nuclide_header.F90
@@ -36,22 +36,20 @@ module nuclide_header
     real(8), allocatable :: energy(:)     ! energy values corresponding to xs
   end type EnergyGrid
 
+  ! Positions for first dimension of Nuclide % xs
   integer, parameter :: &
-       SUM_XS_TOTAL      = 1, &
-       SUM_XS_ELASTIC    = 2, &
-       SUM_XS_FISSION    = 3, &
-       SUM_XS_NU_FISSION = 4, &
-       SUM_XS_ABSORPTION = 5, &
-       SUM_XS_HEATING    = 6
+       XS_TOTAL      = 1, &
+       XS_ELASTIC    = 2, &
+       XS_FISSION    = 3, &
+       XS_NU_FISSION = 4, &
+       XS_ABSORPTION = 5, &
+       XS_HEATING    = 6
 
+  ! The array within SumXS is of shape (6, n_energy) where the first dimension
+  ! corresponds to the following values: 1) total, 2) elastic scattering, 3)
+  ! fission, 4) neutron production, 5) absorption (MT > 100), 6) heating
   type SumXS
     real(8), allocatable :: value(:,:)
-!!$    real(8), allocatable :: total(:)      ! total cross section
-!!$    real(8), allocatable :: elastic(:)    ! elastic scattering
-!!$    real(8), allocatable :: fission(:)    ! fission
-!!$    real(8), allocatable :: nu_fission(:) ! neutron production
-!!$    real(8), allocatable :: absorption(:) ! absorption (MT > 100)
-!!$    real(8), allocatable :: heating(:)    ! heating
   end type SumXS
 
   type :: Nuclide
@@ -70,7 +68,7 @@ module nuclide_header
     type(EnergyGrid), allocatable :: grid(:)
 
     ! Microscopic cross sections
-    type(SumXS), allocatable :: sum_xs(:)
+    type(SumXS), allocatable :: xs(:)
 
     ! Resonance scattering info
     logical              :: resonant = .false. ! resonant scatterer?
@@ -575,13 +573,13 @@ contains
     type(VectorInt) :: MTs
 
     n_temperature = size(this % kTs)
-    allocate(this % sum_xs(n_temperature))
+    allocate(this % xs(n_temperature))
     this % reaction_index(:) = 0
     do i = 1, n_temperature
       ! Allocate and initialize derived cross sections
       n_grid = size(this % grid(i) % energy)
-      allocate(this % sum_xs(i) % value(6,n_grid))
-      this % sum_xs(i) % value(:,:) = ZERO
+      allocate(this % xs(i) % value(6,n_grid))
+      this % xs(i) % value(:,:) = ZERO
     end do
 
     i_fission = 0
@@ -609,17 +607,17 @@ contains
           n = size(rx % xs(t) % value)
 
           ! Copy elastic
-          if (rx % MT == ELASTIC) this % sum_xs(t) % value(SUM_XS_ELASTIC,:) = &
+          if (rx % MT == ELASTIC) this % xs(t) % value(XS_ELASTIC,:) = &
                rx % xs(t) % value
 
           ! Add contribution to total cross section
-          this % sum_xs(t) % value(SUM_XS_TOTAL,j:j+n-1) = this % sum_xs(t) % &
-               value(SUM_XS_TOTAL,j:j+n-1) + rx % xs(t) % value
+          this % xs(t) % value(XS_TOTAL,j:j+n-1) = this % xs(t) % &
+               value(XS_TOTAL,j:j+n-1) + rx % xs(t) % value
 
           ! Add contribution to absorption cross section
           if (is_disappearance(rx % MT)) then
-            this % sum_xs(t) % value(SUM_XS_ABSORPTION,j:j+n-1) = this % sum_xs(t) % &
-                 value(SUM_XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value
+            this % xs(t) % value(XS_ABSORPTION,j:j+n-1) = this % xs(t) % &
+                 value(XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value
           end if
 
           ! Information about fission reactions
@@ -635,12 +633,12 @@ contains
           ! Add contribution to fission cross section
           if (is_fission(rx % MT)) then
             this % fissionable = .true.
-            this % sum_xs(t) % value(SUM_XS_FISSION,j:j+n-1) = this % sum_xs(t) % &
-                 value(SUM_XS_FISSION,j:j+n-1) + rx % xs(t) % value
+            this % xs(t) % value(XS_FISSION,j:j+n-1) = this % xs(t) % &
+                 value(XS_FISSION,j:j+n-1) + rx % xs(t) % value
 
             ! Also need to add fission cross sections to absorption
-            this % sum_xs(t) % value(SUM_XS_ABSORPTION,j:j+n-1) = this % sum_xs(t) % &
-                 value(SUM_XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value
+            this % xs(t) % value(XS_ABSORPTION,j:j+n-1) = this % xs(t) % &
+                 value(XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value
 
             ! If total fission reaction is present, there's no need to store the
             ! reaction cross-section since it was copied to this % fission
@@ -692,8 +690,9 @@ contains
     do t = 1, n_temperature
       if (this % fissionable) then
         do i = 1, n_grid
-          this % sum_xs(t) % value(SUM_XS_NU_FISSION,i) = this % nu(this % grid(t) % energy(i), &
-               EMISSION_TOTAL) * this % sum_xs(t) % value(SUM_XS_FISSION,i)
+          this % xs(t) % value(XS_NU_FISSION,i) = &
+               this % nu(this % grid(t) % energy(i), EMISSION_TOTAL) * &
+               this % xs(t) % value(XS_FISSION,i)
         end do
       end if
     end do

From df40af793b94db0580ca42567403e1562a8bb909 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 12 Jan 2018 15:09:09 -0600
Subject: [PATCH 3/8] Minor simplification on multipole branch in
 calculate_nuclide_xs

---
 src/cross_section.F90 | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/cross_section.F90 b/src/cross_section.F90
index 1dfa7f52b..0a944f211 100644
--- a/src/cross_section.F90
+++ b/src/cross_section.F90
@@ -181,14 +181,13 @@ contains
         call multipole_eval(nuc % multipole, E, sqrtkT, sig_t, sig_a, sig_f)
 
         micro_xs(i_nuclide) % total = sig_t
-        micro_xs(i_nuclide) % absorption = sig_a
         micro_xs(i_nuclide) % elastic = sig_t - sig_a
+        micro_xs(i_nuclide) % absorption = sig_a
+        micro_xs(i_nuclide) % fission = sig_f
 
         if (nuc % fissionable) then
-          micro_xs(i_nuclide) % fission = sig_f
           micro_xs(i_nuclide) % nu_fission = sig_f * nuc % nu(E, EMISSION_TOTAL)
         else
-          micro_xs(i_nuclide) % fission    = ZERO
           micro_xs(i_nuclide) % nu_fission = ZERO
         end if
 

From 21585e937cbcae426324f497ce8492155807104a Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 12 Jan 2018 15:44:30 -0600
Subject: [PATCH 4/8] Remove elastic from sum_xs and caches

---
 src/cross_section.F90  | 29 ++++++++++++++++-------------
 src/nuclide_header.F90 | 16 +++++-----------
 src/physics.F90        | 12 ++++++++++++
 src/tallies/tally.F90  | 14 --------------
 4 files changed, 33 insertions(+), 38 deletions(-)

diff --git a/src/cross_section.F90 b/src/cross_section.F90
index 0a944f211..82e3c00e8 100644
--- a/src/cross_section.F90
+++ b/src/cross_section.F90
@@ -42,7 +42,6 @@ contains
 
     ! Set all material macroscopic cross sections to zero
     material_xs % total          = ZERO
-    material_xs % elastic        = ZERO
     material_xs % absorption     = ZERO
     material_xs % fission        = ZERO
     material_xs % nu_fission     = ZERO
@@ -115,10 +114,6 @@ contains
         material_xs % total = material_xs % total + &
              atom_density * micro_xs(i_nuclide) % total
 
-        ! Add contributions to material macroscopic scattering cross section
-        material_xs % elastic = material_xs % elastic + &
-             atom_density * micro_xs(i_nuclide) % elastic
-
         ! Add contributions to material macroscopic absorption cross section
         material_xs % absorption = material_xs % absorption + &
              atom_density * micro_xs(i_nuclide) % absorption
@@ -162,6 +157,7 @@ contains
     real(8) :: sig_t, sig_a, sig_f ! Intermediate multipole variables
 
     ! Initialize cached cross sections to zero
+    micro_xs(i_nuclide) % elastic         = ZERO
     micro_xs(i_nuclide) % thermal         = ZERO
     micro_xs(i_nuclide) % thermal_elastic = ZERO
 
@@ -181,7 +177,6 @@ contains
         call multipole_eval(nuc % multipole, E, sqrtkT, sig_t, sig_a, sig_f)
 
         micro_xs(i_nuclide) % total = sig_t
-        micro_xs(i_nuclide) % elastic = sig_t - sig_a
         micro_xs(i_nuclide) % absorption = sig_a
         micro_xs(i_nuclide) % fission = sig_f
 
@@ -268,9 +263,9 @@ contains
           micro_xs(i_nuclide) % total = (ONE - f) * xs % value(XS_TOTAL,i_grid) &
                + f * xs % value(XS_TOTAL,i_grid + 1)
 
-          ! Calculate microscopic nuclide elastic cross section
-          micro_xs(i_nuclide) % elastic = (ONE - f) * xs % value(XS_ELASTIC,i_grid) &
-               + f * xs % value(XS_ELASTIC,i_grid + 1)
+          ! Calculate microscopic nuclide absorption cross section
+          micro_xs(i_nuclide) % absorption = (ONE - f) * xs % value(XS_ABSORPTION, &
+               i_grid) + f * xs % value(XS_ABSORPTION,i_grid + 1)
 
           if (nuc % fissionable) then
             ! Calculate microscopic nuclide total cross section
@@ -284,10 +279,6 @@ contains
             micro_xs(i_nuclide) % fission         = ZERO
             micro_xs(i_nuclide) % nu_fission      = ZERO
           end if
-
-          ! Calculate microscopic nuclide absorption cross section
-          micro_xs(i_nuclide) % absorption = (ONE - f) * xs % value(XS_ABSORPTION, &
-               i_grid) + f * xs % value(XS_ABSORPTION,i_grid + 1)
         end associate
 
         ! Depletion-related reactions
@@ -450,6 +441,18 @@ contains
     micro_xs(i_nuclide) % thermal = sab_frac * (elastic + inelastic)
     micro_xs(i_nuclide) % thermal_elastic = sab_frac * elastic
 
+    ! Calculate free atom elastic cross section
+    f = micro_xs(i_nuclide) % interp_factor
+    i_grid = micro_xs(i_nuclide) % index_grid
+    i_temp = micro_xs(i_nuclide) % index_temp
+    if (i_temp > 0) then
+      associate (xs => nuclides(i_nuclide) % reactions(1) % xs(i_temp) % value)
+        micro_xs(i_nuclide) % elastic = (ONE - f)*xs(i_grid) + f*xs(i_grid + 1)
+      end associate
+    else
+      micro_xs(i_nuclide) % elastic = ZERO
+    end if
+
     ! Correct total and elastic cross sections
     micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total &
          + micro_xs(i_nuclide) % thermal &
diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90
index 88e934856..23e6e5825 100644
--- a/src/nuclide_header.F90
+++ b/src/nuclide_header.F90
@@ -39,11 +39,9 @@ module nuclide_header
   ! Positions for first dimension of Nuclide % xs
   integer, parameter :: &
        XS_TOTAL      = 1, &
-       XS_ELASTIC    = 2, &
+       XS_ABSORPTION = 2, &
        XS_FISSION    = 3, &
-       XS_NU_FISSION = 4, &
-       XS_ABSORPTION = 5, &
-       XS_HEATING    = 6
+       XS_NU_FISSION = 4
 
   ! The array within SumXS is of shape (6, n_energy) where the first dimension
   ! corresponds to the following values: 1) total, 2) elastic scattering, 3)
@@ -120,11 +118,12 @@ module nuclide_header
   type NuclideMicroXS
     ! Microscopic cross sections in barns
     real(8) :: total
-    real(8) :: elastic          ! If sab_frac is not 1 or 0, then this value is
-                                !   averaged over bound and non-bound nuclei
     real(8) :: absorption       ! absorption (disappearance)
     real(8) :: fission          ! fission
     real(8) :: nu_fission       ! neutron production from fission
+
+    real(8) :: elastic          ! If sab_frac is not 1 or 0, then this value is
+                                !   averaged over bound and non-bound nuclei
     real(8) :: thermal          ! Bound thermal elastic & inelastic scattering
     real(8) :: thermal_elastic  ! Bound thermal elastic scattering
 
@@ -155,7 +154,6 @@ module nuclide_header
 
   type MaterialMacroXS
     real(8) :: total         ! macroscopic total xs
-    real(8) :: elastic       ! macroscopic elastic scattering xs
     real(8) :: absorption    ! macroscopic absorption xs
     real(8) :: fission       ! macroscopic fission xs
     real(8) :: nu_fission    ! macroscopic production xs
@@ -606,10 +604,6 @@ contains
           j = rx % xs(t) % threshold
           n = size(rx % xs(t) % value)
 
-          ! Copy elastic
-          if (rx % MT == ELASTIC) this % xs(t) % value(XS_ELASTIC,:) = &
-               rx % xs(t) % value
-
           ! Add contribution to total cross section
           this % xs(t) % value(XS_TOTAL,j:j+n-1) = this % xs(t) % &
                value(XS_TOTAL,j:j+n-1) + rx % xs(t) % value
diff --git a/src/physics.F90 b/src/physics.F90
index d13313e05..5d31d87d7 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -338,6 +338,18 @@ contains
          micro_xs(i_nuclide) % absorption)
     sampled = .false.
 
+    ! Calculate elastic cross section if need be
+    if (micro_xs(i_nuclide) % elastic == ZERO) then
+      if (i_temp > 0) then
+        associate (xs => nuc % reactions(1) % xs(i_temp) % value)
+          micro_xs(i_nuclide) % elastic = (ONE - f)*xs(i_grid) + f*xs(i_grid + 1)
+        end associate
+      else
+        micro_xs(i_nuclide) % elastic = micro_xs(i_nuclide) % total - &
+             micro_xs(i_nuclide) % absorption
+      end if
+    end if
+
     prob = micro_xs(i_nuclide) % elastic - micro_xs(i_nuclide) % thermal
     if (prob > cutoff) then
       ! =======================================================================
diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index 7684a36aa..d3dbf888d 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -1009,20 +1009,6 @@ contains
         ! Simply count number of scoring events
         score = ONE
 
-      case (ELASTIC)
-        if (t % estimator == ESTIMATOR_ANALOG) then
-          ! Check if event MT matches
-          if (p % event_MT /= ELASTIC) cycle SCORE_LOOP
-          score = p % last_wgt * flux
-
-        else
-          if (i_nuclide > 0) then
-            score = micro_xs(i_nuclide) % elastic * atom_density * flux
-          else
-            score = material_xs % elastic * flux
-          end if
-        end if
-
       case (SCORE_FISS_Q_PROMPT)
         score = ZERO
 

From 0cf1f918f5386ffc77654b181ea5caecab713a82 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 15 Jan 2018 16:44:18 -0600
Subject: [PATCH 5/8] Put automatic arrays on stack when using gfortran (for
 multipole performance)

---
 CMakeLists.txt | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index c5150745b..51bdcc4d9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -116,9 +116,9 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
   endif()
 
   # GCC compiler options
-  list(APPEND f90flags -cpp -std=f2008ts -fbacktrace -O2)
+  list(APPEND f90flags -cpp -std=f2008ts -fbacktrace -O2 -fstack-arrays)
   if(debug)
-    list(REMOVE_ITEM f90flags -O2)
+    list(REMOVE_ITEM f90flags -O2 -fstack-arrays)
     list(APPEND f90flags -g -Wall -Wno-unused-dummy-argument -pedantic
       -fbounds-check -ffpe-trap=invalid,overflow,underflow)
     list(APPEND ldflags -g)

From 7d7453b8ffbbfa4e8c82b152c9e00a5b644e29bd Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 15 Jan 2018 17:25:21 -0600
Subject: [PATCH 6/8] Make sure elastic scattering cross section is
 precalculated for use in URR

---
 src/cross_section.F90  | 44 +++++++++++++++++++++++++++++++-----------
 src/nuclide_header.F90 |  8 ++++----
 src/physics.F90        | 15 ++++----------
 3 files changed, 41 insertions(+), 26 deletions(-)

diff --git a/src/cross_section.F90 b/src/cross_section.F90
index 82e3c00e8..22598dc42 100644
--- a/src/cross_section.F90
+++ b/src/cross_section.F90
@@ -157,7 +157,7 @@ contains
     real(8) :: sig_t, sig_a, sig_f ! Intermediate multipole variables
 
     ! Initialize cached cross sections to zero
-    micro_xs(i_nuclide) % elastic         = ZERO
+    micro_xs(i_nuclide) % elastic         = -ONE
     micro_xs(i_nuclide) % thermal         = ZERO
     micro_xs(i_nuclide) % thermal_elastic = ZERO
 
@@ -442,16 +442,7 @@ contains
     micro_xs(i_nuclide) % thermal_elastic = sab_frac * elastic
 
     ! Calculate free atom elastic cross section
-    f = micro_xs(i_nuclide) % interp_factor
-    i_grid = micro_xs(i_nuclide) % index_grid
-    i_temp = micro_xs(i_nuclide) % index_temp
-    if (i_temp > 0) then
-      associate (xs => nuclides(i_nuclide) % reactions(1) % xs(i_temp) % value)
-        micro_xs(i_nuclide) % elastic = (ONE - f)*xs(i_grid) + f*xs(i_grid + 1)
-      end associate
-    else
-      micro_xs(i_nuclide) % elastic = ZERO
-    end if
+    call calculate_elastic_xs(i_nuclide)
 
     ! Correct total and elastic cross sections
     micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total &
@@ -581,6 +572,7 @@ contains
 
       ! Multiply by smooth cross-section if needed
       if (urr % multiply_smooth) then
+        call calculate_elastic_xs(i_nuclide)
         elastic = elastic * micro_xs(i_nuclide) % elastic
         capture = capture * (micro_xs(i_nuclide) % absorption - &
              micro_xs(i_nuclide) % fission)
@@ -609,6 +601,36 @@ contains
 
   end subroutine calculate_urr_xs
 
+!===============================================================================
+! CALCULATE_ELASTIC_XS precalculates the free atom elastic scattering cross
+! section. Normally it is not needed until a collision actually occurs in a
+! material. However, in the thermal and unresolved resonance regions, we have to
+! calculate it early to adjust the total cross section correctly.
+!===============================================================================
+
+  subroutine calculate_elastic_xs(i_nuclide)
+    integer, intent(in) :: i_nuclide
+
+    integer :: i_temp
+    integer :: i_grid
+    real(8) :: f
+
+    ! Get temperature index, grid index, and interpolation factor
+    i_temp =  micro_xs(i_nuclide) % index_temp
+    i_grid =  micro_xs(i_nuclide) % index_grid
+    f      =  micro_xs(i_nuclide) % interp_factor
+
+    if (i_temp > 0) then
+      associate (xs => nuclides(i_nuclide) % reactions(1) % xs(i_temp) % value)
+        micro_xs(i_nuclide) % elastic = (ONE - f)*xs(i_grid) + f*xs(i_grid + 1)
+      end associate
+    else
+      ! For multipole, elastic is total - absorption
+      micro_xs(i_nuclide) % elastic = micro_xs(i_nuclide) % total - &
+           micro_xs(i_nuclide) % absorption
+    end if
+  end subroutine calculate_elastic_xs
+
 !===============================================================================
 ! MULTIPOLE_EVAL evaluates the windowed multipole equations for cross
 ! sections in the resolved resonance regions
diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90
index 23e6e5825..328fe123a 100644
--- a/src/nuclide_header.F90
+++ b/src/nuclide_header.F90
@@ -43,9 +43,9 @@ module nuclide_header
        XS_FISSION    = 3, &
        XS_NU_FISSION = 4
 
-  ! The array within SumXS is of shape (6, n_energy) where the first dimension
-  ! corresponds to the following values: 1) total, 2) elastic scattering, 3)
-  ! fission, 4) neutron production, 5) absorption (MT > 100), 6) heating
+  ! The array within SumXS is of shape (4, n_energy) where the first dimension
+  ! corresponds to the following values: 1) total, 2) absorption (MT > 100), 3)
+  ! fission, 4) neutron production
   type SumXS
     real(8), allocatable :: value(:,:)
   end type SumXS
@@ -576,7 +576,7 @@ contains
     do i = 1, n_temperature
       ! Allocate and initialize derived cross sections
       n_grid = size(this % grid(i) % energy)
-      allocate(this % xs(i) % value(6,n_grid))
+      allocate(this % xs(i) % value(4,n_grid))
       this % xs(i) % value(:,:) = ZERO
     end do
 
diff --git a/src/physics.F90 b/src/physics.F90
index 5d31d87d7..cb2877824 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -2,7 +2,7 @@ module physics
 
   use algorithm,              only: binary_search
   use constants
-  use cross_section,          only: elastic_xs_0K
+  use cross_section,          only: elastic_xs_0K, calculate_elastic_xs
   use endf,                   only: reaction_name
   use error,                  only: fatal_error, warning, write_message
   use material_header,        only: Material, materials
@@ -338,16 +338,9 @@ contains
          micro_xs(i_nuclide) % absorption)
     sampled = .false.
 
-    ! Calculate elastic cross section if need be
-    if (micro_xs(i_nuclide) % elastic == ZERO) then
-      if (i_temp > 0) then
-        associate (xs => nuc % reactions(1) % xs(i_temp) % value)
-          micro_xs(i_nuclide) % elastic = (ONE - f)*xs(i_grid) + f*xs(i_grid + 1)
-        end associate
-      else
-        micro_xs(i_nuclide) % elastic = micro_xs(i_nuclide) % total - &
-             micro_xs(i_nuclide) % absorption
-      end if
+    ! Calculate elastic cross section if it wasn't precalculated
+    if (micro_xs(i_nuclide) % elastic < ZERO) then
+      call calculate_elastic_xs(i_nuclide)
     end if
 
     prob = micro_xs(i_nuclide) % elastic - micro_xs(i_nuclide) % thermal

From f74c2fae563ff944902059c0d786940723a9a885 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 15 Jan 2018 17:54:34 -0600
Subject: [PATCH 7/8] Fix elastic scattering rate tallies

---
 src/tallies/tally.F90 | 33 ++++++++++++++++++++++++++++++++-
 1 file changed, 32 insertions(+), 1 deletion(-)

diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index d3dbf888d..eb17070c1 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -4,7 +4,7 @@ module tally
 
   use algorithm,        only: binary_search
   use constants
-  use cross_section,    only: multipole_deriv_eval
+  use cross_section,    only: multipole_deriv_eval, calculate_elastic_xs
   use dict_header,      only: EMPTY
   use error,            only: fatal_error
   use geometry_header
@@ -1009,6 +1009,37 @@ contains
         ! Simply count number of scoring events
         score = ONE
 
+      case (ELASTIC)
+        if (t % estimator == ESTIMATOR_ANALOG) then
+          ! Check if event MT matches
+          if (p % event_MT /= ELASTIC) cycle SCORE_LOOP
+          score = p % last_wgt * flux
+
+        else
+          if (i_nuclide > 0) then
+            if (micro_xs(i_nuclide) % elastic < ZERO) then
+              call calculate_elastic_xs(i_nuclide)
+            end if
+            score = micro_xs(i_nuclide) % elastic * atom_density * flux
+          else
+            score = ZERO
+            if (p % material /= MATERIAL_VOID) then
+              do l = 1, materials(p % material) % n_nuclides
+                ! Get atom density
+                atom_density_ = materials(p % material) % atom_density(l)
+
+                ! Get index in nuclides array
+                i_nuc = materials(p % material) % nuclide(l)
+                if (micro_xs(i_nuc) % elastic < ZERO) then
+                  call calculate_elastic_xs(i_nuc)
+                end if
+
+                score = score + micro_xs(i_nuc) % elastic * atom_density_ * flux
+              end do
+            end if
+          end if
+        end if
+
       case (SCORE_FISS_Q_PROMPT)
         score = ZERO
 

From e300c84db5618d74fa636b0a61185a16245398b1 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 17 Jan 2018 12:02:01 -0600
Subject: [PATCH 8/8] Replace -ONE with an arbitrary value in named constant
 CACHE_INVALID

---
 src/cross_section.F90  | 2 +-
 src/nuclide_header.F90 | 4 ++++
 src/physics.F90        | 2 +-
 src/tallies/tally.F90  | 4 ++--
 4 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/cross_section.F90 b/src/cross_section.F90
index 22598dc42..47e1c4fc5 100644
--- a/src/cross_section.F90
+++ b/src/cross_section.F90
@@ -157,7 +157,7 @@ contains
     real(8) :: sig_t, sig_a, sig_f ! Intermediate multipole variables
 
     ! Initialize cached cross sections to zero
-    micro_xs(i_nuclide) % elastic         = -ONE
+    micro_xs(i_nuclide) % elastic         = CACHE_INVALID
     micro_xs(i_nuclide) % thermal         = ZERO
     micro_xs(i_nuclide) % thermal_elastic = ZERO
 
diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90
index 328fe123a..1ca584914 100644
--- a/src/nuclide_header.F90
+++ b/src/nuclide_header.F90
@@ -115,6 +115,10 @@ module nuclide_header
 ! nuclide at the current energy
 !===============================================================================
 
+  ! Arbitrary value to indicate invalid cache state for elastic scattering
+  ! (NuclideMicroXS % elastic)
+  real(8), parameter :: CACHE_INVALID = dble(Z"FFE0000000000000")
+
   type NuclideMicroXS
     ! Microscopic cross sections in barns
     real(8) :: total
diff --git a/src/physics.F90 b/src/physics.F90
index cb2877824..fe242dbb7 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -339,7 +339,7 @@ contains
     sampled = .false.
 
     ! Calculate elastic cross section if it wasn't precalculated
-    if (micro_xs(i_nuclide) % elastic < ZERO) then
+    if (micro_xs(i_nuclide) % elastic == CACHE_INVALID) then
       call calculate_elastic_xs(i_nuclide)
     end if
 
diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index eb17070c1..027e335e9 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -1017,7 +1017,7 @@ contains
 
         else
           if (i_nuclide > 0) then
-            if (micro_xs(i_nuclide) % elastic < ZERO) then
+            if (micro_xs(i_nuclide) % elastic == CACHE_INVALID) then
               call calculate_elastic_xs(i_nuclide)
             end if
             score = micro_xs(i_nuclide) % elastic * atom_density * flux
@@ -1030,7 +1030,7 @@ contains
 
                 ! Get index in nuclides array
                 i_nuc = materials(p % material) % nuclide(l)
-                if (micro_xs(i_nuc) % elastic < ZERO) then
+                if (micro_xs(i_nuc) % elastic == CACHE_INVALID) then
                   call calculate_elastic_xs(i_nuc)
                 end if