Added most logic for unresolved resonance probability tables. Does not work yet.

src/constants.F90

@@ -208,6 +208,15 @@ module constants
        ASCII  = 1, & ! ASCII cross section file
        BINARY = 2    ! Binary cross section file
 
+  ! Probability table parameters
+  integer, parameter :: &
+       URR_CUM_PROB = 1, &
+       URR_TOTAL    = 2, &
+       URR_ELASTIC  = 3, &
+       URR_FISSION  = 4, &
+       URR_N_GAMMA  = 5, &
+       URR_HEATING  = 6
+
   ! Maximum number of partial fission reactions
   integer, parameter :: PARTIAL_FISSION_MAX = 4
 

src/cross_section.F90

@@ -1107,6 +1107,25 @@ contains
        nuc % urr_data % multiply_smooth = .true.
     end if
 
+    ! if the inelastic competition flag indicates that the inelastic cross
+    ! section should be determined from a normal reaction cross section, we need
+    ! to set up a pointer to that reaction
+    nuc % urr_inelastic = NONE
+    if (nuc % urr_data % inelastic_flag > 0) then
+       do i = 1, nuc % n_reaction
+          if (nuc % reactions(i) % MT == nuc % urr_data % inelastic_flag) then
+             nuc % urr_inelastic = i
+          end if
+       end do
+
+       ! Abort if no corresponding inelastic reaction was found
+       if (nuc % urr_inelastic == NONE) then
+          message = "Could not find inelastic reaction specified on " &
+               // "unresolved resonance probability table."
+          call fatal_error()
+       end if
+    end if
+
     ! allocate incident energies and probability tables
     N = nuc % urr_data % n_energy
     M = nuc % urr_data % n_prob

src/cross_section_header.F90

@@ -113,6 +113,7 @@ module cross_section_header
 
      ! Unresolved resonance data
      logical                :: urr_present
+     integer                :: urr_inelastic
      type(UrrData), pointer :: urr_data => null()
 
      ! Reactions

src/global.F90

@@ -73,6 +73,9 @@ module global
   integer              :: n_grid    ! number of points on unionized grid
   real(8), allocatable :: e_grid(:) ! energies on unionized grid
 
+  ! Unreoslved resonance probablity tables
+  logical :: urr_ptables_on = .false.
+
   ! ============================================================================
   ! TALLY-RELATED VARIABLES
 

src/physics.F90

@@ -1,7 +1,7 @@
 module physics
 
   use constants
-  use cross_section_header, only: Nuclide, Reaction, DistEnergy
+  use cross_section_header, only: Nuclide, Reaction, DistEnergy, UrrData
   use endf,                 only: reaction_name, is_fission, is_scatter
   use error,                only: fatal_error, warning
   use fission,              only: nu_total, nu_prompt, nu_delayed
@@ -279,6 +279,16 @@ contains
 
     if (index_sab > 0) call calculate_sab_xs(p, index_nuclide, index_sab)
 
