diff --git a/src/ace.F90 b/src/ace.F90
index 584254149d..cb5beecb52 100644
--- a/src/ace.F90
+++ b/src/ace.F90
@@ -704,7 +704,8 @@ contains
     integer :: LXS       ! location of cross-section locators
     integer :: LOCA      ! location of cross-section for given MT
     integer :: IE        ! reaction's starting index on energy grid
-    integer :: NE        ! number of energies for reaction
+    integer :: NE        ! number of energies
+    integer :: NR        ! number of interpolation regions
     type(Reaction), pointer :: rxn => null()
     type(ListInt) :: MTs
 
@@ -755,18 +756,40 @@ contains
       rxn % multiplicity  = abs(nint(XSS(JXS5 + i - 1)))
       rxn % scatter_in_cm = (nint(XSS(JXS5 + i - 1)) < 0)
 
-      ! If multiplicity is energy-dependent (absolute value > 100), set it based
-      ! on the MT value
+      ! Read energy-dependent multiplicities
       if (rxn % multiplicity > 100) then
-        if (any(rxn%MT == [11, 16, 24, 30, 41])) then
-          rxn % multiplicity = 2
-        elseif (any(rxn%MT == [17, 25, 42])) then
-          rxn % multiplicity = 3
-        elseif (rxn%MT == 37) then
-          rxn % multiplicity = 4
-        else
-          rxn % multiplicity = 1
+        ! Set flag and allocate space for Tab1 to store yield
+        rxn % multiplicity_with_E = .true.
+        allocate(rxn % multiplicity_E)
+        
+        XSS_index = JXS(11) + rxn % multiplicity - 101
+        NR = nint(XSS(XSS_index))
+        rxn % multiplicity_E % n_regions = NR
+        
+        ! allocate space for ENDF interpolation parameters
+        if (NR > 0) then
+          allocate(rxn % multiplicity_E % nbt(NR))
+          allocate(rxn % multiplicity_E % int(NR))
         end if
+        
+        ! read ENDF interpolation parameters
+        XSS_index = XSS_index + 1
+        if (NR > 0) then
+          rxn % multiplicity_E % nbt = get_int(NR)
+          rxn % multiplicity_E % int = get_int(NR)
+        end if
+        
+        ! allocate space for yield data
+        XSS_index = XSS_index + 2*NR
+        NE = nint(XSS(XSS_index))
+        rxn % multiplicity_E % n_pairs = NE
+        allocate(rxn % multiplicity_E % x(NE))
+        allocate(rxn % multiplicity_E % y(NE))
+        
+        ! read yield data
+        XSS_index = XSS_index + 1
+        rxn % multiplicity_E % x = get_real(NE)
+        rxn % multiplicity_E % y = get_real(NE)
       end if
 
       ! read starting energy index
diff --git a/src/ace_header.F90 b/src/ace_header.F90
index afed796e29..f1cdd057f0 100644
--- a/src/ace_header.F90
+++ b/src/ace_header.F90
@@ -51,8 +51,10 @@ module ace_header
     integer :: MT                      ! ENDF MT value
     real(8) :: Q_value                 ! Reaction Q value
     integer :: multiplicity            ! Number of secondary particles released
+    type(Tab1), pointer :: multiplicity_E => null() ! Energy-dependent neutron yield
     integer :: threshold               ! Energy grid index of threshold
     logical :: scatter_in_cm           ! scattering system in center-of-mass?
+    logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity
     real(8), allocatable :: sigma(:)   ! Cross section values
     logical :: has_angle_dist          ! Angle distribution present?
     logical :: has_energy_dist         ! Energy distribution present?
diff --git a/src/physics.F90 b/src/physics.F90
index 860a20745c..b472b36ca0 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -1291,6 +1291,7 @@ contains
     real(8) :: E_in        ! incoming energy
     real(8) :: E_cm        ! outgoing energy in center-of-mass
     real(8) :: Q           ! Q-value of reaction
+    real(8) :: yield       ! neutron yield
 
     ! copy energy of neutron
     E_in = E
@@ -1330,8 +1331,13 @@ contains
     ! change direction of particle
     uvw = rotate_angle(uvw, mu)
 
-    ! change weight of particle based on multiplicity
-    wgt = rxn % multiplicity * wgt
+    ! change weight of particle based on yield
+    if (rxn % multiplicity_with_E) then
+      yield = interpolate_tab1(rxn % multiplicity_E, E_in)
+    else
+      yield = rxn % multiplicity
+    end if
+    wgt = yield * wgt
 
   end subroutine inelastic_scatter