Inelastic scattering index list

src/nuclide_header.F90

@@ -100,6 +100,7 @@ module nuclide_header
 
     ! Reactions
     type(Reaction), allocatable :: reactions(:)
+    integer, allocatable :: index_inelastic_scatter(:)
 
     ! Array that maps MT values to index in reactions; used at tally-time. Note
     ! that ENDF-102 does not assign any MT values above 891.
@@ -302,6 +303,7 @@ contains
     real(8) :: temp_actual
     type(VectorInt) :: MTs
     type(VectorInt) :: temps_to_read
+    type(VectorInt) :: index_inelastic_scatter_vector
 
     ! Get name of nuclide from group
     name_len = len(this % name)
@@ -469,10 +471,26 @@ contains
         end if
       end if
 
+      ! Add the reaction index to the scattering array if this is an inelastic
+      ! scatter reaction
+      if (MTs % data(i) /= N_FISSION .and. MTs % data(i) /= N_F .and. &
+          MTs % data(i) /= N_NF .and. MTs % data(i) /= N_2NF .and. &
+          MTs % data(i) /= N_3NF .and. MTs % data(i) < 200 .and. &
+          MTs % data(i) /= N_LEVEL .and. MTs % data(i) /= ELASTIC) then
+
+        call index_inelastic_scatter_vector % push_back(i)
+      end if
+
       call close_group(rx_group)
     end do
     call close_group(rxs_group)
 
+    ! Recast to a regular array to save space
+    allocate(this % index_inelastic_scatter( &
+         index_inelastic_scatter_vector % size()))
+    this % index_inelastic_scatter = index_inelastic_scatter_vector % data(:)
+    call index_inelastic_scatter_vector % clear()
+
     ! Read unresolved resonance probability tables if present
     if (object_exists(group_id, 'urr')) then
       this % urr_present = .true.

src/physics.F90

@@ -310,6 +310,7 @@ contains
     integer,        intent(in)    :: i_nuc_mat
 
     integer :: i
+    integer :: j
     integer :: i_temp
     integer :: i_grid
     real(8) :: f
@@ -378,11 +379,10 @@ contains
       ! =======================================================================
       ! INELASTIC SCATTERING
 
-      ! note that indexing starts from 2 since nuc % reactions(1) is elastic
-      ! scattering
-      i = 1
+      j = 0
       do while (prob < cutoff)
-        i = i + 1
+        j = j + 1
+        i = nuc % index_inelastic_scatter(j)
 
         ! Check to make sure inelastic scattering reaction sampled
         if (i > size(nuc % reactions)) then
@@ -391,24 +391,14 @@ contains
                &// trim(nuc % name))
         end if
 
-        associate (rx => nuc % reactions(i))
-          ! Skip fission reactions
-          if (rx % MT == N_FISSION .or. rx % MT == N_F .or. rx % MT == N_NF &
-               .or. rx % MT == N_2NF .or. rx % MT == N_3NF) cycle
-
-          ! Some materials have gas production cross sections with MT > 200 that
-          ! are duplicates. Also MT=4 is total level inelastic scattering which
-          ! should be skipped
-          if (rx % MT >= 200 .or. rx % MT == N_LEVEL) cycle
-
-          associate (xs => rx % xs(i_temp))
+        associate (rx => nuc % reactions(i), &
+                   xs => nuc % reactions(i) % xs(i_temp))
             ! if energy is below threshold for this reaction, skip it
             if (i_grid < xs % threshold) cycle
 
             ! add to cumulative probability
             prob = prob + ((ONE - f)*xs % value(i_grid - xs % threshold + 1) &
                  + f*(xs % value(i_grid - xs % threshold + 2)))
-          end associate
         end associate
       end do