+    ! if the particle is in the unresolved resonance range and there are
+    ! probability tables, we need to determine cross sections from the table
+
+    if (urr_ptables_on .and. nuc % urr_present) then
+       if (p % E > nuc % urr_data % energy(1) .and. &
+            p % E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then
+          call calculate_urr_xs(p, index_nuclide)
+       end if
+    end if
+
     ! Set last evaluated energy
     micro_xs(index_nuclide) % last_E = p % E
 
@@ -368,6 +378,113 @@ contains
 
   end subroutine calculate_sab_xs
 
+!===============================================================================
+! CALCULATE_URR_XS determines cross sections in the unresolved resonance range
+! from probability tables
+!===============================================================================
+
+  subroutine calculate_urr_xs(p, index_nuclide)
+
+    type(Particle), pointer :: p
+    integer, intent(in)     :: index_nuclide ! index into nuclides array
+
+    integer :: i_energy   ! index for energy
+    integer :: i_table    ! index for table
+    real(8) :: f          ! interpolation factor
+    real(8) :: r          ! pseudo-random number
+    real(8) :: elastic    ! smooth elastic cross section
+    real(8) :: absorption ! smooth absorption cross section
+    real(8) :: fission    ! smooth fission cross section
+    real(8) :: inelastic  ! smooth inelastic cross section
+    type(UrrData),  pointer :: urr => null()
+    type(Nuclide),  pointer :: nuc => null()
+    type(Reaction), pointer :: rxn => null()
+
+    ! copy cross-sections already calculated
+    elastic    = micro_xs(index_nuclide) % elastic
+    absorption = micro_xs(index_nuclide) % absorption
+    fission    = micro_xs(index_nuclide) % fission
+
+    ! get pointer to probability table
+    nuc => nuclides(index_nuclide)
+    urr => nuc % urr_data
+
+    ! determine energy table
+    i_energy = 1
+    do
+       if (p % E > urr % energy(i_energy)) exit
+       i_energy = i_energy + 1
+    end do
+
+    ! determine interpolation factor on table
+    f = (p % E - urr % energy(i_energy)) / &
+         (urr % energy(i_energy + 1) - urr % energy(i_energy))
+
+    ! sample probability table using the cumulative distribution
+    r = prn()
+    i_table = 1
+    do
+       if (urr % prob(i_energy, URR_CUM_PROB, i_table) > r) exit
+       i_table = i_table + 1
+    end do
+
+    ! determine elastic, fission, and capture cross sections from probability
+    ! table
+    if (urr % interp == LINEAR_LINEAR) then
+       micro_xs(index_nuclide) % elastic = &
+            (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_table) + &
+            f * urr % prob(i_energy + 1, URR_ELASTIC, i_table)
+       micro_xs(index_nuclide) % fission = &
+            (ONE - f) * urr % prob(i_energy, URR_FISSION, i_table) + &
+            f * urr % prob(i_energy + 1, URR_FISSION, i_table)
+       micro_xs(index_nuclide) % absorption = &
+            micro_xs(index_nuclide) % fission + &
+            (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_table) + &
+            f * urr % prob(i_energy + 1, URR_N_GAMMA, i_table)
+    elseif (urr % interp == LOG_LOG) then
+       message = "Log-log interpolation on probability table not yet supported."
+       call fatal_error()
+    end if
+
+    ! Determine treatment of inelastic scattering
+    inelastic = ZERO
+    if (urr % inelastic_flag > 0) then
+       ! Get pointer to inelastic scattering reaction
+       rxn => nuc % reactions(nuc % urr_inelastic)
+
+       ! Get index on energy grid and interpolation factor
+       i_energy = micro_xs(index_nuclide) % index_grid
+       f = micro_xs(index_nuclide) % interp_factor
+
+       ! Determine inelastic scattering cross section
+       if (i_energy >= rxn % IE) then
+          inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%IE + 1) + &
+               f * rxn % sigma(i_energy - rxn%IE + 2)
+       end if
+    end if
+
+    ! Multiply by smooth cross-section if needed
+    if (urr % multiply_smooth) then
+       micro_xs(index_nuclide) % elastic = elastic * &
+            micro_xs(index_nuclide) % elastic
+       micro_xs(index_nuclide) % absorption = absorption * &
+            micro_xs(index_nuclide) % absorption
+       micro_xs(index_nuclide) % fission = fission * &
+            micro_xs(index_nuclide) % fission
+    end if
+
+    ! Determine nu-fission cross section
+    if (nuc % fissionable) then
+       micro_xs(index_nuclide) % nu_fission = nu_total(nuc, p % E) * &
+            micro_xs(index_nuclide) % fission
+    end if
+
+    ! Calculate total cross section as sum of partials
+    micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % elastic + &
+         inelastic + micro_xs(index_nuclide) % absorption
+
+  end subroutine calculate_urr_xs
+
 !===============================================================================
 ! FIND_ENERGY_INDEX determines the index on the union energy grid and the
 ! interpolation factor for a particle at a certain energy