diff --git a/src/ace_header.F90 b/src/ace_header.F90 index c881a559dc..b0fec99b6c 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -117,7 +117,6 @@ module ace_header character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c character(16) :: scheme ! target velocity sampling scheme integer :: n_grid_0K ! number of 0K energy grid points - integer, allocatable :: grid_index_0K(:) ! pointers to union grid real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section @@ -372,9 +371,6 @@ module ace_header if (allocated(this % grid_index)) & deallocate(this % grid_index) - - if (allocated(this % grid_index_0K)) & - deallocate(this % grid_index_0K) if (allocated(this % energy)) & deallocate(this % energy, this % total, this % elastic, & diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 1c6ed8ae8f..0620cf11f8 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -527,37 +527,26 @@ contains real(8) :: xs_out ! 0K xs at trial energy ! Determine index on nuclide energy grid - select case (grid_method) - case (GRID_UNION) - ! If we're using the unionized grid with pointers, finding the index on - ! the nuclide energy grid is as simple as looking up the pointer - - call find_energy_index(E) - i_grid = nuc % grid_index_0K(union_grid_index) - - case (GRID_NUCLIDE) - ! If we're not using the unionized grid, we have to do a binary search on - ! the nuclide energy grid in order to determine which points to - ! interpolate between - - if (E < nuc % energy(1)) then - i_grid = 1 - elseif (E > nuc % energy(nuc % n_grid)) then - i_grid = nuc % n_grid - 1 - else - i_grid = binary_search(nuc % energy, nuc % n_grid, E) - end if - - end select + if (E < nuc % energy_0K(1)) then + i_grid = 1 + elseif (E > nuc % energy_0K(nuc % n_grid_0K)) then + i_grid = nuc % n_grid_0K - 1 + else + i_grid = binary_search(nuc % energy_0K, nuc % n_grid_0K, E) + end if ! check for rare case where two energy points are the same - if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) i_grid = i_grid + 1 + if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) then + i_grid = i_grid + 1 + end if ! calculate interpolation factor - f = (E - nuc%energy_0K(i_grid))/(nuc%energy_0K(i_grid+1) - nuc%energy_0K(i_grid)) + f = (E - nuc % energy_0K(i_grid)) & + & / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid)) ! Calculate microscopic nuclide elastic cross section - xs_out = (ONE - f) * nuc % elastic_0K(i_grid) + f * nuc % elastic_0K(i_grid+1) + xs_out = (ONE - f) * nuc % elastic_0K(i_grid) & + & + f * nuc % elastic_0K(i_grid + 1) end function elastic_xs_0K diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 642b86b100..a38a305d3e 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -147,24 +147,6 @@ contains end if nuc % grid_index(j) = index_e - 1 end do - - ! set pointers for 0K energy grid to the unionized grid - if (nuc % resonant) then - allocate(nuc % grid_index_0K(n_grid)) - - index_e = 1 - energy = nuc % energy_0K(index_e) - - do j = 1, n_grid - union_energy = e_grid(j) - if (union_energy >= energy .and. index_e < nuc % n_grid_0K) then - index_e = index_e + 1 - energy = nuc % energy_0K(index_e) - end if - nuc % grid_index_0K(j) = index_e - 1 - end do - end if - end do end subroutine grid_pointers diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 9b8b1a3c11..27995597ce 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -842,7 +842,7 @@ contains end if ! check that E_min is non-negative - if (E_min < ZERO) then + if (nuclides_0K(i) % E_min < ZERO) then message = "Lower resonance scattering energy bound is negative" call fatal_error() end if @@ -853,7 +853,7 @@ contains end if ! check that E_max is not less than E_min - if (E_max < E_min) then + if (nuclides_0K(i) % E_max < nuclides_0K(i) % E_min) then message = "Lower resonance scattering energy bound exceeds upper" call fatal_error() end if diff --git a/src/physics.F90 b/src/physics.F90 index 32d3175673..e0fbf80d25 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -835,13 +835,26 @@ contains E_red = sqrt((awr * E) / kT) E_low = (((E_red - 4.0_8)**2) * kT) / awr E_up = (((E_red + 4.0_8)**2) * kT) / awr - + ! find lower and upper energy bound indices - call find_energy_index(E_low) - i_E_low = nuc % grid_index_0K(union_grid_index) - call find_energy_index(E_up) - i_E_up = nuc % grid_index_0K(union_grid_index) - + ! lower index + if (E_low < nuc % energy_0K(1)) then + i_E_low = 1 + elseif (E_low > nuc % energy_0K(nuc % n_grid_0K)) then + i_E_low = nuc % n_grid_0K - 1 + else + i_E_low = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_low) + end if + + ! upper index + if (E_up < nuc % energy_0K(1)) then + i_E_up = 1 + elseif (E_up > nuc % energy_0K(nuc % n_grid_0K)) then + i_E_up = nuc % n_grid_0K - 1 + else + i_E_up = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_up) + end if + ! interpolate xs since we're not exactly at the energy indices xs_low = nuc % elastic_0K(i_E_low) m = (nuc % elastic_0K(i_E_low + 1) - xs_low) & @@ -878,10 +891,23 @@ contains E_up = (((E_red + 4.0_8)**2) * kT) / awr ! find lower and upper energy bound indices - call find_energy_index(E_low) - i_E_low = nuc % grid_index_0K(union_grid_index) - call find_energy_index(E_up) - i_E_up = nuc % grid_index_0K(union_grid_index) + ! lower index + if (E_low < nuc % energy_0K(1)) then + i_E_low = 1 + elseif (E_low > nuc % energy_0K(nuc % n_grid_0K)) then + i_E_low = nuc % n_grid_0K - 1 + else + i_E_low = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_low) + end if + + ! upper index + if (E_up < nuc % energy_0K(1)) then + i_E_up = 1 + elseif (E_up > nuc % energy_0K(nuc % n_grid_0K)) then + i_E_up = nuc % n_grid_0K - 1 + else + i_E_up = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_up) + end if ! interpolate xs CDF since we're not exactly at the energy indices ! cdf value at lower bound attainable energy