diff --git a/src/ace.F90 b/src/ace.F90 index dbed66757..4cd53a251 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -715,7 +715,6 @@ contains integer :: IE ! reaction's starting index on energy grid integer :: NE ! number of energies integer :: NR ! number of interpolation regions - type(Reaction), pointer :: rxn type(ListInt) :: MTs LMT = JXS(3) @@ -733,14 +732,15 @@ contains ! Store elastic scattering cross-section on reaction one -- note that the ! sigma array is not allocated or stored for elastic scattering since it is ! already stored in nuc % elastic - rxn => nuc % reactions(1) - rxn % MT = 2 - rxn % Q_value = ZERO - rxn % multiplicity = 1 - rxn % threshold = 1 - rxn % scatter_in_cm = .true. - rxn % has_angle_dist = .false. - rxn % has_energy_dist = .false. + associate (rxn => nuc % reactions(1)) + rxn % MT = 2 + rxn % Q_value = ZERO + rxn % multiplicity = 1 + rxn % threshold = 1 + rxn % scatter_in_cm = .true. + rxn % has_angle_dist = .false. + rxn % has_energy_dist = .false. + end associate ! Add contribution of elastic scattering to total cross section nuc % total = nuc % total + nuc % elastic @@ -753,64 +753,64 @@ contains i_fission = 0 do i = 1, NMT - rxn => nuc % reactions(i+1) + associate (rxn => nuc % reactions(i+1)) + ! set defaults + rxn % has_angle_dist = .false. + rxn % has_energy_dist = .false. - ! set defaults - rxn % has_angle_dist = .false. - rxn % has_energy_dist = .false. + ! read MT number, Q-value, and neutrons produced + rxn % MT = int(XSS(LMT + i - 1)) + rxn % Q_value = XSS(JXS4 + i - 1) + rxn % multiplicity = abs(nint(XSS(JXS5 + i - 1))) + rxn % scatter_in_cm = (nint(XSS(JXS5 + i - 1)) < 0) - ! read MT number, Q-value, and neutrons produced - rxn % MT = int(XSS(LMT + i - 1)) - rxn % Q_value = XSS(JXS4 + i - 1) - rxn % multiplicity = abs(nint(XSS(JXS5 + i - 1))) - rxn % scatter_in_cm = (nint(XSS(JXS5 + i - 1)) < 0) + ! Read energy-dependent multiplicities + if (rxn % multiplicity > 100) then + ! Set flag and allocate space for Tab1 to store yield + rxn % multiplicity_with_E = .true. + allocate(rxn % multiplicity_E) - ! Read energy-dependent multiplicities - if (rxn % multiplicity > 100) then - ! 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 - 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 - ! allocate space for ENDF interpolation parameters - if (NR > 0) then - allocate(rxn % multiplicity_E % nbt(NR)) - allocate(rxn % multiplicity_E % int(NR)) + ! 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 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 + ! read starting energy index + LOCA = int(XSS(LXS + i - 1)) + IE = int(XSS(JXS7 + LOCA - 1)) + rxn % threshold = IE - ! 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 - LOCA = int(XSS(LXS + i - 1)) - IE = int(XSS(JXS7 + LOCA - 1)) - rxn % threshold = IE - - ! read number of energies cross section values - NE = int(XSS(JXS7 + LOCA)) - allocate(rxn % sigma(NE)) - XSS_index = JXS7 + LOCA + 1 - rxn % sigma = get_real(NE) + ! read number of energies cross section values + NE = int(XSS(JXS7 + LOCA)) + allocate(rxn % sigma(NE)) + XSS_index = JXS7 + LOCA + 1 + rxn % sigma = get_real(NE) + end associate end do ! Create set of MT values @@ -821,56 +821,57 @@ contains ! Create total, absorption, and fission cross sections do i = 2, size(nuc % reactions) - rxn => nuc % reactions(i) - IE = rxn % threshold - NE = size(rxn % sigma) + associate (rxn => nuc % reactions(i)) + IE = rxn % threshold + NE = size(rxn % sigma) - ! Skip total inelastic level scattering, gas production cross sections - ! (MT=200+), etc. - if (rxn % MT == N_LEVEL) cycle - if (rxn % MT > N_5N2P .and. rxn % MT < N_P0) cycle + ! Skip total inelastic level scattering, gas production cross sections + ! (MT=200+), etc. + if (rxn % MT == N_LEVEL) cycle + if (rxn % MT > N_5N2P .and. rxn % MT < N_P0) cycle - ! Skip level cross sections if total is available - if (rxn % MT >= N_P0 .and. rxn % MT <= N_PC .and. MTs % contains(N_P)) cycle - if (rxn % MT >= N_D0 .and. rxn % MT <= N_DC .and. MTs % contains(N_D)) cycle - if (rxn % MT >= N_T0 .and. rxn % MT <= N_TC .and. MTs % contains(N_T)) cycle - if (rxn % MT >= N_3HE0 .and. rxn % MT <= N_3HEC .and. MTs % contains(N_3HE)) cycle - if (rxn % MT >= N_A0 .and. rxn % MT <= N_AC .and. MTs % contains(N_A)) cycle - if (rxn % MT >= N_2N0 .and. rxn % MT <= N_2NC .and. MTs % contains(N_2N)) cycle + ! Skip level cross sections if total is available + if (rxn % MT >= N_P0 .and. rxn % MT <= N_PC .and. MTs % contains(N_P)) cycle + if (rxn % MT >= N_D0 .and. rxn % MT <= N_DC .and. MTs % contains(N_D)) cycle + if (rxn % MT >= N_T0 .and. rxn % MT <= N_TC .and. MTs % contains(N_T)) cycle + if (rxn % MT >= N_3HE0 .and. rxn % MT <= N_3HEC .and. MTs % contains(N_3HE)) cycle + if (rxn % MT >= N_A0 .and. rxn % MT <= N_AC .and. MTs % contains(N_A)) cycle + if (rxn % MT >= N_2N0 .and. rxn % MT <= N_2NC .and. MTs % contains(N_2N)) cycle - ! Add contribution to total cross section - nuc % total(IE:IE+NE-1) = nuc % total(IE:IE+NE-1) + rxn % sigma + ! Add contribution to total cross section + nuc % total(IE:IE+NE-1) = nuc % total(IE:IE+NE-1) + rxn % sigma - ! Add contribution to absorption cross section - if (is_disappearance(rxn % MT)) then - nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma - end if + ! Add contribution to absorption cross section + if (is_disappearance(rxn % MT)) then + nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma + end if - ! Information about fission reactions - if (rxn % MT == N_FISSION) then - allocate(nuc % index_fission(1)) - elseif (rxn % MT == N_F) then - allocate(nuc % index_fission(PARTIAL_FISSION_MAX)) - nuc % has_partial_fission = .true. - end if + ! Information about fission reactions + if (rxn % MT == N_FISSION) then + allocate(nuc % index_fission(1)) + elseif (rxn % MT == N_F) then + allocate(nuc % index_fission(PARTIAL_FISSION_MAX)) + nuc % has_partial_fission = .true. + end if - ! Add contribution to fission cross section - if (is_fission(rxn % MT)) then - nuc % fissionable = .true. - nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma + ! Add contribution to fission cross section + if (is_fission(rxn % MT)) then + nuc % fissionable = .true. + nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma - ! Also need to add fission cross sections to absorption - nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma + ! Also need to add fission cross sections to absorption + nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma - ! If total fission reaction is present, there's no need to store the - ! reaction cross-section since it was copied to nuc % fission - if (rxn % MT == N_FISSION) deallocate(rxn % sigma) + ! If total fission reaction is present, there's no need to store the + ! reaction cross-section since it was copied to nuc % fission + if (rxn % MT == N_FISSION) deallocate(rxn % sigma) - ! Keep track of this reaction for easy searching later - i_fission = i_fission + 1 - nuc % index_fission(i_fission) = i - nuc % n_fission = nuc % n_fission + 1 - end if + ! Keep track of this reaction for easy searching later + i_fission = i_fission + 1 + nuc % index_fission(i_fission) = i + nuc % n_fission = nuc % n_fission + 1 + end if + end associate end do ! Clear MTs set @@ -895,7 +896,6 @@ contains integer :: i ! index in reactions array integer :: j ! index over incoming energies integer :: length ! length of data array to allocate - type(Reaction), pointer :: rxn JXS8 = JXS(8) JXS9 = JXS(9) @@ -903,71 +903,72 @@ contains ! loop over all reactions with secondary neutrons -- NXS(5) does not include ! elastic scattering do i = 1, NXS(5) + 1 - rxn => nuc%reactions(i) + associate (rxn => nuc%reactions(i)) - ! find location of angular distribution - LOCB = int(XSS(JXS8 + i - 1)) - if (LOCB == -1) then - ! Angular distribution data are specified through LAWi = 44 in the DLW - ! block - cycle - elseif (LOCB == 0) then - ! No angular distribution data are given for this reaction, isotropic - ! scattering is asssumed (in CM if TY < 0 and in LAB if TY > 0) - cycle - end if - rxn % has_angle_dist = .true. - - ! allocate space for incoming energies and locations - NE = int(XSS(JXS9 + LOCB - 1)) - rxn % adist % n_energy = NE - allocate(rxn % adist % energy(NE)) - allocate(rxn % adist % type(NE)) - allocate(rxn % adist % location(NE)) - - ! read incoming energy grid and location of nucs - XSS_index = JXS9 + LOCB - rxn % adist % energy = get_real(NE) - rxn % adist % location = get_int(NE) - - ! determine dize of data block - length = 0 - do j = 1, NE - LC = rxn % adist % location(j) - if (LC == 0) then - ! isotropic - rxn % adist % type(j) = ANGLE_ISOTROPIC - elseif (LC > 0) then - ! 32 equiprobable bins - rxn % adist % type(j) = ANGLE_32_EQUI - length = length + 33 - elseif (LC < 0) then - ! tabular distribution - rxn % adist % type(j) = ANGLE_TABULAR - NP = int(XSS(JXS9 + abs(LC))) - length = length + 2 + 3*NP + ! find location of angular distribution + LOCB = int(XSS(JXS8 + i - 1)) + if (LOCB == -1) then + ! Angular distribution data are specified through LAWi = 44 in the DLW + ! block + cycle + elseif (LOCB == 0) then + ! No angular distribution data are given for this reaction, isotropic + ! scattering is asssumed (in CM if TY < 0 and in LAB if TY > 0) + cycle end if - end do + rxn % has_angle_dist = .true. - ! allocate angular distribution data and read - allocate(rxn % adist % data(length)) + ! allocate space for incoming energies and locations + NE = int(XSS(JXS9 + LOCB - 1)) + rxn % adist % n_energy = NE + allocate(rxn % adist % energy(NE)) + allocate(rxn % adist % type(NE)) + allocate(rxn % adist % location(NE)) - ! read angular distribution -- currently this does not actually parse the - ! angular distribution tables for each incoming energy, that must be done - ! on-the-fly - XSS_index = JXS9 + LOCB + 2 * NE - rxn % adist % data = get_real(length) + ! read incoming energy grid and location of nucs + XSS_index = JXS9 + LOCB + rxn % adist % energy = get_real(NE) + rxn % adist % location = get_int(NE) - ! change location pointers since they are currently relative to JXS(9) - LC = LOCB + 2 * NE + 1 - do j = 1, NE - ! For consistency, leave location as 0 if type is isotropic. - ! This is not necessary for current correctness, but can avoid - ! future issues - if (rxn % adist % location(j) /= 0) then - rxn % adist % location(j) = abs(rxn % adist % location(j)) - LC - end if - end do + ! determine dize of data block + length = 0 + do j = 1, NE + LC = rxn % adist % location(j) + if (LC == 0) then + ! isotropic + rxn % adist % type(j) = ANGLE_ISOTROPIC + elseif (LC > 0) then + ! 32 equiprobable bins + rxn % adist % type(j) = ANGLE_32_EQUI + length = length + 33 + elseif (LC < 0) then + ! tabular distribution + rxn % adist % type(j) = ANGLE_TABULAR + NP = int(XSS(JXS9 + abs(LC))) + length = length + 2 + 3*NP + end if + end do + + ! allocate angular distribution data and read + allocate(rxn % adist % data(length)) + + ! read angular distribution -- currently this does not actually parse the + ! angular distribution tables for each incoming energy, that must be done + ! on-the-fly + XSS_index = JXS9 + LOCB + 2 * NE + rxn % adist % data = get_real(length) + + ! change location pointers since they are currently relative to JXS(9) + LC = LOCB + 2 * NE + 1 + do j = 1, NE + ! For consistency, leave location as 0 if type is isotropic. + ! This is not necessary for current correctness, but can avoid + ! future issues + if (rxn % adist % location(j) /= 0) then + rxn % adist % location(j) = abs(rxn % adist % location(j)) - LC + end if + end do + end associate end do end subroutine read_angular_dist @@ -983,23 +984,23 @@ contains integer :: LED ! location of energy distribution locators integer :: LOCC ! location of energy distributions for given MT integer :: i ! loop index - type(Reaction), pointer :: rxn LED = JXS(10) ! Loop over all reactions do i = 1, NXS(5) - rxn => nuc % reactions(i+1) ! skip over elastic scattering - rxn % has_energy_dist = .true. + associate (rxn => nuc % reactions(i+1)) ! skip over elastic scattering + rxn % has_energy_dist = .true. - ! find location of energy distribution data - LOCC = int(XSS(LED + i - 1)) + ! find location of energy distribution data + LOCC = int(XSS(LED + i - 1)) - ! allocate energy distribution - allocate(rxn % edist) + ! allocate energy distribution + allocate(rxn % edist) - ! read data for energy distribution - call get_energy_dist(rxn % edist, LOCC) + ! read data for energy distribution + call get_energy_dist(rxn % edist, LOCC) + end associate end do end subroutine read_energy_dist diff --git a/src/ace_header.F90 b/src/ace_header.F90 index a10c6fe87..9eadca594 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -146,7 +146,7 @@ module ace_header ! Reactions integer :: n_reaction ! # of reactions - type(Reaction), pointer :: reactions(:) => null() + type(Reaction), allocatable :: reactions(:) type(DictIntInt) :: reaction_index ! map MT values to index in reactions ! array; used at tally-time @@ -341,11 +341,10 @@ module ace_header end if if (associated(this % urr_data)) then - call this % urr_data % clear() deallocate(this % urr_data) end if - if (associated(this % reactions)) then + if (allocated(this % reactions)) then do i = 1, size(this % reactions) call this % reactions(i) % clear() end do diff --git a/src/cross_section.F90 b/src/cross_section.F90 index b9c76b503..03bcefca8 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -379,7 +379,6 @@ contains logical :: same_nuc ! do we know the xs for this nuclide at this energy? type(UrrData), pointer :: urr type(Nuclide), pointer :: nuc - type(Reaction), pointer :: rxn micro_xs(i_nuclide) % use_ptable = .true. @@ -475,18 +474,17 @@ contains ! 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(i_nuclide) % index_grid f = micro_xs(i_nuclide) % interp_factor ! Determine inelastic scattering cross section - if (i_energy >= rxn % threshold) then - inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + & - f * rxn % sigma(i_energy - rxn%threshold + 2) - end if + associate (rxn => nuc % reactions(nuc % urr_inelastic)) + if (i_energy >= rxn % threshold) then + inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + & + f * rxn % sigma(i_energy - rxn%threshold + 2) + end if + end associate end if ! Multiply by smooth cross-section if needed diff --git a/src/output.F90 b/src/output.F90 index 2aa90985a..19ae6268b 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -330,7 +330,6 @@ contains integer :: size_energy ! memory used for a energy distributions (bytes) integer :: size_urr ! memory used for probability tables (bytes) character(11) :: law ! secondary energy distribution law - type(Reaction), pointer :: rxn type(UrrData), pointer :: urr ! set default unit for writing information @@ -359,32 +358,32 @@ contains ! Information on each reaction write(unit_,*) ' Reaction Q-value COM Law IE size(angle) size(energy)' do i = 1, nuc % n_reaction - rxn => nuc % reactions(i) + associate (rxn => nuc % reactions(i)) + ! Determine size of angle distribution + if (rxn % has_angle_dist) then + size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8 + else + size_angle = 0 + end if - ! Determine size of angle distribution - if (rxn % has_angle_dist) then - size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8 - else - size_angle = 0 - end if + ! Determine size of energy distribution and law + if (rxn % has_energy_dist) then + size_energy = size(rxn % edist % data) * 8 + law = to_str(rxn % edist % law) + else + size_energy = 0 + law = 'None' + end if - ! Determine size of energy distribution and law - if (rxn % has_energy_dist) then - size_energy = size(rxn % edist % data) * 8 - law = to_str(rxn % edist % law) - else - size_energy = 0 - law = 'None' - end if + write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') & + reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, & + law(1:4), rxn % threshold, size_angle, size_energy - write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') & - reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, & - law(1:4), rxn % threshold, size_angle, size_energy - - ! Accumulate data size - size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8 - size_angle_total = size_angle_total + size_angle - size_energy_total = size_energy_total + size_energy + ! Accumulate data size + size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8 + size_angle_total = size_angle_total + size_angle + size_energy_total = size_energy_total + size_energy + end associate end do ! Add memory required for summary reactions (total, absorption, fission, diff --git a/src/physics.F90 b/src/physics.F90 index c06919758..a01a3a30c 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -185,7 +185,6 @@ contains !=============================================================================== subroutine sample_fission(i_nuclide, i_reaction) - integer, intent(in) :: i_nuclide ! index in nuclides array integer, intent(out) :: i_reaction ! index in nuc % reactions array @@ -195,7 +194,6 @@ contains real(8) :: prob real(8) :: cutoff type(Nuclide), pointer :: nuc - type(Reaction), pointer :: rxn ! Get pointer to nuclide nuc => nuclides(i_nuclide) @@ -220,14 +218,15 @@ contains FISSION_REACTION_LOOP: do i = 1, nuc % n_fission i_reaction = nuc % index_fission(i) - rxn => nuc % reactions(i_reaction) - ! if energy is below threshold for this reaction, skip it - if (i_grid < rxn % threshold) cycle + associate (rxn => nuc % reactions(i_reaction)) + ! if energy is below threshold for this reaction, skip it + if (i_grid < rxn % threshold) cycle - ! add to cumulative probability - prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & - + f*(rxn%sigma(i_grid - rxn%threshold + 2))) + ! add to cumulative probability + prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + + f*(rxn%sigma(i_grid - rxn%threshold + 2))) + end associate ! Create fission bank sites if fission occurs if (prob > cutoff) exit FISSION_REACTION_LOOP @@ -312,11 +311,10 @@ contains real(8) :: f real(8) :: prob real(8) :: cutoff - type(Nuclide), pointer :: nuc - type(Reaction), pointer :: rxn real(8) :: uvw_new(3) ! outgoing uvw for iso-in-lab scattering real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering real(8) :: phi ! azimuthal angle for iso-in-lab scattering + type(Nuclide), pointer :: nuc ! copy incoming direction uvw_old(:) = p % coord(1) % uvw @@ -343,11 +341,8 @@ contains p % E, p % coord(1) % uvw, p % mu) else - ! get pointer to elastic scattering reaction - rxn => nuc % reactions(1) - ! Perform collision physics for elastic scattering - call elastic_scatter(i_nuclide, rxn, & + call elastic_scatter(i_nuclide, nuc % reactions(1), & p % E, p % coord(1) % uvw, p % mu, p % wgt) end if @@ -370,28 +365,28 @@ contains &// trim(nuc % name)) end if - rxn => nuc % reactions(i) + associate (rxn => nuc % reactions(i)) + ! Skip fission reactions + if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF & + .or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) cycle - ! Skip fission reactions - if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF & - .or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) cycle + ! some materials have gas production cross sections with MT > 200 that + ! are duplicates. Also MT=4 is total level inelastic scattering which + ! should be skipped + if (rxn % MT >= 200 .or. rxn % MT == N_LEVEL) cycle - ! some materials have gas production cross sections with MT > 200 that - ! are duplicates. Also MT=4 is total level inelastic scattering which - ! should be skipped - if (rxn % MT >= 200 .or. rxn % MT == N_LEVEL) cycle + ! if energy is below threshold for this reaction, skip it + if (i_grid < rxn % threshold) cycle - ! if energy is below threshold for this reaction, skip it - if (i_grid < rxn % threshold) cycle - - ! add to cumulative probability - prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & - + f*(rxn%sigma(i_grid - rxn%threshold + 2))) + ! add to cumulative probability + prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + + f*(rxn%sigma(i_grid - rxn%threshold + 2))) + end associate end do ! Perform collision physics for inelastic scattering - call inelastic_scatter(nuc, rxn, p) - p % event_MT = rxn % MT + call inelastic_scatter(nuc, nuc%reactions(i), p) + p % event_MT = nuc%reactions(i)%MT end if @@ -1090,11 +1085,9 @@ contains real(8) :: weight ! weight adjustment for ufs method logical :: in_mesh ! source site in ufs mesh? type(Nuclide), pointer :: nuc - type(Reaction), pointer :: rxn ! Get pointers nuc => nuclides(i_nuclide) - rxn => nuc % reactions(i_reaction) ! TODO: Heat generation from fission @@ -1165,7 +1158,8 @@ contains ! Sample secondary energy distribution for fission reaction and set energy ! in fission bank - fission_bank(i) % E = sample_fission_energy(nuc, rxn, p) + fission_bank(i) % E = sample_fission_energy(nuc, nuc%reactions(& + i_reaction), p) ! Set the delayed group of the neutron fission_bank(i) % delayed_group = p % delayed_group diff --git a/src/tally.F90 b/src/tally.F90 index 2f178fe46..4790a9900 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -66,9 +66,6 @@ contains real(8) :: macro_scatt ! material macro scatt xs real(8) :: uvw(3) ! particle direction real(8) :: E ! particle energy - type(Material), pointer :: mat - type(Reaction), pointer :: rxn - type(Nuclide), pointer :: nuc i = 0 SCORE_LOOP: do q = 1, t % n_user_score_bins @@ -220,24 +217,20 @@ contains ! of one. score = p % last_wgt else - do m = 1, nuclides(p % event_nuclide) % n_reaction - ! Check if this is the desired MT - if (p % event_MT == nuclides(p % event_nuclide) % reactions(m) % MT) then - ! Found the reaction, set our pointer and move on with life - rxn => nuclides(p % event_nuclide) % reactions(m) - exit - end if - end do + m = nuclides(p%event_nuclide)%reaction_index% & + get_key(p % event_MT) ! Get multiplicity and apply to score - if (rxn % multiplicity_with_E) then - ! Then the multiplicity was already incorporated in to p % wgt - ! per the scattering routine, - score = p % wgt - else - ! Grab the multiplicity from the rxn - score = p % last_wgt * rxn % multiplicity - end if + associate (rxn => nuclides(p%event_nuclide)%reactions(m)) + if (rxn % multiplicity_with_E) then + ! Then the multiplicity was already incorporated in to p % wgt + ! per the scattering routine, + score = p % wgt + else + ! Grab the multiplicity from the rxn + score = p % last_wgt * rxn % multiplicity + end if + end associate end if @@ -257,24 +250,20 @@ contains ! of one. score = p % last_wgt else - do m = 1, nuclides(p % event_nuclide) % n_reaction - ! Check if this is the desired MT - if (p % event_MT == nuclides(p % event_nuclide) % reactions(m) % MT) then - ! Found the reaction, set our pointer and move on with life - rxn => nuclides(p % event_nuclide) % reactions(m) - exit - end if - end do + m = nuclides(p%event_nuclide)%reaction_index% & + get_key(p % event_MT) ! Get multiplicity and apply to score - if (rxn % multiplicity_with_E) then - ! Then the multiplicity was already incorporated in to p % wgt - ! per the scattering routine, - score = p % wgt - else - ! Grab the multiplicity from the rxn - score = p % last_wgt * rxn % multiplicity - end if + associate (rxn => nuclides(p%event_nuclide)%reactions(m)) + if (rxn % multiplicity_with_E) then + ! Then the multiplicity was already incorporated in to p % wgt + ! per the scattering routine, + score = p % wgt + else + ! Grab the multiplicity from the rxn + score = p % last_wgt * rxn % multiplicity + end if + end associate end if @@ -294,24 +283,20 @@ contains ! of one. score = p % last_wgt else - do m = 1, nuclides(p % event_nuclide) % n_reaction - ! Check if this is the desired MT - if (p % event_MT == nuclides(p % event_nuclide) % reactions(m) % MT) then - ! Found the reaction, set our pointer and move on with life - rxn => nuclides(p % event_nuclide) % reactions(m) - exit - end if - end do + m = nuclides(p%event_nuclide)%reaction_index% & + get_key(p % event_MT) ! Get multiplicity and apply to score - if (rxn % multiplicity_with_E) then - ! Then the multiplicity was already incorporated in to p % wgt - ! per the scattering routine, - score = p % wgt - else - ! Grab the multiplicity from the rxn - score = p % last_wgt * rxn % multiplicity - end if + associate (rxn => nuclides(p%event_nuclide)%reactions(m)) + if (rxn % multiplicity_with_E) then + ! Then the multiplicity was already incorporated in to p % wgt + ! per the scattering routine, + score = p % wgt + else + ! Grab the multiplicity from the rxn + score = p % last_wgt * rxn % multiplicity + end if + end associate end if @@ -466,9 +451,6 @@ contains ! delayed-nu-fission if (micro_xs(p % event_nuclide) % absorption > ZERO) then - ! Get the event nuclide - nuc => nuclides(p % event_nuclide) - ! Check if the delayed group filter is present if (dg_filter > 0) then @@ -480,11 +462,11 @@ contains d = t % filters(dg_filter) % int_bins(d_bin) ! Compute the yield for this delayed group - yield = yield_delayed(nuc, E, d) + yield = yield_delayed(nuclides(p % event_nuclide), E, d) ! Compute the score and tally to bin score = p % absorb_wgt * yield * micro_xs(p % event_nuclide) & - % fission * nu_delayed(nuc, E) / & + % fission * nu_delayed(nuclides(p % event_nuclide), E) / & micro_xs(p % event_nuclide) % absorption call score_fission_delayed_dg(t, d_bin, score, score_index) end do @@ -494,7 +476,7 @@ contains ! by multiplying the absorbed weight by the fraction of the ! delayed-nu-fission xs to the absorption xs score = p % absorb_wgt * micro_xs(p % event_nuclide) & - % fission * nu_delayed(nuc, E) / & + % fission * nu_delayed(nuclides(p % event_nuclide), E) / & micro_xs(p % event_nuclide) % absorption end if end if @@ -535,9 +517,6 @@ contains ! Check if tally is on a single nuclide if (i_nuclide > 0) then - ! Get the nuclide of interest - nuc => nuclides(i_nuclide) - ! Check if the delayed group filter is present if (dg_filter > 0) then @@ -548,11 +527,19 @@ contains d = t % filters(dg_filter) % int_bins(d_bin) ! Compute the yield for this delayed group +<<<<<<< HEAD yield = yield_delayed(nuc, E, d) ! Compute the score and tally to bin score = micro_xs(i_nuclide) % fission * yield & * nu_delayed(nuc, E) * atom_density * flux +======= + yield = yield_delayed(nuclides(i_nuclide), p % E, d) + + ! Compute the score and tally to bin + score = micro_xs(i_nuclide) % fission * yield & + * nu_delayed(nuclides(i_nuclide), p % E) * atom_density * flux +>>>>>>> Make Nuclide%reactions allocatable by using associate constructs call score_fission_delayed_dg(t, d_bin, score, score_index) end do cycle SCORE_LOOP @@ -560,27 +547,29 @@ contains ! If the delayed group filter is not present, compute the score ! by multiplying the delayed-nu-fission macro xs by the flux +<<<<<<< HEAD score = micro_xs(i_nuclide) % fission * nu_delayed(nuc, E)& * atom_density * flux +======= + score = micro_xs(i_nuclide) % fission * & + nu_delayed(nuclides(i_nuclide), p % E) * atom_density * flux +>>>>>>> Make Nuclide%reactions allocatable by using associate constructs end if ! Tally is on total nuclides else - ! Get pointer to current material - mat => materials(p % material) - ! Check if the delayed group filter is present if (dg_filter > 0) then ! Loop over all nuclides in the current material - do l = 1, mat % n_nuclides + do l = 1, materials(p % material) % n_nuclides ! Get atom density - atom_density_ = mat % atom_density(l) + atom_density_ = materials(p % material) % atom_density(l) ! Get index in nuclides array - i_nuc = mat % nuclide(l) + i_nuc = materials(p % material) % nuclide(l) ! Loop over all delayed group bins and tally to them individually do d_bin = 1, t % filters(dg_filter) % n_bins @@ -588,15 +577,20 @@ contains ! Get the delayed group for this bin d = t % filters(dg_filter) % int_bins(d_bin) - ! Get the current nuclide - nuc => nuclides(i_nuc) - ! Get the yield for the desired nuclide and delayed group +<<<<<<< HEAD yield = yield_delayed(nuc, E, d) ! Compute the score and tally to bin score = micro_xs(i_nuc) % fission * yield & * nu_delayed(nuc, E) * atom_density_ * flux +======= + yield = yield_delayed(nuclides(i_nuc), p % E, d) + + ! Compute the score and tally to bin + score = micro_xs(i_nuc) % fission * yield & + * nu_delayed(nuclides(i_nuc), p % E) * atom_density_ * flux +>>>>>>> Make Nuclide%reactions allocatable by using associate constructs call score_fission_delayed_dg(t, d_bin, score, score_index) end do end do @@ -606,13 +600,13 @@ contains score = ZERO ! Loop over all nuclides in the current material - do l = 1, mat % n_nuclides + do l = 1, materials(p % material) % n_nuclides ! Get atom density - atom_density_ = mat % atom_density(l) + atom_density_ = materials(p % material) % atom_density(l) ! Get index in nuclides array - i_nuc = mat % nuclide(l) + i_nuc = materials(p % material) % nuclide(l) ! Accumulate the contribution from each nuclide score = score + micro_xs(i_nuc) % fission & @@ -693,42 +687,41 @@ contains if (i_nuclide > 0) then if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then m = nuclides(i_nuclide)%reaction_index%get_key(score_bin) - rxn => nuclides(i_nuclide) % reactions(m) - - ! Retrieve index on nuclide energy grid and interpolation - ! factor - i_energy = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - if (i_energy >= rxn % threshold) then - score = ((ONE - f) * rxn % sigma(i_energy - & - rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density * flux - end if - end if - - else - ! Get pointer to current material - mat => materials(p % material) - do l = 1, mat % n_nuclides - ! Get atom density - atom_density_ = mat % atom_density(l) - - ! Get index in nuclides array - i_nuc = mat % nuclide(l) - - if (nuclides(i_nuc)%reaction_index%has_key(score_bin)) then - m = nuclides(i_nuc)%reaction_index%get_key(score_bin) - rxn => nuclides(i_nuc) % reactions(m) + associate (rxn => nuclides(i_nuclide) % reactions(m)) ! Retrieve index on nuclide energy grid and interpolation ! factor - i_energy = micro_xs(i_nuc) % index_grid - f = micro_xs(i_nuc) % interp_factor + i_energy = micro_xs(i_nuclide) % index_grid + f = micro_xs(i_nuclide) % interp_factor if (i_energy >= rxn % threshold) then - score = score + ((ONE - f) * rxn % sigma(i_energy - & + score = ((ONE - f) * rxn % sigma(i_energy - & rxn%threshold + 1) + f * rxn % sigma(i_energy - & - rxn%threshold + 2)) * atom_density_ * flux + rxn%threshold + 2)) * atom_density * flux end if + end associate + end if + + else + do l = 1, materials(p % material) % n_nuclides + ! Get atom density + atom_density_ = materials(p % material) % atom_density(l) + + ! Get index in nuclides array + i_nuc = materials(p % material) % nuclide(l) + + if (nuclides(i_nuc)%reaction_index%has_key(score_bin)) then + m = nuclides(i_nuc)%reaction_index%get_key(score_bin) + associate (rxn => nuclides(i_nuc) % reactions(m)) + ! Retrieve index on nuclide energy grid and interpolation + ! factor + i_energy = micro_xs(i_nuc) % index_grid + f = micro_xs(i_nuc) % interp_factor + if (i_energy >= rxn % threshold) then + score = score + ((ONE - f) * rxn % sigma(i_energy - & + rxn%threshold + 1) + f * rxn % sigma(i_energy - & + rxn%threshold + 2)) * atom_density_ * flux + end if + end associate end if end do end if