diff --git a/src/cross_section.f90 b/src/cross_section.f90 index 82cbfd02c0..7e63eeed1e 100644 --- a/src/cross_section.f90 +++ b/src/cross_section.f90 @@ -1,7 +1,8 @@ module cross_section use constants - use cross_section_header, only: Nuclide, Reaction, SAB_Table, XsListing + use cross_section_header, only: Nuclide, Reaction, SAB_Table, XsListing, & + DistEnergy use datatypes, only: dict_create, dict_add_key, dict_get_key, & dict_has_key, dict_delete, dict_keys use datatypes_header, only: DictionaryCI, ListKeyValueCI @@ -841,21 +842,11 @@ contains type(Nuclide), pointer :: nuc integer :: LED ! location of energy distribution locators - integer :: LDIS ! location of all energy distributions integer :: LOCC ! location of energy distributions for given MT - integer :: LNW ! location of next energy distribution if multiple - integer :: LAW ! secondary energy distribution law - integer :: NR ! number of interpolation regions - integer :: NE ! number of incoming energies - integer :: IDAT ! location of first energy distribution for given MT - integer :: loc ! locator - integer :: length ! length of data to allocate - integer :: length_interp_data ! length of interpolation data - integer :: i ! loop index - type(Reaction), pointer :: rxn => null() + integer :: i ! loop index + type(Reaction), pointer :: rxn => null() LED = JXS(10) - LDIS = JXS(11) ! Loop over all reactions do i = 1, NXS(5) @@ -864,55 +855,94 @@ contains ! find location of energy distribution data LOCC = XSS(LED + i - 1) - - LNW = XSS(LDIS + LOCC - 1) - LAW = XSS(LDIS + LOCC) - IDAT = XSS(LDIS + LOCC + 1) - NR = XSS(LDIS + LOCC + 2) - rxn % edist % law = LAW - rxn % edist % n_interp = NR - ! allocate space for ENDF interpolation parameters - if (NR > 0) then - allocate(rxn % edist % nbt(NR)) - allocate(rxn % edist % int(NR)) - end if - - ! read ENDF interpolation parameters - XSS_index = LDIS + LOCC + 3 - if (NR > 0) then - rxn % edist % nbt = get_real(NR) - rxn % edist % int = get_real(NR) - end if + ! allocate energy distribution + allocate(rxn % edist) - ! allocate space for law validity data - NE = XSS(LDIS + LOCC + 3 + 2*NR) - allocate(rxn % edist % energy(NE)) - allocate(rxn % edist % pvalid(NE)) - - length_interp_data = 5 + 2*(NR + NE) - - ! read law validity data - XSS_index = LDIS + LOCC + 4 + 2*NR - rxn % edist % energy = get_real(NE) - rxn % edist % pvalid = get_real(NE) - - ! Set index to beginning of IDAT array - loc = LDIS + IDAT - 2 - - ! determine length of energy distribution - length = length_energy_dist(loc, LAW, LOCC, length_interp_data) - - ! allocate secondary energy distribution array - allocate(rxn % edist % data(length)) - - ! read secondary energy distribution - XSS_index = loc + 1 - rxn % edist % data = get_real(length) + ! read data for energy distribution + call get_energy_dist(rxn % edist, LOCC) end do end subroutine read_energy_dist +!=============================================================================== +! GET_ENERGY_DIST reads in data for a single law for an energy distribution and +! calls itself recursively if there are multiple energy distributions for a +! single reaction +!=============================================================================== + + recursive subroutine get_energy_dist(edist, loc_law) + + type(DistEnergy), pointer :: edist ! energy distribution + integer, intent(in) :: loc_law ! locator for data + + integer :: LDIS ! location of all energy distributions + integer :: LNW ! location of next energy distribution if multiple + integer :: LAW ! secondary energy distribution law + integer :: NR ! number of interpolation regions + integer :: NE ! number of incoming energies + integer :: IDAT ! location of first energy distribution for given MT + integer :: loc ! locator + integer :: length ! length of data to allocate + integer :: length_interp_data ! length of interpolation data + + LDIS = JXS(11) + + ! locator for next law and information on this law + LNW = XSS(LDIS + loc_law - 1) + LAW = XSS(LDIS + loc_law) + IDAT = XSS(LDIS + loc_law + 1) + NR = XSS(LDIS + loc_law + 2) + edist % law = LAW + edist % n_interp = NR + + ! allocate space for ENDF interpolation parameters + if (NR > 0) then + allocate(edist % nbt(NR)) + allocate(edist % int(NR)) + end if + + ! read ENDF interpolation parameters + XSS_index = LDIS + loc_law + 3 + if (NR > 0) then + edist % nbt = get_real(NR) + edist % int = get_real(NR) + end if + + ! allocate space for law validity data + NE = XSS(LDIS + loc_law + 3 + 2*NR) + edist % n_energy = NE + allocate(edist % energy(NE)) + allocate(edist % pvalid(NE)) + + length_interp_data = 5 + 2*(NR + NE) + + ! read law validity data + XSS_index = LDIS + loc_law + 4 + 2*NR + edist % energy = get_real(NE) + edist % pvalid = get_real(NE) + + ! Set index to beginning of IDAT array + loc = LDIS + IDAT - 2 + + ! determine length of energy distribution + length = length_energy_dist(loc, LAW, loc_law, length_interp_data) + + ! allocate secondary energy distribution array + allocate(edist % data(length)) + + ! read secondary energy distribution + XSS_index = loc + 1 + edist % data = get_real(length) + + ! read next energy distribution if present + if (LNW > 0) then + allocate(edist % next) + call get_energy_dist(edist % next, LNW) + end if + + end subroutine get_energy_dist + !=============================================================================== ! LENGTH_ENERGY_DIST determines how many values are contained in an LDAT energy ! distribution array based on the secondary energy law and location in XSS diff --git a/src/cross_section_header.f90 b/src/cross_section_header.f90 index 3ea2a2963e..ccbb19e8ca 100644 --- a/src/cross_section_header.f90 +++ b/src/cross_section_header.f90 @@ -31,6 +31,10 @@ module cross_section_header real(8), allocatable :: energy(:) ! energy grid for law validity real(8), allocatable :: pvalid(:) ! probability of law validity real(8), allocatable :: data(:) ! energy distribution data + + ! For reactions that may have multiple energy distributions such as (n.2n), + ! this pointer allows multiple laws to be stored + type(DistEnergy), pointer :: next => null() end type DistEnergy !=============================================================================== @@ -39,15 +43,15 @@ module cross_section_header !=============================================================================== type Reaction - integer :: MT ! ENDF MT value - real(8) :: Q_value ! Reaction Q value - integer :: TY ! Number of neutrons released - integer :: IE ! Starting energy grid index - real(8), allocatable :: sigma(:) ! Cross section values - logical :: has_angle_dist ! Angle distribution present? - logical :: has_energy_dist ! Energy distribution present? - type(DistAngle) :: adist ! Secondary angular distribution - type(DistEnergy) :: edist ! Secondary energy distribution + integer :: MT ! ENDF MT value + real(8) :: Q_value ! Reaction Q value + integer :: TY ! Number of neutrons released + integer :: IE ! Starting energy grid index + real(8), allocatable :: sigma(:) ! Cross section values + logical :: has_angle_dist ! Angle distribution present? + logical :: has_energy_dist ! Energy distribution present? + type(DistAngle) :: adist ! Secondary angular distribution + type(DistEnergy), pointer :: edist ! Secondary energy distribution end type Reaction !=============================================================================== @@ -103,7 +107,7 @@ module cross_section_header integer :: n_precursor real(8), allocatable :: nu_d_data(:) real(8), allocatable :: nu_d_precursor_data(:) - type(DistEnergy), allocatable :: nu_d_edist(:) + type(DistEnergy), pointer :: nu_d_edist(:) => null() ! Unresolved resonance data logical :: urr_present diff --git a/src/physics.f90 b/src/physics.f90 index 7bcd4b806d..7c78232b85 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -1025,7 +1025,8 @@ contains real(8) :: yield ! delayed neutron precursor yield real(8) :: prob ! cumulative probability logical :: actual_event ! did fission actually occur? (no survival biasing) - type(Nuclide), pointer :: nuc + type(Nuclide), pointer :: nuc + type(DistEnergy), pointer :: edist => null() ! Get pointer to nuclide nuc => nuclides(index_nuclide) @@ -1130,11 +1131,12 @@ contains ! sample from energy distribution for group j law = nuc % nu_d_edist(j) % law + edist => nuc % nu_d_edist(j) do if (law == 44 .or. law == 61) then - call sample_energy(nuc%nu_d_edist(j), E, E_out, mu) + call sample_energy(edist, E, E_out, mu) else - call sample_energy(nuc%nu_d_edist(j), E, E_out) + call sample_energy(edist, E, E_out) end if ! resample if energy is >= 20 MeV if (E_out < 20) exit @@ -1425,7 +1427,7 @@ contains subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q) - type(DistEnergy), intent(inout) :: edist + type(DistEnergy), pointer :: edist real(8), intent(in) :: E_in real(8), intent(out) :: E_out real(8), intent(inout), optional :: mu_out