diff --git a/ChangeLog b/ChangeLog index 6c2b1d731..ccadd40f5 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,18 @@ +2011-03-25 Paul Romano + + * ace.f90: Added ability to read delayed neutron energy + distributions and precursor yields/constants. Created function + length_energy_dist from inner section of read_energy_dist so that + it could used for delayed neutrons as well. Added many comments + where they were previously missing. + * fileio.f90: Added many comments. + * global.f90: Added lattice types + * physics.f90: Added sampling of delayed neutrons (delayed/prompt, + precursor group, outgoing energy). Added many comments. + * types.f90: Moved all delayed neutron precursor data into one + array. Made delayed neutron energy distribution allocatable (since + there are 6 groups, each with a dist). Added many comments. + 2011-03-24 Paul Romano * geometry.f90: Fixed significant bug in how quadratic surfaces diff --git a/src/ace.f90 b/src/ace.f90 index 164b08c3b..a77164f82 100644 --- a/src/ace.f90 +++ b/src/ace.f90 @@ -9,10 +9,10 @@ module ace & dict_get_key, dict_delete use endf, only: reaction_name - integer :: NXS(16) - integer :: JXS(32) - real(8), allocatable :: XSS(:) - integer :: XSS_index + integer :: NXS(16) ! Descriptors for ACE XSS tables + integer :: JXS(32) ! Pointers into ACE XSS tables + real(8), allocatable :: XSS(:) ! Cross section data + integer :: XSS_index ! current index in XSS data private :: NXS private :: JXS @@ -27,18 +27,19 @@ contains subroutine read_xs() + integer :: i ! index in materials array + integer :: j ! index over isotopes in material + integer :: index ! index in xsdatas array + integer :: n ! length of key + integer :: index_continuous ! index in xs_continuous + integer :: index_thermal ! index in xs_thermal + character(10) :: key ! name of isotope, e.g. 92235.03c + character(250) :: msg ! output/error message type(Material), pointer :: mat => null() type(xsData), pointer :: iso => null() type(AceContinuous), pointer :: ace_cont => null() type(AceThermal), pointer :: ace_thermal => null() - integer :: i, j - integer :: index - character(10) :: key - character(250) :: msg - integer :: n - integer :: index_continuous - integer :: index_thermal - type(DictionaryCI), pointer :: temp_dict => null() + type(DictionaryCI), pointer :: temp_dict => null() call dict_create(ace_dict) @@ -114,24 +115,24 @@ contains subroutine read_ACE_continuous(index_table, index) - integer, intent(in) :: index_table - integer, intent(in) :: index + integer, intent(in) :: index_table ! index in xs_continuous array + integer, intent(in) :: index ! index in xsdatas array + integer :: in = 7 ! unit to read from + integer :: ioError ! error status for file access + integer :: words_per_line ! number of words per line (data) + integer :: lines ! number of lines (data + integer :: n ! number of data values + real(8) :: kT ! ACE table temperature + logical :: file_exists ! does ACE library exist? + logical :: found_xs ! did we find table in library? + character(7) :: readable ! is ACE library readable? + character(250) :: msg ! output/error message + character(250) :: line ! single line to read + character(32) :: words(max_words) ! words on a line + character(100) :: filename ! name of ACE library file + character(10) :: tablename ! name of cross section table type(AceContinuous), pointer :: table => null() - integer :: i - integer :: in = 7 - integer :: ioError - integer :: words_per_line - integer :: lines - integer :: n - logical :: file_exists - logical :: found_xs - character(7) :: readable - character(250) :: msg, line - character(32) :: words(max_words) - character(100) :: filename - character(10) :: tablename - real(8) :: kT filename = xsdatas(index)%path tablename = xsdatas(index)%id @@ -230,7 +231,8 @@ contains type(AceContinuous), pointer :: table - integer :: NE + integer :: NE ! number of energy points for total and elastic + ! cross sections ! determine number of energy points NE = NXS(3) @@ -246,7 +248,7 @@ contains ! read data from XSS -- right now the total, absorption and ! elastic scattering are read in to these special arrays, but in ! reality, it should be necessary to only store elastic scattering - ! and possible total cross-section for total material xs + ! and possibly total cross-section for total material xs ! generation. XSS_index = 1 table%energy = get_real(NE) @@ -267,30 +269,36 @@ contains type(AceContinuous), pointer :: table - integer :: JXS2 - integer :: KNU ! Location for nu data - integer :: LNU ! Type of nu data (polynomial or tabular) - integer :: NC ! Number of polynomial coefficients - integer :: NR ! Number of interpolation regions - integer :: NE ! Number of energies - integer :: length ! Length of data block to allocate + integer :: JXS2 ! location for fission nu data + integer :: JXS24 ! location for delayed neutron data + integer :: KNU ! location for nu data + integer :: LNU ! type of nu data (polynomial or tabular) + integer :: NC ! number of polynomial coefficients + integer :: NR ! number of interpolation regions + integer :: NE ! number of energies + integer :: NPCR ! number of delayed neutron precursor groups + 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 :: 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 :: JXS24 - integer :: DEC1 - integer :: DEC2 - - JXS2 = JXS(2) + JXS2 = JXS(2) JXS24 = JXS(24) if (JXS2 == 0) then ! ============================================================= - ! No prompt/total nu data is present + ! NO PROMPT/TOTAL NU DATA table % nu_t_type = NU_NONE table % nu_p_type = NU_NONE elseif (XSS(JXS2) > 0) then ! ============================================================= - ! Prompt or total nu data is present + ! PROMPT OR TOTAL NU DATA KNU = JXS2 LNU = int(XSS(KNU)) if (LNU == 1) then @@ -324,7 +332,7 @@ contains elseif (XSS(JXS2) < 0) then ! ============================================================= - ! Prompt and total nu data is present -- read prompt data first + ! PROMPT AND TOTAL NU DATA -- read prompt data first KNU = JXS2 + 1 LNU = XSS(KNU) if (LNU == 1) then @@ -383,7 +391,7 @@ contains if (JXS24 > 0) then ! ============================================================= - ! Delayed nu data is present + ! DELAYED NU DATA table % nu_d_type = NU_TABULAR KNU = JXS24 @@ -399,9 +407,85 @@ contains ! read delayed nu data XSS_index = KNU + 1 table % nu_d_data = get_real(length) + + ! ============================================================= + ! DELAYED NEUTRON ENERGY DISTRIBUTION + + ! Allocate space for secondary energy distribution + NPCR = NXS(8) + table % n_precursor = NPCR + allocate(table % nu_d_edist(NPCR)) + + LED = JXS(26) + LDIS = JXS(27) + + ! Loop over all delayed neutron precursor groups + do i = 1, NPCR + ! 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) + table % nu_d_edist(i) % law = LAW + + ! allocate space for ENDF interpolation parameters + if (NR > 0) then + allocate(table % nu_d_edist(i) % nbt(NR)) + allocate(table % nu_d_edist(i) % int(NR)) + end if + + ! read ENDF interpolation parameters + XSS_index = LDIS + LOCC + 3 + table % nu_d_edist(i) % nbt = get_real(NR) + table % nu_d_edist(i) % int = get_real(NR) + + ! allocate space for law validity data + NE = XSS(LDIS + LOCC + 3 + 2*NR) + allocate(table % nu_d_edist(i) % energy(NE)) + allocate(table % nu_d_edist(i) % pvalid(NE)) + + length_interp_data = 5 + 2*(NR + NE) + + ! read law validity data + XSS_index = LDIS + LOCC + 4 + 2*NR + table % nu_d_edist(i) % energy = get_real(NE) + table % nu_d_edist(i) % 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(table % nu_d_edist(i) % data(length)) + + ! read secondary energy distribution + XSS_index = loc + 1 + table % nu_d_edist(i) % data = get_real(length) + end do + + ! ============================================================= + ! DELAYED NEUTRON PRECUSOR YIELDS AND CONSTANTS + + ! determine length of all precursor constants/yields/interp data + length = 0 + loc = JXS(25) + do i = 1, NPCR + NR = XSS(loc + length + 1) + NE = XSS(loc + length + 2 + 2*NR) + length = length + 3 + 2*NR + 2*NE + end do + + ! allocate space for precusor data + allocate(table % nu_d_precursor_data(length)) + + ! read delayed neutron precursor data + XSS_index = loc + table % nu_d_precursor_data = get_real(length) - ! TODO: Read secondary energy distribution - ! TODO: Read precursor data else table % nu_d_type = NU_NONE end if @@ -418,22 +502,22 @@ contains type(AceContinuous), pointer :: table - type(AceReaction), pointer :: rxn => null() integer :: LMT ! index of MT list in XSS integer :: NMT ! Number of reactions integer :: JXS4 ! index of Q values in XSS integer :: JXS5 ! index of neutron multiplicities in XSS integer :: JXS7 ! index of reactions cross-sections in XSS - integer :: LXS ! - integer :: LOCA ! + integer :: LXS ! location of cross-section locators + integer :: LOCA ! location of cross-section for given MT integer :: NE ! number of energies for reaction + type(AceReaction), pointer :: rxn => null() - LMT = JXS(3) + LMT = JXS(3) JXS4 = JXS(4) JXS5 = JXS(5) LXS = JXS(6) JXS7 = JXS(7) - NMT = NXS(4) + NMT = NXS(4) ! allocate array of reactions. Add one since we need to include an ! elastic scattering channel @@ -474,22 +558,23 @@ contains !===================================================================== ! READ_ANGULAR_DIST parses the angular distribution for each reaction +! with secondary neutrons !===================================================================== subroutine read_angular_dist(table) type(AceContinuous), pointer :: table + integer :: JXS8 ! location of angular distribution locators + integer :: JXS9 ! location of angular distributions + integer :: LOCB ! location of angular distribution for given MT + integer :: NE ! number of incoming energies + integer :: NP ! number of points for cosine distribution + integer :: LC ! locator + integer :: i ! index in reactions array + integer :: j ! index over incoming energies + integer :: length ! length of data array to allocate type(AceReaction), pointer :: rxn => null() - integer :: JXS8 - integer :: JXS9 - integer :: NMT - integer :: LOCB - integer :: NE - integer :: NP - integer :: LC - integer :: i, j - integer :: length JXS8 = JXS(8) JXS9 = JXS(9) @@ -564,29 +649,29 @@ contains end subroutine read_angular_dist !===================================================================== -! READ_ENERGY_DIST +! READ_ENERGY_DIST parses the secondary energy distribution for each +! reaction with seconary neutrons (except elastic scattering) !===================================================================== subroutine read_energy_dist(table) type(AceContinuous), pointer :: table + 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, j, k, l ! indices type(AceReaction), pointer :: rxn => null() - integer :: LOCC - integer :: LED ! location of LDLW block - integer :: LDIS ! location of DLW block - integer :: LNW ! location of next law - integer :: LAW - integer :: NR, NE - integer :: NMU - integer :: NP - integer :: NP2 - integer :: NRa, NEa, NRb, NEb - integer :: IDAT - integer :: start, length, length_interp_data - integer :: i, j, k, l - LED = JXS(10) + LED = JXS(10) LDIS = JXS(11) ! Loop over all reactions @@ -627,135 +712,163 @@ contains rxn % edist % pvalid = get_real(NE) ! Set index to beginning of IDAT array - start = LDIS + IDAT - 2 + loc = LDIS + IDAT - 2 - ! Determine size of LDAT array based on which secondary energy - ! law it is - length = 0 - select case (LAW) - case (1) - ! Tabular equiprobable energy bins - NR = XSS(start + 1) - NE = XSS(start + 2 + 2*NR) - NET = XSS(start + 3 + 2*NR + NE) - length = 3 + 2*NR + NE + 3*NET*NE - - case (2) - ! Discrete photon energy - length = 2 - - case (3) - ! Level scattering - length = 2 - - case (4) - ! Continuous tabular distribution - NR = XSS(start + 1) - NE = XSS(start + 2 + 2*NR) - length = length + 2 + 2*NR + 2*NE - do j = 1,NE - ! determine length - NP = XSS(start + length + 2) - length = length + 2 + 3*NP - - ! adjust location for this block - k = start + 2 + 2*NR + NE + j - XSS(k) = XSS(k) - LOCC - length_interp_data - end do - k = start + 2 + 2*NR + NE - - case (5) - ! General evaporation spectrum - NR = XSS(start + 1) - NE = XSS(start + 2 + 2*NR) - NET = XSS(start + 3 + 2*NR + 2*NE) - length = 3 + 2*NR + 2*NE + NET - - case (7) - ! Maxwell fission spectrum - NR = XSS(start + 1) - NE = XSS(start + 2 + 2*NR) - length = 3 + 2*NR + 2*NE - - case (9) - ! Evaporation spectrum - NR = XSS(start + 1) - NE = XSS(start + 2 + 2*NR) - length = 3 + 2*NR + 2*NE - - case (11) - ! Watt spectrum - NRa = XSS(start + 1) - NEa = XSS(start + 2 + 2*NRa) - NRb = XSS(start + 3 + 2*(NRa+NEa)) - NEb = XSS(start + 4 + 2*(NRa+NEa+NRb)) - length = 5 + 2*(NRa + NEa + NRb + NEb) - - case (44) - ! Kalbach-Mann correlated scattering - NR = XSS(start + 1) - NE = XSS(start + 2 + 2*NR) - length = length + 2 + 2*NR + 2*NE - do j = 1,NE - NP = XSS(start + length + 2) - length = length + 2 + 5*NP - - ! adjust location for this block - k = start + 2 + 2*NR + NE + j - XSS(k) = XSS(k) - LOCC - length_interp_data - end do - - case (61) - ! Correlated energy and angle distribution - NR = XSS(start + 1) - NE = XSS(start + 2 + 2*NR) - length = length + 2 + 2*NR + 2*NE - do j = 1,NE - ! outgoing energy distribution - NP = XSS(start + length + 2) - - ! adjust locators for angular distribution - do k = 1, NP - l = start + length + 2 + 3*NP + k - if (XSS(l) /= 0) XSS(l) = XSS(l) - LOCC - length_interp_data - end do - - length = length + 2 + 4*NP - do k = 1, NP - ! outgoing angle distribution - NP2 = XSS(start + length + 2) - length = length + 2 + 3*NP2 - end do - - ! adjust locators for energy distribution - k = start + 2 + 2*NR + NE + j - XSS(k) = XSS(k) - LOCC - length_interp_data - end do - - case (66) - ! N-body phase space distribution - length = 2 - - case (67) - ! Laboratory energy-angle law - NR = XSS(start + 1) - NE = XSS(start + 2 + 2*NR) - NMU = XSS(start + 4 + 2*NR + 2*NE) - length = 4 + 2*(NR + NE + NMU) - - end select + ! 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 = start + 1 + XSS_index = loc + 1 rxn % edist % data = get_real(length) - end do end subroutine read_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 +!===================================================================== + + function length_energy_dist(loc, law, LOCC, lid) result(length) + + integer, intent(in) :: loc ! location in XSS array + integer, intent(in) :: law ! energy distribution law + integer, intent(in) :: LOCC ! location of energy distribution + integer, intent(in) :: lid ! length of interpolation data + integer :: length ! length of energy distribution (LDAT) + + integer :: i,j,k + integer :: NR ! number of interpolation regions + integer :: NE ! number of incoming energies + integer :: NP ! number of points in outgoing energy distribution + integer :: NMU ! number of points in outgoing cosine distribution + integer :: NRa ! number of interpolation regions for Watt 'a' + integer :: NEa ! number of energies for Watt 'a' + integer :: NRb ! number of interpolation regions for Watt 'b' + integer :: NEb ! number of energies for Watt 'b' + + ! initialize length + length = 0 + + select case (law) + case (1) + ! Tabular equiprobable energy bins + NR = XSS(loc + 1) + NE = XSS(loc + 2 + 2*NR) + NP = XSS(loc + 3 + 2*NR + NE) + length = 3 + 2*NR + NE + 3*NP*NE + + case (2) + ! Discrete photon energy + length = 2 + + case (3) + ! Level scattering + length = 2 + + case (4) + ! Continuous tabular distribution + NR = XSS(loc + 1) + NE = XSS(loc + 2 + 2*NR) + length = length + 2 + 2*NR + 2*NE + do i = 1,NE + ! determine length + NP = XSS(loc + length + 2) + length = length + 2 + 3*NP + + ! adjust location for this block + j = loc + 2 + 2*NR + NE + i + XSS(j) = XSS(j) - LOCC - lid + end do + + case (5) + ! General evaporation spectrum + NR = XSS(loc + 1) + NE = XSS(loc + 2 + 2*NR) + NP = XSS(loc + 3 + 2*NR + 2*NE) + length = 3 + 2*NR + 2*NE + NP + + case (7) + ! Maxwell fission spectrum + NR = XSS(loc + 1) + NE = XSS(loc + 2 + 2*NR) + length = 3 + 2*NR + 2*NE + + case (9) + ! Evaporation spectrum + NR = XSS(loc + 1) + NE = XSS(loc + 2 + 2*NR) + length = 3 + 2*NR + 2*NE + + case (11) + ! Watt spectrum + NRa = XSS(loc + 1) + NEa = XSS(loc + 2 + 2*NRa) + NRb = XSS(loc + 3 + 2*(NRa+NEa)) + NEb = XSS(loc + 4 + 2*(NRa+NEa+NRb)) + length = 5 + 2*(NRa + NEa + NRb + NEb) + + case (44) + ! Kalbach-Mann correlated scattering + NR = XSS(loc + 1) + NE = XSS(loc + 2 + 2*NR) + length = length + 2 + 2*NR + 2*NE + do i = 1,NE + NP = XSS(loc + length + 2) + length = length + 2 + 5*NP + + ! adjust location for this block + j = loc + 2 + 2*NR + NE + i + XSS(j) = XSS(j) - LOCC - lid + end do + + case (61) + ! Correlated energy and angle distribution + NR = XSS(loc + 1) + NE = XSS(loc + 2 + 2*NR) + length = length + 2 + 2*NR + 2*NE + do i = 1,NE + ! outgoing energy distribution + NP = XSS(loc + length + 2) + + ! adjust locators for angular distribution + do j = 1, NP + k = loc + length + 2 + 3*NP + j + if (XSS(k) /= 0) XSS(k) = XSS(k) - LOCC - lid + end do + + length = length + 2 + 4*NP + do j = 1, NP + ! outgoing angle distribution -- NMU here is actually + ! referred to as NP in the MCNP documentation + NMU = XSS(loc + length + 2) + length = length + 2 + 3*NMU + end do + + ! adjust locators for energy distribution + j = loc + 2 + 2*NR + NE + i + XSS(j) = XSS(j) - LOCC - lid + end do + + case (66) + ! N-body phase space distribution + length = 2 + + case (67) + ! Laboratory energy-angle law + NR = XSS(loc + 1) + NE = XSS(loc + 2 + 2*NR) + NMU = XSS(loc + 4 + 2*NR + 2*NE) + length = 4 + 2*(NR + NE + NMU) + + end select + + end function length_energy_dist + !===================================================================== ! GET_INT returns an array of integers read from the current position ! in the XSS array @@ -763,8 +876,8 @@ contains function get_int(n_values) result(array) - integer, intent(in) :: n_values - integer :: array(n_values) + integer, intent(in) :: n_values ! number of values to read + integer :: array(n_values) ! array of values array = int(XSS(XSS_index:XSS_index + n_values - 1)) XSS_index = XSS_index + n_values @@ -778,8 +891,8 @@ contains function get_real(n_values) result(array) - integer, intent(in) :: n_values - real(8) :: array(n_values) + integer, intent(in) :: n_values ! number of values to read + real(8) :: array(n_values) ! array of values array = XSS(XSS_index:XSS_index + n_values - 1) XSS_index = XSS_index + n_values diff --git a/src/endf.f90 b/src/endf.f90 index f232a39e2..9321f20cb 100644 --- a/src/endf.f90 +++ b/src/endf.f90 @@ -18,6 +18,8 @@ contains string = '(n,total)' case (2) string = '(elastic)' + case (4) + string = '(n,level)' case (11) string = '(n,2nd)' case (16) diff --git a/src/fileio.f90 b/src/fileio.f90 index 0b1f51054..ba4b20b1e 100644 --- a/src/fileio.f90 +++ b/src/fileio.f90 @@ -36,13 +36,13 @@ contains subroutine read_command_line() - integer :: argc - character(250) :: msg - logical :: file_exists - character(7) :: readable + integer :: argc ! number of command line arguments + logical :: file_exists ! does specified input file exist? + character(250) :: msg ! error message + character(7) :: readable ! is input file readable? argc = COMMAND_ARGUMENT_COUNT() - if ( argc > 0 ) then + if (argc > 0) then call GET_COMMAND_ARGUMENT(1, path_input) else msg = "No input file specified!" @@ -75,16 +75,17 @@ contains character(*), intent(in) :: filename - character(250) :: line, msg - character(32) :: words(max_words) - integer :: in = 7 - integer :: ioError - integer :: n - integer :: count - integer :: universe_num - type(ListKeyValueII), pointer :: key_list - type(Universe), pointer :: univ - integer :: index + integer :: in = 7 ! unit # for input file + integer :: index ! index in universes array + integer :: ioError ! error status for file access + integer :: n ! number of words on a line + integer :: count ! number of cells in a universe + integer :: universe_num ! user-specified universe # + character(250) :: line ! a line of words in input file + character(250) :: msg ! output/error message + character(32) :: words(max_words) ! words on a line + type(ListKeyValueII), pointer :: key_list => null() + type(Universe), pointer :: univ => null() msg = "First pass through input file..." call message(msg, 5) @@ -206,17 +207,18 @@ contains character(*), intent(in) :: filename - character(250) :: line, msg - character(32) :: words(max_words) - integer :: in = 7 - integer :: ioError - integer :: n, i - integer :: index_cell - integer :: index_universe - integer :: index_surface - integer :: index_lattice - integer :: index_material - integer :: index_source + integer :: in = 7 ! unit # for input file + integer :: ioError ! error status for file access + integer :: n ! number of words on a line + integer :: index_cell ! index in cells array + integer :: index_universe ! index in universes array + integer :: index_surface ! index in surfaces array + integer :: index_lattice ! index in lattices array + integer :: index_material ! index in materials array + integer :: index_source ! index in source array (?) + character(250) :: line ! a line of words + character(250) :: msg ! output/error message + character(32) :: words(max_words) ! words on a single line msg = "Second pass through input file..." call message(msg, 5) @@ -237,7 +239,7 @@ contains do call get_next_line(in, words, n, ioError) - if ( ioError /= 0 ) exit + if (ioError /= 0) exit ! skip comments if (n==0 .or. words(1)(1:1) == '#') cycle @@ -294,7 +296,7 @@ contains end subroutine read_input !===================================================================== -! ADJUST_INDICES changes the boundary_list values for each cell and +! ADJUST_INDICES changes the values for 'surfaces' for each cell and ! the material index assigned to each to the indices in the surfaces ! and material array rather than the unique IDs assigned to each ! surface and material @@ -302,12 +304,12 @@ contains subroutine adjust_indices() - type(Cell), pointer :: c => null() - - integer :: i, j - integer :: index - integer :: surf_num - character(250) :: msg + integer :: i ! index in cells array + integer :: j ! index over surface list + integer :: index ! index in surfaces/materials array + integer :: surf_num ! user-specified surface number + character(250) :: msg ! output/error message + type(Cell), pointer :: c => null() ! cell pointer do i = 1, n_cells ! adjust boundary list @@ -327,9 +329,7 @@ contains end do ! adjust material indices - if (c % fill /= 0) then - ! cell is filled with universe -- do nothing - else + if (c % fill == 0) then index = dict_get_key(material_dict, c%material) if (index == DICT_NULL) then msg = "Could not find material " // trim(int_to_str(c%material)) // & @@ -349,12 +349,12 @@ contains recursive subroutine build_universe(univ, parent, level) - type(Universe), pointer :: univ - integer, intent(in) :: parent - integer, intent(in) :: level + type(Universe), pointer :: univ ! univese pointer + integer, intent(in) :: parent ! cell containing universe + integer, intent(in) :: level ! level of universe - integer :: i - integer :: index + integer :: i ! index for cells in universe + integer :: index ! index in cells array type(Cell), pointer :: c => null() type(Universe), pointer :: subuniverse => null() @@ -367,8 +367,8 @@ contains c => cells(index) c%parent = parent - ! if this cell is filled with another universe, recursive call - ! this subroutine + ! if this cell is filled with another universe, recursively + ! call this subroutine if (c % fill > 0) then subuniverse => universes(c % fill) call build_universe(subuniverse, index, level + 1) @@ -383,16 +383,16 @@ contains subroutine read_cell(index, words, n_words) - integer, intent(in) :: index - character(*), intent(in) :: words(max_words) - integer, intent(in) :: n_words + integer, intent(in) :: index ! index in cells array + character(*), intent(in) :: words(n_words) ! words on cell card entry + integer, intent(in) :: n_words ! number of words - integer :: ioError - integer :: i - integer :: universe_num - integer :: n_surfaces - character(250) :: msg - character(32) :: word + integer :: ioError ! error status for file access + integer :: i ! index for surface list in a cell + integer :: universe_num ! user-specified universe number + integer :: n_surfaces ! number of surfaces in a cell + character(250) :: msg ! output/error message + character(32) :: word ! single word type(Cell), pointer :: c => null() c => cells(index) @@ -417,7 +417,7 @@ contains if (trim(words(4)) == 'fill') then c % type = CELL_FILL c % material = 0 - n_surfaces = n_words - 5 + n_surfaces = n_words - 5 ! find universe universe_num = str_to_int(words(5)) @@ -471,69 +471,69 @@ contains subroutine read_surface(index, words, n_words) - integer, intent(in) :: index - character(*), intent(in) :: words(max_words) - integer, intent(in) :: n_words + integer, intent(in) :: index ! index in surfaces array + character(*), intent(in) :: words(n_words) ! words in surface card entry + integer, intent(in) :: n_words ! number of words - integer :: ioError - integer :: i - integer :: coeffs_reqd - character(250) :: msg - character(32) :: word - type(Surface), pointer :: this_surface => null() + integer :: ioError ! error status for file access + integer :: i ! index for surface coefficients + integer :: coeffs_reqd ! number of coefficients are required + character(250) :: msg ! output/error message + character(32) :: word ! single word + type(Surface), pointer :: surf => null() - this_surface => surfaces(index) + surf => surfaces(index) ! Read surface identifier - read(words(2), FMT='(I8)', IOSTAT=ioError) this_surface%uid + read(words(2), FMT='(I8)', IOSTAT=ioError) surf % uid if (ioError > 0) then msg = "Invalid surface name: " // words(2) call error(msg) end if - call dict_add_key(surface_dict, this_surface%uid, index) + call dict_add_key(surface_dict, surf % uid, index) ! Read surface type word = words(3) call lower_case(word) select case (trim(word)) case ('px') - this_surface%type = SURF_PX + surf % type = SURF_PX coeffs_reqd = 1 case ('py') - this_surface%type = SURF_PY + surf % type = SURF_PY coeffs_reqd = 1 case ('pz') - this_surface%type = SURF_PZ + surf % type = SURF_PZ coeffs_reqd = 1 case ('plane') - this_surface%type = SURF_PLANE + surf % type = SURF_PLANE coeffs_reqd = 4 case ('cylx') - this_surface%type = SURF_CYL_X + surf % type = SURF_CYL_X coeffs_reqd = 3 case ('cyly') - this_surface%type = SURF_CYL_Y + surf % type = SURF_CYL_Y coeffs_reqd = 3 case ('cylz') - this_surface%type = SURF_CYL_Z + surf % type = SURF_CYL_Z coeffs_reqd = 3 case ('sph') - this_surface%type = SURF_SPHERE + surf % type = SURF_SPHERE coeffs_reqd = 4 case ('boxx') - this_surface%type = SURF_BOX_X + surf % type = SURF_BOX_X coeffs_reqd = 4 case ('boxy') - this_surface%type = SURF_BOX_Y + surf % type = SURF_BOX_Y coeffs_reqd = 4 case ('boxz') - this_surface%type = SURF_BOX_Z + surf % type = SURF_BOX_Z coeffs_reqd = 4 case ('box') - this_surface%type = SURF_BOX + surf % type = SURF_BOX coeffs_reqd = 6 case ('gq') - this_surface%type = SURF_GQ + surf % type = SURF_GQ coeffs_reqd = 10 case default msg = "Invalid surface type: " // words(3) @@ -547,9 +547,9 @@ contains end if ! Read list of surfaces - allocate(this_surface%coeffs(n_words-3)) + allocate(surf%coeffs(n_words-3)) do i = 1, n_words-3 - this_surface%coeffs(i) = str_to_real(words(i+3)) + surf % coeffs(i) = str_to_real(words(i+3)) end do end subroutine read_surface @@ -560,9 +560,9 @@ contains subroutine read_lattice(index, words, n_words) - integer, intent(in) :: index - character(*), intent(in) :: words(max_words) - integer, intent(in) :: n_words + integer, intent(in) :: index + character(*), intent(in) :: words(n_words) + integer, intent(in) :: n_words end subroutine read_lattice @@ -572,21 +572,21 @@ contains subroutine read_source(words, n_words) - character(*), intent(in) :: words(max_words) - integer, intent(in) :: n_words + character(*), intent(in) :: words(n_words) ! words on source entry + integer, intent(in) :: n_words ! number of words - character(250) :: msg - character(32) :: word - integer :: ioError - integer :: values_reqd - integer :: i + integer :: i ! index in values list + integer :: ioError ! error status for file access + integer :: values_reqd ! number of values required to specify source + character(250) :: msg ! output/error message + character(32) :: word ! single word ! Read source type word = words(2) call lower_case(word) select case (trim(word)) case ('box') - external_source%type = SRC_BOX + external_source % type = SRC_BOX values_reqd = 6 case default msg = "Invalid source type: " // words(2) @@ -599,10 +599,10 @@ contains call error(msg) end if - ! Read list of surfaces + ! Read values allocate(external_source%values(n_words-2)) do i = 1, n_words-2 - external_source%values(i) = str_to_real(words(i+2)) + external_source % values(i) = str_to_real(words(i+2)) end do end subroutine read_source @@ -614,27 +614,29 @@ contains subroutine read_material(index, words, n_words) - integer, intent(in) :: index - character(*), intent(in) :: words(n_words) - integer, intent(in) :: n_words + integer, intent(in) :: index ! index in materials array + character(*), intent(in) :: words(n_words) ! words on material entry + integer, intent(in) :: n_words ! number of words - character(100) :: msg - integer :: i - integer :: ioError - integer :: n_isotopes + integer :: i ! index over isotopes + integer :: ioError ! error status for file access + integer :: n_isotopes ! number of isotopes in material + character(250) :: msg ! output/error message type(Material), pointer :: mat => null() + ! Check for correct number of arguments if (mod(n_words,2) == 0 .or. n_words < 5) then msg = "Invalid number of arguments for material: " // words(2) call error(msg) end if + ! Determine and set number of isotopes n_isotopes = (n_words-3)/2 mat => materials(index) - mat%n_isotopes = n_isotopes + mat % n_isotopes = n_isotopes ! Read surface identifier - read(words(2), FMT='(I8)', IOSTAT=ioError) mat%uid + read(words(2), FMT='(I8)', IOSTAT=ioError) mat % uid if (ioError > 0) then msg = "Invalid surface name: " // words(2) call error(msg) @@ -642,7 +644,7 @@ contains call dict_add_key(material_dict, mat%uid, index) ! Read atom density - mat%atom_density = str_to_real(words(3)) + mat % atom_density = str_to_real(words(3)) ! allocate isotope and density list allocate(mat%names(n_isotopes)) @@ -652,8 +654,8 @@ contains ! read names and percentage do i = 1, n_isotopes - mat%names(i) = words(2*i+2) - mat%atom_percent(i) = str_to_real(words(2*i+3)) + mat % names(i) = words(2*i+2) + mat % atom_percent(i) = str_to_real(words(2*i+3)) end do end subroutine read_material @@ -665,18 +667,19 @@ contains subroutine normalize_ao() + integer :: index ! index used for several purposes + integer :: i ! index in materials array + integer :: j ! index over isotopes in material + real(8) :: sum_percent ! + real(8) :: awr ! atomic weight ratio + real(8) :: w ! weight percent + real(8) :: x ! atom percent + logical :: percent_in_atom ! isotopes specified in atom percent? + logical :: density_in_atom ! density specified in atom/b-cm? + character(10) :: key ! name of isotopes, e.g. 92235.03c + character(100) :: msg ! output/error message type(xsData), pointer :: iso => null() type(Material), pointer :: mat => null() - integer :: index - integer :: i, j - real(8) :: sum_percent - real(8) :: awr ! atomic weight ratio - real(8) :: w ! weight percent - real(8) :: x ! atom percent - logical :: percent_in_atom - logical :: density_in_atom - character(10) :: key - character(100) :: msg ! first find the index in the xsdata array for each isotope in ! each material @@ -685,31 +688,31 @@ contains ! Check to make sure either all atom percents or all weight ! percents are given - if (.not. (all(mat%atom_percent > 0) .or. & - & all(mat%atom_percent < 0))) then + if (.not. (all(mat%atom_percent > ZERO) .or. & + & all(mat%atom_percent < ZERO))) then msg = "Cannot mix atom and weight percents in material " // & & int_to_str(mat%uid) call error(msg) end if - percent_in_atom = (mat%atom_percent(1) > 0) - density_in_atom = (mat%atom_density > 0) + percent_in_atom = (mat%atom_percent(1) > ZERO) + density_in_atom = (mat%atom_density > ZERO) - sum_percent = 0 - do j = 1, mat%n_isotopes + sum_percent = ZERO + do j = 1, mat % n_isotopes ! Set indices for isotopes - key = mat%names(j) + key = mat % names(j) index = dict_get_key(xsdata_dict, key) - mat%isotopes(j) = index + mat % isotopes(j) = index ! determine atomic weight ratio - awr = xsdatas(index)%awr + awr = xsdatas(index) % awr ! if given weight percent, convert all values so that they ! are divided by awr. thus, when a sum is done over the ! values, it's actually sum(w/awr) if (.not. percent_in_atom) then - mat%atom_percent(j) = -mat%atom_percent(j) / awr + mat % atom_percent(j) = -mat % atom_percent(j) / awr end if end do @@ -717,20 +720,20 @@ contains ! this is straightforward. if given weight percents, the value ! is w/awr and is divided by sum(w/awr) sum_percent = sum(mat%atom_percent) - mat%atom_percent = mat%atom_percent / sum_percent + mat % atom_percent = mat % atom_percent / sum_percent ! Change density in g/cm^3 to atom/b-cm. Since all values are ! now in atom percent, the sum needs to be re-evaluated as ! 1/sum(x*awr) if (.not. density_in_atom) then - sum_percent = 0 - do j = 1, mat%n_isotopes - index = mat%isotopes(j) - awr = xsdatas(index)%awr - x = mat%atom_percent(j) + sum_percent = ZERO + do j = 1, mat % n_isotopes + index = mat % isotopes(j) + awr = xsdatas(index) % awr + x = mat % atom_percent(j) sum_percent = sum_percent + x*awr end do - sum_percent = 1.0_8 / sum_percent + sum_percent = ONE / sum_percent mat%atom_density = -mat%atom_density * N_AVOGADRO & & / MASS_NEUTRON * sum_percent end if @@ -743,10 +746,10 @@ contains ! specifies what cross section library should be used by default !===================================================================== - subroutine read_xs_library( words, n_words ) + subroutine read_xs_library(words, n_words) - character(*), intent(in) :: words(max_words) - integer, intent(in) :: n_words + character(*), intent(in) :: words(n_words) + integer, intent(in) :: n_words end subroutine read_xs_library @@ -757,34 +760,33 @@ contains ! cycle. !===================================================================== - subroutine read_criticality( words, n_words ) + subroutine read_criticality(words, n_words) - character(*), intent(in) :: words(max_words) - integer, intent(in) :: n_words + character(*), intent(in) :: words(n_words) ! words on criticality card + integer, intent(in) :: n_words ! number of words - character(250) :: msg - integer :: ioError + character(250) :: msg ! output/error message ! Set problem type to criticality problem_type = PROB_CRITICALITY ! Read number of cycles - read(words(2), FMT='(I8)', IOSTAT=ioError) n_cycles - if (ioError > 0) then + n_cycles = str_to_int(words(2)) + if (n_cycles == ERROR_CODE) then msg = "Invalid number of cycles: " // words(2) call error(msg) end if ! Read number of inactive cycles - read(words(3), FMT='(I8)', IOSTAT=ioError) n_inactive - if (ioError > 0) then + n_inactive = str_to_int(words(3)) + if (n_inactive == ERROR_CODE) then msg = "Invalid number of inactive cycles: " // words(2) call error(msg) end if ! Read number of particles - read(words(4), FMT='(I8)', IOSTAT=ioError) n_particles - if (ioError > 0) then + n_particles = str_to_int(words(4)) + if (n_particles == ERROR_CODE) then msg = "Invalid number of particles: " // words(2) call error(msg) end if @@ -797,9 +799,9 @@ contains subroutine read_line(unit, line, ioError) - integer, intent(in) :: unit - character(max_line), intent(out) :: line - integer, intent(out) :: ioError + integer, intent(in) :: unit ! unit to read from + character(max_line), intent(out) :: line ! line to return + integer, intent(out) :: ioError ! error status read(UNIT=unit, FMT='(A250)', IOSTAT=ioError) line @@ -814,30 +816,41 @@ contains subroutine get_next_line(unit, words, n, ioError) - integer, intent(in) :: unit - character(*), intent(out) :: words(max_words) - integer, intent(out) :: n - integer, intent(out) :: ioError + integer, intent(in) :: unit ! unit to read from + character(*), intent(out) :: words(max_words) ! words read + integer, intent(out) :: n ! number of words + integer, intent(out) :: ioError ! error status - character(250) :: line - character(32) :: local_words(max_words) - integer :: index + character(250) :: line ! single line + character(32) :: local_words(max_words) ! words on one line + integer :: index ! index of words index = 0 - do + do + ! read line from file read(UNIT=unit, FMT='(A250)', IOSTAT=ioError) line + + ! if we're at the end of the file, return if (ioError /= 0) return + + ! split a single line into words call split_string_wl(line, local_words, n) + + ! if there are no words, we're done if (n == 0) exit + + ! Check whether there is a continuation line if (local_words(n) == '&') then words(index+1:index+n-1) = local_words(1:n-1) - index = index + n-1 + index = index + n - 1 else words(index+1:index+n) = local_words(1:n) index = index + n exit end if end do + + ! set total number of words n = index end subroutine get_next_line @@ -848,12 +861,12 @@ contains subroutine skip_lines(unit, n_lines, ioError) - integer, intent(in) :: unit - integer, intent(in) :: n_lines - integer, intent(out) :: ioError + integer, intent(in) :: unit ! unit to read from + integer, intent(in) :: n_lines ! number of lines to skip + integer, intent(out) :: ioError ! error status - integer :: i - character(max_line) :: tmp + integer :: i ! index for number of lines + character(max_line) :: tmp ! single line do i = 1, n_lines read(UNIT=unit, FMT='(A250)', IOSTAT=ioError) tmp @@ -862,63 +875,56 @@ contains end subroutine skip_lines !===================================================================== -! READ_DATA_INT reads integer data into an arrya from a file open +! READ_DATA_INT reads integer data into an array from a file open !===================================================================== - subroutine read_data_int(in, array, n, lines, words_per_line) + subroutine read_data_int(unit, array, n, lines, words_per_line) - integer, intent(in) :: in - integer, intent(out) :: array(n) - integer, intent(in) :: n - integer, intent(in) :: lines - integer, intent(in) :: words_per_line + integer, intent(in) :: unit ! unit to read from + integer, intent(out) :: array(n) ! ints read from file + integer, intent(in) :: n ! total number of ints + integer, intent(in) :: lines ! total number of lines + integer, intent(in) :: words_per_line ! number of words per line - integer :: i, j - character(250) :: line - character(32) :: words(max_words) - integer :: ioError - integer :: n_words - integer :: index - integer :: n_last + integer :: i ! line index + integer :: loc ! locator for array + integer :: ioError ! error status - index = 1 + loc = 0 do i = 1, lines if (i == lines) then - read(in,*) array(index:n) + read(UNIT=unit,FMT=*) array(loc+1:n) else - read(in,*) array(index:index+words_per_line-1) - index = index + words_per_line + read(UNIT=unit,FMT=*) array(loc+1:loc+words_per_line) + loc = loc + words_per_line end if end do end subroutine read_data_int !===================================================================== -! READ_DATA_REAL reads real(8) data into an arrya from a file open +! READ_DATA_REAL reads real(8) data into an array from a file open !===================================================================== - subroutine read_data_real(in, array, n, lines, words_per_line) + subroutine read_data_real(unit, array, n, lines, words_per_line) - integer, intent(in) :: in - real(8), intent(out) :: array(n) - integer, intent(in) :: n - integer, intent(in) :: lines - integer, intent(in) :: words_per_line + integer, intent(in) :: unit ! unit to read from + real(8), intent(out) :: array(n) ! real(8)s read from file + integer, intent(in) :: n ! total number of ints + integer, intent(in) :: lines ! total number of lines + integer, intent(in) :: words_per_line ! number of words per line - integer :: i, j - character(250) :: line - character(32) :: words(max_words) - integer :: ioError - integer :: n_words - integer :: index + integer :: i ! line index + integer :: loc ! locator for array + integer :: ioError ! error status - index = 1 + loc = 0 do i = 1, lines if (i == lines) then - read(in,*) array(index:n) + read(UNIT=unit,FMT=*) array(loc+1:n) else - read(in,*) array(index:index+words_per_line-1) - index = index + words_per_line + read(UNIT=unit,FMT=*) array(loc+1:loc+words_per_line) + loc = loc + words_per_line end if end do diff --git a/src/global.f90 b/src/global.f90 index f68feb91c..7db0e7dcb 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -107,6 +107,11 @@ module global & CELL_FILL = 2, & ! Cell filled by a separate universe & CELL_LATTICE = 3 ! Cell filled with a lattice + ! Lattice types + integer, parameter :: & + & LATTICE_RECT = 1, & + & LATTICE_HEX = 2 + ! array index of universe 0 integer :: BASE_UNIVERSE @@ -247,13 +252,13 @@ contains integer, intent(in) :: num character(10) :: str - write ( str, '(I10)' ) num + write (str, '(I10)') num str = adjustl(str) end function int_to_str !===================================================================== -! INT_TO_STR converts a string to an integer. +! STR_TO_INT converts a string to an integer. !===================================================================== function str_to_int(str) result(num) diff --git a/src/physics.f90 b/src/physics.f90 index 3c4c810ff..a39252839 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -21,18 +21,15 @@ contains type(Particle), pointer :: p - character(250) :: msg - integer :: i - integer :: surf - logical :: alive - logical :: found_cell - - real(8) :: d_to_boundary - real(8) :: d_to_collision - real(8) :: distance - real(8) :: Sigma ! total cross-section - real(8) :: f ! interpolation factor - integer :: IE ! index on energy grid + integer :: surf ! surface which particle is on + integer :: IE ! index on energy grid + real(8) :: d_to_boundary ! distance to nearest boundary + real(8) :: d_to_collision ! sampled distance to collision + real(8) :: distance ! distance particle travels + real(8) :: Sigma ! total cross-section + real(8) :: f ! interpolation factor + logical :: found_cell ! found cell which particle is in? + character(250) :: msg ! output/error message type(Universe), pointer :: univ if (p % cell == 0) then @@ -48,13 +45,12 @@ contains end if if (verbosity >= 9) then - msg = "=== Particle " // trim(int_to_str(p % uid)) // " ===" + msg = "Simulating Particle " // trim(int_to_str(p % uid)) call message(msg, 9) end if if (verbosity >= 10) then - i = cells(p % cell)%uid - msg = " Born in cell " // trim(int_to_str(i)) + msg = " Born in cell " // trim(int_to_str(cells(p%cell)%uid)) call message(msg, 10) end if @@ -107,9 +103,9 @@ contains type(Particle), pointer :: p - integer :: IE - real(8) :: E - real(8) :: interp + integer :: IE ! index on union energy grid + real(8) :: E ! energy of particle + real(8) :: interp ! interpolation factor ! copy particle's energy E = p % E @@ -137,26 +133,30 @@ contains end subroutine find_energy_index !===================================================================== -! COLLISION +! COLLISION samples a nuclide and reaction and then calls the +! appropriate routine for that reaction !===================================================================== subroutine collision(p) type(Particle), pointer :: p + integer :: i ! index over isotopes in a material + integer :: n_isotopes ! number of isotopes in material + integer :: IE ! index on nuclide energy grid + real(8) :: f ! interpolation factor + real(8) :: sigma ! microscopic total xs for nuclide + real(8) :: total ! total macroscopic xs for material + real(8) :: density ! atom density of material + real(8) :: density_i ! atom density of nuclide + real(8) :: prob ! cumulative probability + real(8) :: r1 ! random number + real(8), allocatable :: Sigma_t(:) ! macroscopic xs for each nuclide + character(250) :: msg ! output/error message type(AceContinuous), pointer :: table type(AceReaction), pointer :: rxn - real(8) :: r1 - character(250) :: msg - integer :: i,j - integer :: n_isotopes - integer :: IE - real(8) :: f, Sigma, total - real(8) :: density, density_i - real(8) :: prob - real(8), allocatable :: Sigma_t(:) - ! tallies + ! TODO: tallies density = cMaterial%atom_density @@ -172,14 +172,15 @@ contains IE = table%grid_index(p % IE) f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE)) - Sigma = density_i*((ONE-f)*table%Sigma_t(IE) + f*(table%Sigma_t(IE+1))) - Sigma_t(i) = Sigma + ! calculate macroscopic xs for this nuclide + Sigma_t(i) = density_i*((ONE-f)*table%Sigma_t(IE) + & + & f*(table%Sigma_t(IE+1))) end do + total = sum(Sigma_t) ! sample nuclide r1 = rang() prob = ZERO - total = sum(Sigma_t) do i = 1, n_isotopes prob = prob + Sigma_t(i) / total if (r1 < prob) exit @@ -188,7 +189,7 @@ contains ! Get table, total xs, interpolation factor density_i = cMaterial%atom_percent(i)*density table => xs_continuous(cMaterial%table(i)) - Sigma = Sigma_t(i) / density_i + sigma = Sigma_t(i) / density_i IE = table%grid_index(p % IE) f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE)) @@ -196,16 +197,25 @@ contains deallocate(Sigma_t) ! sample reaction channel - r1 = rang()*Sigma + r1 = rang()*sigma prob = ZERO do i = 1, table%n_reaction rxn => table%reactions(i) + + ! 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 == 4) cycle + + ! if energy is below threshold for this reaction, skip it if (IE < rxn%IE) cycle + + ! add to cumulative probability prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) & & + f*(rxn%sigma(IE-rxn%IE+2))) if (r1 < prob) exit end do + if (verbosity >= 10 .or. p % uid == 4593) then msg = " " // trim(reaction_name(rxn%MT)) // " with nuclide " // & & trim(table%name) @@ -213,7 +223,7 @@ contains end if ! call appropriate subroutine - select case (rxn%MT) + select case (rxn % MT) case (2) call elastic_scatter(p, table, rxn) case (11, 16, 17, 24, 25, 30, 37, 41, 42) @@ -225,7 +235,8 @@ contains case (102:117) call n_absorption(p) case default - print *, 'warning: cant simulate reaction: ', rxn % MT + msg = "Cannot simulate reaction with MT " // int_to_str(rxn%MT) + call warning(msg) end select ! check for very low energy @@ -241,7 +252,9 @@ contains end subroutine collision !===================================================================== -! ELASTIC_SCATTER +! ELASTIC_SCATTER treats the elastic scattering of a neutron with a +! target. Currently this assumes target-at-rest kinematics -- +! obviously will need to be fixed !===================================================================== subroutine elastic_scatter(p, table, rxn) @@ -250,22 +263,27 @@ contains type(AceContinuous), pointer :: table type(AceReaction), pointer :: rxn - real(8) :: awr + real(8) :: awr ! atomic weight ratio of target real(8) :: phi ! azimuthal angle real(8) :: mu ! cosine of polar angle - real(8) :: vx, vy, vz - real(8) :: vcx, vcy ,vcz - real(8) :: vel - real(8) :: u, v, w - real(8) :: E - integer :: IE + real(8) :: vx ! velocity of neutron in x-direction + real(8) :: vy ! velocity of neutron in y-direction + real(8) :: vz ! velocity of neutron in z-direction + real(8) :: vcx ! velocity of CM in x-direction + real(8) :: vcy ! velocity of CM in y-direction + real(8) :: vcz ! velocity of CM in z-direction + real(8) :: vel ! magnitude of velocity + real(8) :: u ! x-direction + real(8) :: v ! y-direction + real(8) :: w ! z-direction + real(8) :: E ! energy vel = sqrt(p % E) awr = table % awr - vx = vel*p % uvw(1) - vy = vel*p % uvw(2) - vz = vel*p % uvw(3) + vx = vel * p%uvw(1) + vy = vel * p%uvw(2) + vz = vel * p%uvw(3) vcx = vx/(awr + ONE) vcy = vy/(awr + ONE) @@ -281,7 +299,7 @@ contains ! Sample scattering angle mu = sample_angle(rxn, p % E) - ! Determine direction cosines + ! Determine direction cosines in CM u = vx/vel v = vy/vel w = vz/vel @@ -322,23 +340,27 @@ contains type(AceReaction), pointer :: rxn integer :: i ! loop index - integer :: j ! index on nu energy grid + integer :: j ! index on nu energy grid / precursor group + integer :: k ! index on precursor yield grid integer :: loc ! index before start of energies/nu values integer :: NC ! number of polynomial coefficients integer :: NR ! number of interpolation regions integer :: NE ! number of energies tabulated + integer :: nu ! actual number of neutrons produced + integer :: law ! energy distribution law real(8) :: E ! incoming energy of neutron real(8) :: E_out ! outgoing energy of fission neutron - real(8) :: mu_out ! outgoing angle of fission neutron (if needed) real(8) :: c ! polynomial coefficient real(8) :: f ! interpolation factor real(8) :: nu_total ! total nu real(8) :: nu_prompt ! prompt nu real(8) :: nu_delay ! delayed nu - integer :: nu ! actual number of neutrons produced real(8) :: mu ! fission neutron angular cosine real(8) :: phi ! fission neutron azimuthal angle real(8) :: beta ! delayed neutron fraction + real(8) :: xi ! random number + real(8) :: yield ! delayed neutron precursor yield + real(8) :: prob ! cumulative probability character(250) :: msg ! error message ! copy energy of neutron @@ -502,14 +524,81 @@ contains ! sample cosine of angle mu = sample_angle(rxn, E) + ! sample between delayed and prompt neutrons + if (rang() < beta) then + ! ========================================================== + ! DELAYED NEUTRON SAMPLED + + ! sampled delayed precursor group + xi = rang() + loc = 1 + prob = ZERO + do j = 1, table % n_precursor + ! determine number of interpolation regions and energies + NR = table % nu_d_precursor_data(loc + 1) + NE = table % nu_d_precursor_data(loc + 2 + 2*NR) + if (NR > 0) then + msg = "Multiple interpolation regions not supported while & + &sampling delayed neutron precursor yield." + call error(msg) + end if + + ! interpolate on energy grid + loc = loc + 2 + 2*NR + if (E < table%nu_d_precursor_data(loc+1)) then + k = 1 + f = ZERO + elseif (E > table%nu_d_precursor_data(loc+NE)) then + k = NE - 1 + f = ONE + else + k = binary_search(table%nu_d_precursor_data(loc+1), NE, E) + f = (E - table%nu_d_precursor_data(loc+k)) / & + & (table%nu_d_precursor_data(loc+k+1) - & + & table%nu_d_precursor_data(loc+k)) + end if + + ! determine delayed neutron precursor yield for group j + loc = loc + NE + yield = table%nu_d_precursor_data(loc+k) + f * & + (table%nu_d_precursor_data(loc+k+1) - & + & table%nu_d_precursor_data(loc+k)) + prob = prob + yield + if (xi < prob) exit + + ! advance pointer + loc = loc + NE + 1 + end do + + ! sample from energy distribution for group j + law = table % nu_d_edist(j) % law + if (law == 44 .or. law == 61) then + call sample_energy(table%nu_d_edist(j), E, E_out, mu) + else + call sample_energy(table%nu_d_edist(j), E, E_out) + end if + + else + ! ========================================================== + ! PROMPT NEUTRON SAMPLED + + ! sample from prompt neutron energy distribution + law = rxn % edist % law + if (law == 44 .or. law == 61) then + call sample_energy(rxn%edist, E, E_out, prob) + else + call sample_energy(rxn%edist, E, E_out) + end if + + end if + ! Sample azimuthal angle uniformly in [0,2*pi) phi = TWO*PI*rang() fission_bank(i) % uvw(1) = mu - fission_bank(i) % uvw(2) = sqrt(ONE - mu**2) * cos(phi) - fission_bank(i) % uvw(3) = sqrt(ONE - mu**2) * sin(phi) + fission_bank(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi) + fission_bank(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi) - ! determine energy of fission neutron - call sample_energy(rxn, E, E_out, mu_out) + ! set energy of fission neutron fission_bank(i) % E = E_out end do @@ -522,7 +611,9 @@ contains end subroutine n_fission !===================================================================== -! INELASTIC_SCATTER +! INELASTIC_SCATTER handles all reactions with a single secondary +! neutron (other than fission), i.e. level scattering, (n,np), (n,na), +! etc. !===================================================================== subroutine inelastic_scatter(p, table, rxn) @@ -554,9 +645,9 @@ contains ! sample outgoing energy if (law == 44 .or. law == 61) then - call sample_energy(rxn, E_in, E, mu) + call sample_energy(rxn%edist, E_in, E, mu) else - call sample_energy(rxn, E_in, E) + call sample_energy(rxn%edist, E_in, E) end if ! if scattering system is in center-of-mass, transfer cosine of @@ -598,7 +689,7 @@ contains integer :: i ! loop index integer :: n_secondary ! number of secondary particles - integer :: law + integer :: law ! secondary energy distribution law real(8) :: A ! atomic weight ratio of nuclide real(8) :: E_in ! incoming energy real(8) :: mu ! cosine of scattering angle @@ -621,9 +712,9 @@ contains ! sample outgoing energy if (law == 44 .or. law == 61) then - call sample_energy(rxn, E_in, E, mu) + call sample_energy(rxn%edist, E_in, E, mu) else - call sample_energy(rxn, E_in, E) + call sample_energy(rxn%edist, E_in, E) end if ! if scattering system is in center-of-mass, transfer cosine of @@ -654,21 +745,22 @@ contains ! produce secondary particles in source bank n_secondary = abs(rxn % TY) - 1 do i = 1, n_secondary - ! create particle in source bank + ! TODO: create particle in source bank end do end subroutine n_xn !===================================================================== -! N_ABSORPTION +! N_ABSORPTION handles all absorbing reactions, i.e. (n,gamma), (n,p), +! (n,a), etc. !===================================================================== subroutine n_absorption(p) type(Particle), pointer :: p - integer :: cell_num - character(250) :: msg + integer :: cell_num ! user-specified cell number + character(250) :: msg ! output/error message p % alive = .false. if (verbosity >= 10) then @@ -848,12 +940,12 @@ contains ! SAMPLE_ENERGY !===================================================================== - subroutine sample_energy(rxn, E_in, E_out, mu_out) + subroutine sample_energy(edist, E_in, E_out, mu_out) - type(AceReaction), pointer :: rxn - real(8), intent(in) :: E_in - real(8), intent(out) :: E_out - real(8), intent(inout), optional :: mu_out + type(AceDistEnergy), intent(inout):: edist + real(8), intent(in) :: E_in + real(8), intent(out) :: E_out + real(8), intent(inout), optional :: mu_out integer :: i ! index on incoming energy grid integer :: k ! sampled index on outgoing grid @@ -900,23 +992,23 @@ contains ! law ! Check for multiple interpolation regions - if (rxn % edist % n_interp > 0) then + if (edist % n_interp > 0) then msg = "Multiple interpolation regions not supported while & &attempting to sampling secondary energy distribution." call error(msg) end if ! Determine which secondary energy distribution law to use - select case (rxn % edist % law) + select case (edist % law) case (1) ! ============================================================= ! TABULAR EQUIPROBABLE ENERGY BINS ! read number of interpolation regions, incoming energies, and ! outgoing energies - NR = rxn % edist % data(1) - NE = rxn % edist % data(2 + 2*NR) - NET = rxn % edist % data(3 + 2*NR + NE) + NR = edist % data(1) + NE = edist % data(2 + 2*NR) + NET = edist % data(3 + 2*NR + NE) if (NR > 0) then msg = "Multiple interpolation regions not supported while & &attempting to sample equiprobable energy bins." @@ -926,9 +1018,9 @@ contains ! determine index on incoming energy grid and interpolation ! factor loc = 2 + 2*NR - i = binary_search(rxn % edist % data(loc+1), NE, E_in) - r = (E_in - rxn%edist%data(loc+i)) / & - & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + i = binary_search(edist % data(loc+1), NE, E_in) + r = (E_in - edist%data(loc+i)) / & + & (edist%data(loc+i+1) - edist%data(loc+i)) ! Sample outgoing energy bin xi1 = rang() @@ -943,8 +1035,8 @@ contains end if loc = 3 + 2*NR + NE + (l-1)*NET - E_l_k = rxn % edist % data(loc+k) - E_l_k1 = rxn % edist % data(loc+k+1) + E_l_k = edist % data(loc+k) + E_l_k1 = edist % data(loc+k+1) xi2 = rang() E_out = E_l_k + xi2*(E_l_k1 - E_l_k) @@ -954,7 +1046,7 @@ contains ! ============================================================= ! INELASTIC LEVEL SCATTERING - E_cm = rxn%edist%data(2) * (E_in - rxn%edist%data(1)) + E_cm = edist%data(2) * (E_in - edist%data(1)) E_out = E_cm @@ -963,8 +1055,8 @@ contains ! CONTINUOUS TABULAR DISTRIBUTION ! read number of interpolation regions and incoming energies - NR = rxn % edist % data(1) - NE = rxn % edist % data(2 + 2*NR) + NR = edist % data(1) + NE = edist % data(2 + 2*NR) if (NR > 0) then msg = "Multiple interpolation regions not supported while & &attempting to sample continuous tabular distribution." @@ -975,16 +1067,16 @@ contains ! energy is outside the range of the tabulated energies, choose ! the first or last bins loc = 2 + 2*NR - if (E_in < rxn % edist % data(loc+1)) then + if (E_in < edist % data(loc+1)) then i = 1 r = ZERO - elseif (E_in > rxn % edist % data(loc+NE)) then + elseif (E_in > edist % data(loc+NE)) then i = NE - 1 r = ONE else - i = binary_search(rxn % edist % data(loc+1), NE, E_in) - r = (E_in - rxn%edist%data(loc+i)) / & - & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + i = binary_search(edist % data(loc+1), NE, E_in) + r = (E_in - edist%data(loc+i)) / & + & (edist%data(loc+i+1) - edist%data(loc+i)) end if ! Sample between the ith and (i+1)th bin @@ -996,25 +1088,25 @@ contains end if ! interpolation for energy E1 and EK - loc = rxn%edist%data(2 + 2*NR + NE + i) - NP = rxn%edist%data(loc + 2) - E_i_1 = rxn%edist%data(loc + 2 + 1) - E_i_K = rxn%edist%data(loc + 2 + NP) + loc = edist%data(2 + 2*NR + NE + i) + NP = edist%data(loc + 2) + E_i_1 = edist%data(loc + 2 + 1) + E_i_K = edist%data(loc + 2 + NP) - loc = rxn%edist%data(2 + 2*NR + NE + i + 1) - NP = rxn%edist%data(loc + 2) - E_i1_1 = rxn%edist%data(loc + 2 + 1) - E_i1_K = rxn%edist%data(loc + 2 + NP) + loc = edist%data(2 + 2*NR + NE + i + 1) + NP = edist%data(loc + 2) + E_i1_1 = edist%data(loc + 2 + 1) + E_i1_K = edist%data(loc + 2 + NP) E_1 = E_i_1 + r*(E_i1_1 - E_i_1) E_K = E_i_K + r*(E_i1_K - E_i_K) ! determine location of outgoing energies, pdf, cdf for E(l) - loc = rxn % edist % data(2 + 2*NR + NE + l) + loc = edist % data(2 + 2*NR + NE + l) ! determine type of interpolation and number of discrete lines - INTTp = rxn % edist % data(loc + 1) - NP = rxn % edist % data(loc + 2) + INTTp = edist % data(loc + 1) + NP = edist % data(loc + 2) if (INTTp > 10) then INTT = mod(INTTp,10) ND = (INTTp - INTT)/10 @@ -1033,15 +1125,15 @@ contains ! determine outgoing energy bin xi1 = rang() loc = loc + 2 ! start of EOUT - c_k = rxn % edist % data(loc + 2*NP + 1) + c_k = edist % data(loc + 2*NP + 1) do k = 1, NP-1 - c_k1 = rxn % edist % data(loc + 2*NP + k+1) + c_k1 = edist % data(loc + 2*NP + k+1) if (xi1 < c_k1) exit c_k = c_k1 end do - E_l_k = rxn % edist % data(loc+k) - p_l_k = rxn % edist % data(loc+NP+k) + E_l_k = edist % data(loc+k) + p_l_k = edist % data(loc+NP+k) if (INTT == HISTOGRAM) then ! Histogram interpolation E_out = E_l_k + (xi1 - c_k)/p_l_k @@ -1049,8 +1141,8 @@ contains elseif (INTT == LINEAR_LINEAR) then ! Linear-linear interpolation -- not sure how you come about ! the formula given in the MCNP manual - E_l_k1 = rxn % edist % data(loc+k+1) - p_l_k1 = rxn % edist % data(loc+NP+k+1) + E_l_k1 = edist % data(loc+k+1) + p_l_k1 = edist % data(loc+NP+k+1) frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) if (frac == ZERO) then @@ -1080,8 +1172,8 @@ contains ! MAXWELL FISSION SPECTRUM ! read number of interpolation regions and incoming energies - NR = rxn % edist % data(1) - NE = rxn % edist % data(2 + 2*NR) + NR = edist % data(1) + NE = edist % data(2 + 2*NR) if (NR > 0) then msg = "Multiple interpolation regions not supported while & &attempting to sample Maxwell fission spectrum." @@ -1090,22 +1182,22 @@ contains ! find incident energy bin and calculate interpolation factor loc = 2 + 2*NR - if (E_in < rxn % edist % data(loc+1)) then + if (E_in < edist % data(loc+1)) then i = 1 r = ZERO - elseif (E_in > rxn % edist % data(loc+NE)) then + elseif (E_in > edist % data(loc+NE)) then i = NE - 1 r = ONE else - i = binary_search(rxn % edist % data(loc+1), NE, E_in) - r = (E_in - rxn%edist%data(loc+i)) / & - & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + i = binary_search(edist % data(loc+1), NE, E_in) + r = (E_in - edist%data(loc+i)) / & + & (edist%data(loc+i+1) - edist%data(loc+i)) end if ! determine nuclear temperature from tabulated function loc = loc + NE - T = rxn%edist%data(loc+i) + r * & - & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + T = edist%data(loc+i) + r * & + & (edist%data(loc+i+1) - edist%data(loc+i)) ! sample maxwell fission spectrum E_out = maxwell_spectrum(T) @@ -1115,8 +1207,8 @@ contains ! EVAPORATION SPECTRUM ! read number of interpolation regions and incoming energies - NR = rxn % edist % data(1) - NE = rxn % edist % data(2 + 2*NR) + NR = edist % data(1) + NE = edist % data(2 + 2*NR) if (NR > 0) then msg = "Multiple interpolation regions not supported while & &attempting to sample evaporation spectrum." @@ -1127,26 +1219,26 @@ contains ! energy is outside the range of the tabulated energies, choose ! the first or last bins loc = 2 + 2*NR - if (E_in < rxn % edist % data(loc+1)) then + if (E_in < edist % data(loc+1)) then i = 1 r = ZERO - elseif (E_in > rxn % edist % data(loc+NE)) then + elseif (E_in > edist % data(loc+NE)) then i = NE - 1 r = ONE else - i = binary_search(rxn % edist % data(loc+1), NE, E_in) - r = (E_in - rxn%edist%data(loc+i)) / & - & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + i = binary_search(edist % data(loc+1), NE, E_in) + r = (E_in - edist%data(loc+i)) / & + & (edist%data(loc+i+1) - edist%data(loc+i)) end if ! determine nuclear temperature from tabulated function loc = loc + NE - T = rxn%edist%data(loc+i) + r * & - & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + T = edist%data(loc+i) + r * & + & (edist%data(loc+i+1) - edist%data(loc+i)) ! determine restriction energy loc = loc + NE - U = rxn % edist % data(loc + 1) + U = edist % data(loc + 1) ! sample outgoing energy based on evaporation spectrum ! probability density function @@ -1163,8 +1255,8 @@ contains ! read number of interpolation regions and incoming energies ! for parameter 'a' - NR = rxn % edist % data(1) - NE = rxn % edist % data(2 + 2*NR) + NR = edist % data(1) + NE = edist % data(2 + 2*NR) if (NR > 0) then msg = "Multiple interpolation regions not supported while & &attempting to sample Watt fission spectrum." @@ -1173,28 +1265,28 @@ contains ! find incident energy bin and calculate interpolation factor loc = 2 + 2*NR - if (E_in < rxn % edist % data(loc+1)) then + if (E_in < edist % data(loc+1)) then i = 1 r = ZERO - elseif (E_in > rxn % edist % data(loc+NE)) then + elseif (E_in > edist % data(loc+NE)) then i = NE - 1 r = ONE else - i = binary_search(rxn % edist % data(loc+1), NE, E_in) - r = (E_in - rxn%edist%data(loc+i)) / & - & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + i = binary_search(edist % data(loc+1), NE, E_in) + r = (E_in - edist%data(loc+i)) / & + & (edist%data(loc+i+1) - edist%data(loc+i)) end if ! determine Watt parameter 'a' from tabulated function loc = loc + NE - Watt_a = rxn%edist%data(loc+i) + r * & - & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + Watt_a = edist%data(loc+i) + r * & + & (edist%data(loc+i+1) - edist%data(loc+i)) ! read number of interpolation regions and incoming energies ! for parameter 'b' loc = loc + NE - NR = rxn % edist % data(loc + 1) - NE = rxn % edist % data(loc + 2 + 2*NR) + NR = edist % data(loc + 1) + NE = edist % data(loc + 2 + 2*NR) if (NR > 0) then msg = "Multiple interpolation regions not supported while & &attempting to sample Watt fission spectrum." @@ -1203,22 +1295,22 @@ contains ! find incident energy bin and calculate interpolation factor loc = loc + 2 + 2*NR - if (E_in < rxn % edist % data(loc+1)) then + if (E_in < edist % data(loc+1)) then i = 1 r = ZERO - elseif (E_in > rxn % edist % data(loc+NE)) then + elseif (E_in > edist % data(loc+NE)) then i = NE - 1 r = ONE else - i = binary_search(rxn % edist % data(loc+1), NE, E_in) - r = (E_in - rxn%edist%data(loc+i)) / & - & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + i = binary_search(edist % data(loc+1), NE, E_in) + r = (E_in - edist%data(loc+i)) / & + & (edist%data(loc+i+1) - edist%data(loc+i)) end if ! determine Watt parameter 'b' from tabulated function loc = loc + NE - Watt_b = rxn%edist%data(loc+i) + r * & - & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + Watt_b = edist%data(loc+i) + r * & + & (edist%data(loc+i+1) - edist%data(loc+i)) ! Sample energy-dependent Watt fission spectrum E_out = watt_spectrum(Watt_a, Watt_b) @@ -1233,8 +1325,8 @@ contains end if ! read number of interpolation regions and incoming energies - NR = rxn % edist % data(1) - NE = rxn % edist % data(2 + 2*NR) + NR = edist % data(1) + NE = edist % data(2 + 2*NR) if (NR > 0) then msg = "Multiple interpolation regions not supported while & &attempting to sample Kalbach-Mann distribution." @@ -1245,16 +1337,16 @@ contains ! energy is outside the range of the tabulated energies, choose ! the first or last bins loc = 2 + 2*NR - if (E_in < rxn % edist % data(loc+1)) then + if (E_in < edist % data(loc+1)) then i = 1 r = ZERO - elseif (E_in > rxn % edist % data(loc+NE)) then + elseif (E_in > edist % data(loc+NE)) then i = NE - 1 r = ONE else - i = binary_search(rxn % edist % data(loc+1), NE, E_in) - r = (E_in - rxn%edist%data(loc+i)) / & - & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + i = binary_search(edist % data(loc+1), NE, E_in) + r = (E_in - edist%data(loc+i)) / & + & (edist%data(loc+i+1) - edist%data(loc+i)) end if ! Sample between the ith and (i+1)th bin @@ -1266,26 +1358,26 @@ contains end if ! determine endpoints on grid i - loc = rxn%edist%data(2+2*NR+NE + i) ! start of LDAT for i - NP = rxn%edist%data(loc + 2) - E_i_1 = rxn%edist%data(loc + 2 + 1) - E_i_K = rxn%edist%data(loc + 2 + NP) + loc = edist%data(2+2*NR+NE + i) ! start of LDAT for i + NP = edist%data(loc + 2) + E_i_1 = edist%data(loc + 2 + 1) + E_i_K = edist%data(loc + 2 + NP) ! determine endpoints on grid i+1 - loc = rxn%edist%data(2+2*NR+NE + i+1) ! start of LDAT for i+1 - NP = rxn%edist%data(loc + 2) - E_i1_1 = rxn%edist%data(loc + 2 + 1) - E_i1_K = rxn%edist%data(loc + 2 + NP) + loc = edist%data(2+2*NR+NE + i+1) ! start of LDAT for i+1 + NP = edist%data(loc + 2) + E_i1_1 = edist%data(loc + 2 + 1) + E_i1_K = edist%data(loc + 2 + NP) E_1 = E_i_1 + r*(E_i1_1 - E_i_1) E_K = E_i_K + r*(E_i1_K - E_i_K) ! determine location of outgoing energies, pdf, cdf for E(l) - loc = rxn % edist % data(2 + 2*NR + NE + l) + loc = edist % data(2 + 2*NR + NE + l) ! determine type of interpolation and number of discrete lines - INTTp = rxn % edist % data(loc + 1) - NP = rxn % edist % data(loc + 2) + INTTp = edist % data(loc + 1) + NP = edist % data(loc + 2) if (INTTp > 10) then INTT = mod(INTTp,10) ND = (INTTp - INTT)/10 @@ -1304,28 +1396,28 @@ contains ! determine outgoing energy bin xi1 = rang() loc = loc + 2 ! start of EOUT - c_k = rxn % edist % data(loc + 2*NP + 1) + c_k = edist % data(loc + 2*NP + 1) do k = 1, NP-1 - c_k1 = rxn % edist % data(loc + 2*NP + k+1) + c_k1 = edist % data(loc + 2*NP + k+1) if (xi1 < c_k1) exit c_k = c_k1 end do - E_l_k = rxn % edist % data(loc+k) - p_l_k = rxn % edist % data(loc+NP+k) + E_l_k = edist % data(loc+k) + p_l_k = edist % data(loc+NP+k) if (INTT == HISTOGRAM) then ! Histogram interpolation E_out = E_l_k + (xi1 - c_k)/p_l_k ! Determine Kalbach-Mann parameters - KM_R = rxn % edist % data(loc + 3*NP + k) - KM_A = rxn % edist % data(loc + 4*NP + k) + KM_R = edist % data(loc + 3*NP + k) + KM_A = edist % data(loc + 4*NP + k) elseif (INTT == LINEAR_LINEAR) then ! Linear-linear interpolation -- not sure how you come about ! the formula given in the MCNP manual - E_l_k1 = rxn % edist % data(loc+k+1) - p_l_k1 = rxn % edist % data(loc+NP+k+1) + E_l_k1 = edist % data(loc+k+1) + p_l_k1 = edist % data(loc+NP+k+1) ! Find E prime frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) @@ -1337,10 +1429,10 @@ contains end if ! Determine Kalbach-Mann parameters - R_k = rxn % edist % data(loc + 3*NP + k) - R_k1 = rxn % edist % data(loc + 3*NP + k+1) - A_k = rxn % edist % data(loc + 4*NP + k) - A_k1 = rxn % edist % data(loc + 4*NP + k+1) + R_k = edist % data(loc + 3*NP + k) + R_k1 = edist % data(loc + 3*NP + k+1) + A_k = edist % data(loc + 4*NP + k) + A_k1 = edist % data(loc + 4*NP + k+1) KM_R = R_k + (R_k1 - R_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k) KM_A = A_k + (A_k1 - A_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k) @@ -1376,8 +1468,8 @@ contains end if ! read number of interpolation regions and incoming energies - NR = rxn % edist % data(1) - NE = rxn % edist % data(2 + 2*NR) + NR = edist % data(1) + NE = edist % data(2 + 2*NR) if (NR > 0) then msg = "Multiple interpolation regions not supported while & &attempting to sample correlated energy-angle distribution." @@ -1388,16 +1480,16 @@ contains ! energy is outside the range of the tabulated energies, choose ! the first or last bins loc = 2 + 2*NR - if (E_in < rxn % edist % data(loc+1)) then + if (E_in < edist % data(loc+1)) then i = 1 r = ZERO - elseif (E_in > rxn % edist % data(loc+NE)) then + elseif (E_in > edist % data(loc+NE)) then i = NE - 1 r = ONE else - i = binary_search(rxn % edist % data(loc+1), NE, E_in) - r = (E_in - rxn%edist%data(loc+i)) / & - & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + i = binary_search(edist % data(loc+1), NE, E_in) + r = (E_in - edist%data(loc+i)) / & + & (edist%data(loc+i+1) - edist%data(loc+i)) end if ! Sample between the ith and (i+1)th bin @@ -1409,26 +1501,26 @@ contains end if ! determine endpoints on grid i - loc = rxn%edist%data(2+2*NR+NE + i) ! start of LDAT for i - NP = rxn%edist%data(loc + 2) - E_i_1 = rxn%edist%data(loc + 2 + 1) - E_i_K = rxn%edist%data(loc + 2 + NP) + loc = edist%data(2+2*NR+NE + i) ! start of LDAT for i + NP = edist%data(loc + 2) + E_i_1 = edist%data(loc + 2 + 1) + E_i_K = edist%data(loc + 2 + NP) ! determine endpoints on grid i+1 - loc = rxn%edist%data(2+2*NR+NE + i+1) ! start of LDAT for i+1 - NP = rxn%edist%data(loc + 2) - E_i1_1 = rxn%edist%data(loc + 2 + 1) - E_i1_K = rxn%edist%data(loc + 2 + NP) + loc = edist%data(2+2*NR+NE + i+1) ! start of LDAT for i+1 + NP = edist%data(loc + 2) + E_i1_1 = edist%data(loc + 2 + 1) + E_i1_K = edist%data(loc + 2 + NP) E_1 = E_i_1 + r*(E_i1_1 - E_i_1) E_K = E_i_K + r*(E_i1_K - E_i_K) ! determine location of outgoing energies, pdf, cdf for E(l) - loc = rxn % edist % data(2 + 2*NR + NE + l) + loc = edist % data(2 + 2*NR + NE + l) ! determine type of interpolation and number of discrete lines - INTTp = rxn % edist % data(loc + 1) - NP = rxn % edist % data(loc + 2) + INTTp = edist % data(loc + 1) + NP = edist % data(loc + 2) if (INTTp > 10) then INTT = mod(INTTp,10) ND = (INTTp - INTT)/10 @@ -1447,15 +1539,15 @@ contains ! determine outgoing energy bin xi1 = rang() loc = loc + 2 ! start of EOUT - c_k = rxn % edist % data(loc + 2*NP + 1) + c_k = edist % data(loc + 2*NP + 1) do k = 1, NP-1 - c_k1 = rxn % edist % data(loc + 2*NP + k+1) + c_k1 = edist % data(loc + 2*NP + k+1) if (xi1 < c_k1) exit c_k = c_k1 end do - E_l_k = rxn % edist % data(loc+k) - p_l_k = rxn % edist % data(loc+NP+k) + E_l_k = edist % data(loc+k) + p_l_k = edist % data(loc+NP+k) if (INTT == HISTOGRAM) then ! Histogram interpolation E_out = E_l_k + (xi1 - c_k)/p_l_k @@ -1463,8 +1555,8 @@ contains elseif (INTT == LINEAR_LINEAR) then ! Linear-linear interpolation -- not sure how you come about ! the formula given in the MCNP manual - E_l_k1 = rxn % edist % data(loc+k+1) - p_l_k1 = rxn % edist % data(loc+NP+k+1) + E_l_k1 = edist % data(loc+k+1) + p_l_k1 = edist % data(loc+NP+k+1) ! Find E prime frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) @@ -1487,7 +1579,7 @@ contains end if ! Find location of correlated angular distribution - loc = rxn % edist % data(loc+3*NP+k) + loc = edist % data(loc+3*NP+k) ! Check if angular distribution is isotropic if (loc == 0) then @@ -1496,21 +1588,21 @@ contains end if ! interpolation type and number of points in angular distribution - JJ = rxn % edist % data(loc + 1) - NP = rxn % edist % data(loc + 2) + JJ = edist % data(loc + 1) + NP = edist % data(loc + 2) ! determine outgoing cosine bin xi3 = rang() loc = loc + 2 - c_k = rxn % edist % data(loc + 2*NP + 1) + c_k = edist % data(loc + 2*NP + 1) do k = 1, NP-1 - c_k1 = rxn % edist % data(loc + 2*NP + k+1) + c_k1 = edist % data(loc + 2*NP + k+1) if (xi3 < c_k1) exit c_k = c_k1 end do - p_k = rxn % edist % data(loc + NP + k) - mu_k = rxn % edist % data(loc + k) + p_k = edist % data(loc + NP + k) + mu_k = edist % data(loc + k) if (JJ == HISTOGRAM) then ! Histogram interpolation mu_out = mu_k + (xi3 - c_k)/p_k @@ -1518,8 +1610,8 @@ contains elseif (JJ == LINEAR_LINEAR) then ! Linear-linear interpolation -- not sure how you come about ! the formula given in the MCNP manual - p_k1 = rxn % edist % data(loc + NP + k+1) - mu_k1 = rxn % edist % data(loc + k+1) + p_k1 = edist % data(loc + NP + k+1) + mu_k1 = edist % data(loc + k+1) frac = (p_k1 - p_k)/(mu_k1 - mu_k) if (frac == ZERO) then @@ -1558,6 +1650,8 @@ contains real(8) :: r1, r2, r3, r4 ! random numbers real(8) :: d ! r1^2 + r2^2 + ! TODO: Change this to direct sampling + r1 = rang() do r2 = rang() @@ -1568,7 +1662,7 @@ contains r3 = rang() r4 = rang() - E_out = -T*(r1**2 * log(r3) / d + log(r4)) + E_out = -T*(r1*r1 * log(r3) / d + log(r4)) end function maxwell_spectrum @@ -1581,12 +1675,12 @@ contains function watt_spectrum(a, b) result(E_out) - real(8), intent(in) :: a - real(8), intent(in) :: b - real(8) :: E_out + real(8), intent(in) :: a ! Watt parameter a + real(8), intent(in) :: b ! Watt parameter b + real(8) :: E_out ! energy of emitted neutron - real(8) :: g - real(8) :: r1, r2 + real(8) :: g ! from sampling scheme R28 + real(8) :: r1, r2 ! random numbers g = sqrt((1 + a*b/8)**2 - 1) + (1 + a*b/8) do diff --git a/src/types.f90 b/src/types.f90 index bedeb27cf..decc03908 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -159,7 +159,7 @@ module types integer :: n_interp ! # of interpolation regions integer, allocatable :: nbt(:) ! ENDF interpolation parameters integer, allocatable :: int(:) ! '' - integer :: n_energy ! # of incoming energies + integer :: n_energy ! # of energies for law validity real(8), allocatable :: energy(:) ! energy grid for law validity real(8), allocatable :: pvalid(:) ! probability of law validity real(8), allocatable :: data(:) ! energy distribution data @@ -211,11 +211,10 @@ module types ! Delayed fission neutron emission integer :: nu_d_type + integer :: n_precursor real(8), allocatable :: nu_d_data(:) - type(AceDistEnergy) :: nu_d_edist - real(8), allocatable :: nu_d_precursor_const(:,:) - real(8), allocatable :: nu_d_precursor_energy(:,:) - real(8), allocatable :: nu_d_precursor_prob(:,:) + real(8), allocatable :: nu_d_precursor_data(:) + type(AceDistEnergy), allocatable :: nu_d_edist(:) ! Reactions integer :: n_reaction