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