diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 51e8a1aada..cd1d95fd17 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -369,12 +369,16 @@ contains ! allocate scoring bins allocate(t % score_bins(4)) t % n_score_bins = 4 - + ! allocate scattering order data + allocate(t % scatt_order(4)) + t % scatt_order = 0 + ! set macro_bins - t % score_bins(1) = SCORE_FLUX - t % score_bins(2) = SCORE_TOTAL - t % score_bins(3) = SCORE_SCATTER_1 - t % score_bins(4) = SCORE_DIFFUSION + t % score_bins(1) = SCORE_FLUX + t % score_bins(2) = SCORE_TOTAL + t % score_bins(3) = SCORE_SCATTER_N + t % scatt_order(3) = 1 + t % score_bins(4) = SCORE_DIFFUSION ! Increment the appropriate index and set pointer analog_tallies(n_user_analog_tallies + 1) = i @@ -413,6 +417,9 @@ contains ! allocate macro reactions allocate(t % score_bins(2)) t % n_score_bins = 2 + ! allocate scattering order data + allocate(t % scatt_order(2)) + t % scatt_order = 0 ! set macro_bins t % score_bins(1) = SCORE_NU_SCATTER @@ -454,6 +461,9 @@ contains ! allocate macro reactions allocate(t % score_bins(1)) t % n_score_bins = 1 + ! allocate scattering order data + allocate(t % scatt_order(1)) + t % scatt_order = 0 ! set macro bins t % score_bins(1) = SCORE_CURRENT diff --git a/src/constants.F90 b/src/constants.F90 index e1d8e853a8..eebb90df67 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -251,28 +251,25 @@ module constants EVENT_FISSION = 3 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 20 + integer, parameter :: N_SCORE_TYPES = 17 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate SCORE_SCATTER = -3, & ! scattering rate SCORE_NU_SCATTER = -4, & ! scattering production rate - SCORE_SCATTER_1 = -5, & ! first scattering moment - SCORE_SCATTER_2 = -6, & ! second scattering moment - SCORE_SCATTER_3 = -7, & ! third scattering moment - SCORE_SCATTER_4 = -8, & ! fourth scattering moment - SCORE_SCATTER_5 = -9, & ! fifth scattering moment - SCORE_TRANSPORT = -10, & ! transport reaction rate - SCORE_DIFFUSION = -11, & ! diffusion coefficient - SCORE_N_1N = -12, & ! (n,1n) rate - SCORE_N_2N = -13, & ! (n,2n) rate - SCORE_N_3N = -14, & ! (n,3n) rate - SCORE_N_4N = -15, & ! (n,4n) rate - SCORE_ABSORPTION = -16, & ! absorption rate - SCORE_FISSION = -17, & ! fission rate - SCORE_NU_FISSION = -18, & ! neutron production rate - SCORE_CURRENT = -19, & ! partial current - SCORE_EVENTS = -20 ! number of events + SCORE_SCATTER_N = -5, & ! arbitrary scattering moment + SCORE_SCATTER_PN = -6, & ! system for scoring 0th through nth moment + SCORE_TRANSPORT = -7, & ! transport reaction rate + SCORE_DIFFUSION = -8, & ! diffusion coefficient + SCORE_N_1N = -9, & ! (n,1n) rate + SCORE_N_2N = -10, & ! (n,2n) rate + SCORE_N_3N = -11, & ! (n,3n) rate + SCORE_N_4N = -12, & ! (n,4n) rate + SCORE_ABSORPTION = -13, & ! absorption rate + SCORE_FISSION = -14, & ! fission rate + SCORE_NU_FISSION = -15, & ! neutron production rate + SCORE_CURRENT = -16, & ! partial current + SCORE_EVENTS = -17 ! number of events ! Tally map bin finding integer, parameter :: NO_BIN_FOUND = -1 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index fd5b90ebda..38873162bc 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -19,7 +19,7 @@ module input_xml implicit none type(DictionaryII), pointer :: & ! used to count how many cells each - cells_in_univ_dict => null() ! universe contains + cells_in_univ_dict => null() ! universe contains contains @@ -92,7 +92,7 @@ contains call get_environment_variable("CROSS_SECTIONS", env_variable) if (len_trim(env_variable) == 0) then message = "No cross_sections.xml file was specified in settings.xml & - &or in the CROSS_SECTIONS environment variable." + &or in the CROSS_SECTIONS environment variable." call fatal_error() else path_cross_sections = trim(env_variable) @@ -104,7 +104,7 @@ contains ! Make sure that either criticality or fixed source was specified if (criticality_ % batches == 0 .and. fixed_source_ % batches == 0) then message = "Number of batches on or " & - // "tag was zero." + // "tag was zero." call fatal_error() end if @@ -202,7 +202,7 @@ contains inquire(FILE=path_source, EXIST=file_exists) if (.not. file_exists) then message = "Binary source file '" // trim(path_source) // & - "' does not exist!" + "' does not exist!" call fatal_error() end if @@ -221,7 +221,7 @@ contains coeffs_reqd = 3 case default message = "Invalid spatial distribution for external source: " & - // trim(source_ % space % type) + // trim(source_ % space % type) call fatal_error() end select @@ -235,11 +235,11 @@ contains ! Read parameters for spatial distribution if (n < coeffs_reqd) then message = "Not enough parameters specified for spatial " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > coeffs_reqd) then message = "Too many parameters specified for spatial " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > 0) then allocate(external_source % params_space(n)) @@ -266,7 +266,7 @@ contains external_source % type_angle = SRC_ANGLE_TABULAR case default message = "Invalid angular distribution for external source: " & - // trim(source_ % angle % type) + // trim(source_ % angle % type) call fatal_error() end select @@ -280,11 +280,11 @@ contains ! Read parameters for angle distribution if (n < coeffs_reqd) then message = "Not enough parameters specified for angle " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > coeffs_reqd) then message = "Too many parameters specified for angle " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > 0) then allocate(external_source % params_angle(n)) @@ -314,7 +314,7 @@ contains external_source % type_energy = SRC_ENERGY_TABULAR case default message = "Invalid energy distribution for external source: " & - // trim(source_ % energy % type) + // trim(source_ % energy % type) call fatal_error() end select @@ -328,11 +328,11 @@ contains ! Read parameters for energy distribution if (n < coeffs_reqd) then message = "Not enough parameters specified for energy " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > coeffs_reqd) then message = "Too many parameters specified for energy " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > 0) then allocate(external_source % params_energy(n)) @@ -370,10 +370,10 @@ contains ! Check to make sure enough values were supplied if (size(entropy_(1) % lower_left) /= 3) then message = "Need to specify (x,y,z) coordinates of lower-left corner & - &of Shannon entropy mesh." + &of Shannon entropy mesh." elseif (size(entropy_(1) % upper_right) /= 3) then message = "Need to specify (x,y,z) coordinates of upper-right corner & - &of Shannon entropy mesh." + &of Shannon entropy mesh." end if ! Allocate mesh object and coordinates on mesh @@ -388,7 +388,7 @@ contains ! Check on values provided if (.not. all(entropy_mesh % upper_right > entropy_mesh % lower_left)) then message = "Upper-right coordinate must be greater than lower-left & - &coordinate for Shannon entropy mesh." + &coordinate for Shannon entropy mesh." call fatal_error() end if @@ -399,7 +399,7 @@ contains ! If so, make sure proper number of values were given if (size(entropy_(1) % dimension) /= 3) then message = "Dimension of entropy mesh must be given as three & - &integers." + &integers." call fatal_error() end if @@ -420,13 +420,13 @@ contains ! Check to make sure enough values were supplied if (size(uniform_fs_(1) % lower_left) /= 3) then message = "Need to specify (x,y,z) coordinates of lower-left corner & - &of UFS mesh." + &of UFS mesh." elseif (size(uniform_fs_(1) % upper_right) /= 3) then message = "Need to specify (x,y,z) coordinates of upper-right corner & - &of UFS mesh." + &of UFS mesh." elseif (size(uniform_fs_(1) % dimension) /= 3) then message = "Dimension of UFS mesh must be given as three & - &integers." + &integers." call fatal_error() end if @@ -450,13 +450,13 @@ contains ! Check on values provided if (.not. all(ufs_mesh % upper_right > ufs_mesh % lower_left)) then message = "Upper-right coordinate must be greater than lower-left & - &coordinate for UFS mesh." + &coordinate for UFS mesh." call fatal_error() end if ! Calculate width ufs_mesh % width = (ufs_mesh % upper_right - & - ufs_mesh % lower_left) / ufs_mesh % dimension + ufs_mesh % lower_left) / ufs_mesh % dimension ! Calculate volume fraction of each cell ufs_mesh % volume_frac = ONE/real(product(ufs_mesh % dimension),8) @@ -466,7 +466,7 @@ contains ! Allocate source_frac allocate(source_frac(1, ufs_mesh % dimension(1), & - ufs_mesh % dimension(2), ufs_mesh % dimension(3))) + ufs_mesh % dimension(2), ufs_mesh % dimension(3))) end if ! Check if the user has specified to write state points @@ -488,7 +488,7 @@ contains n_state_points = n_batches / state_point_(1) % interval allocate(statepoint_batch(n_state_points)) statepoint_batch = (/ (state_point_(1) % interval * i, i = 1, & - n_state_points) /) + n_state_points) /) else ! If neither were specified, write state point at last batch n_state_points = 1 @@ -498,7 +498,7 @@ contains ! Check if the user has specified to write binary source file if (trim(state_point_(1) % source_separate) == 'on') & - source_separate = .true. + source_separate = .true. else ! If no tag was present, by default write state point at ! last batch only @@ -629,7 +629,7 @@ contains ! Check to make sure that either material or fill was specified if (c % material == NONE .and. c % fill == NONE) then message = "Neither material nor fill was specified for cell " // & - trim(to_str(c % id)) + trim(to_str(c % id)) call fatal_error() end if @@ -643,7 +643,7 @@ contains ! Check to make sure that surfaces were specified if (.not. associated(cell_(i) % surfaces)) then message = "No surfaces specified for cell " // & - trim(to_str(c % id)) + trim(to_str(c % id)) call fatal_error() end if @@ -659,7 +659,7 @@ contains ! another universe if (c % fill == NONE) then message = "Cannot apply a rotation to cell " // trim(to_str(& - c % id)) // " because it is not filled with another universe" + c % id)) // " because it is not filled with another universe" call fatal_error() end if @@ -667,7 +667,7 @@ contains n = size(cell_(i) % rotation) if (n /= 3) then message = "Incorrect number of rotation parameters on cell " // & - to_str(c % id) + to_str(c % id) call fatal_error() end if @@ -679,13 +679,13 @@ contains ! Calculate rotation matrix based on angles given allocate(c % rotation(3,3)) c % rotation = reshape((/ & - cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), & - -cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), & - cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), & - sin(phi)*cos(theta), & - sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), & - -sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), & - cos(phi)*cos(theta) /), (/ 3,3 /)) + cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), & + -cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), & + cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), & + sin(phi)*cos(theta), & + sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), & + -sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), & + cos(phi)*cos(theta) /), (/ 3,3 /)) end if ! Translation vector @@ -694,7 +694,7 @@ contains ! another universe if (c % fill == NONE) then message = "Cannot apply a translation to cell " // trim(to_str(& - c % id)) // " because it is not filled with another universe" + c % id)) // " because it is not filled with another universe" call fatal_error() end if @@ -702,7 +702,7 @@ contains n = size(cell_(i) % translation) if (n /= 3) then message = "Incorrect number of translation parameters on cell " & - // to_str(c % id) + // to_str(c % id) call fatal_error() end if @@ -814,11 +814,11 @@ contains n = size(surface_(i) % coeffs) if (n < coeffs_reqd) then message = "Not enough coefficients specified for surface: " // & - trim(to_str(s % id)) + trim(to_str(s % id)) call fatal_error() elseif (n > coeffs_reqd) then message = "Too many coefficients specified for surface: " // & - trim(to_str(s % id)) + trim(to_str(s % id)) call fatal_error() else allocate(s % coeffs(n)) @@ -839,7 +839,7 @@ contains s % bc = BC_PERIODIC case default message = "Unknown boundary condition '" // trim(word) // & - "' specified on surface " // trim(to_str(s % id)) + "' specified on surface " // trim(to_str(s % id)) call fatal_error() end select @@ -888,7 +888,7 @@ contains ! Read lattice lower-left location if (size(lattice_(i) % dimension) /= size(lattice_(i) % lower_left)) then message = "Number of entries on must be the same as " // & - "the number of entries on ." + "the number of entries on ." call fatal_error() end if l % x0 = lattice_(i) % lower_left(1) @@ -897,7 +897,7 @@ contains ! Read lattice widths if (size(lattice_(i) % width) /= size(lattice_(i) % lower_left)) then message = "Number of entries on must be the same as " // & - "the number of entries on ." + "the number of entries on ." call fatal_error() end if l % width_x = lattice_(i) % width(1) @@ -998,7 +998,7 @@ contains sum_density = .false. if (val <= ZERO) then message = "Need to specify a positive density on material " // & - trim(to_str(mat % id)) // "." + trim(to_str(mat % id)) // "." call fatal_error() end if @@ -1015,7 +1015,7 @@ contains mat % density = 1.0e-24 * val case default message = "Unkwown units '" // trim(material_(i) % density % units) & - // "' specified on material " // trim(to_str(mat % id)) + // "' specified on material " // trim(to_str(mat % id)) call fatal_error() end select end if @@ -1026,7 +1026,7 @@ contains ! Check to ensure material has at least one nuclide if (.not. associated(material_(i) % nuclides)) then message = "No nuclides specified on material " // & - trim(to_str(mat % id)) + trim(to_str(mat % id)) call fatal_error() end if @@ -1044,7 +1044,7 @@ contains ! Check for empty name on nuclide if (len_trim(nuc % name) == 0) then message = "No name specified on nuclide in material " // & - trim(to_str(mat % id)) + trim(to_str(mat % id)) call fatal_error() end if @@ -1052,7 +1052,7 @@ contains if (len_trim(nuc % xs) == 0) then if (default_xs == '') then message = "No cross section specified for nuclide in material " & - // trim(to_str(mat % id)) + // trim(to_str(mat % id)) call fatal_error() else nuc % xs = default_xs @@ -1066,7 +1066,7 @@ contains ! Check that this nuclide is listed in the cross_sections.xml file if (.not. dict_has_key(xs_listing_dict, name)) then message = "Could not find nuclide " // trim(name) // & - " in cross_sections.xml file!" + " in cross_sections.xml file!" call fatal_error() end if @@ -1074,7 +1074,7 @@ contains n = len_trim(name) if (name(n:n) /= 'c') then message = "Cross-section table " // trim(name) // & - " is not a continuous-energy neutron table." + " is not a continuous-energy neutron table." call fatal_error() end if @@ -1099,11 +1099,11 @@ contains ! weight percents were specified if (nuc % ao == ZERO .and. nuc % wo == ZERO) then message = "No atom or weight percent specified for nuclide " // & - trim(name) + trim(name) call fatal_error() elseif (nuc % ao /= ZERO .and. nuc % wo /= ZERO) then message = "Cannot specify both atom and weight percents for a & - &nuclide: " // trim(name) + &nuclide: " // trim(name) call fatal_error() end if @@ -1118,9 +1118,9 @@ contains ! Check to make sure either all atom percents or all weight percents are ! given if (.not. (all(mat % atom_density > ZERO) .or. & - all(mat % atom_density < ZERO))) then + all(mat % atom_density < ZERO))) then message = "Cannot mix atom and weight percents in material " // & - to_str(mat % id) + to_str(mat % id) call fatal_error() end if @@ -1141,7 +1141,7 @@ contains ! Check that this nuclide is listed in the cross_sections.xml file if (.not. dict_has_key(xs_listing_dict, name)) then message = "Could not find S(a,b) table " // trim(name) // & - " in cross_sections.xml file!" + " in cross_sections.xml file!" call fatal_error() end if mat % has_sab_table = .true. @@ -1205,6 +1205,7 @@ contains logical :: file_exists ! does tallies.xml file exist? character(MAX_LINE_LEN) :: filename character(MAX_WORD_LEN) :: word + character(MAX_WORD_LEN) :: scatter_name type(ListKeyValueCI), pointer :: key_list => null() type(TallyObject), pointer :: t => null() type(StructuredMesh), pointer :: m => null() @@ -1301,7 +1302,7 @@ contains ! Check that dimensions are all greater than zero if (any(mesh_(i) % dimension <= 0)) then message = "All entries on the element for a tally mesh & - &must be positive." + &must be positive." call fatal_error() end if @@ -1311,24 +1312,24 @@ contains ! Read mesh lower-left corner location if (m % n_dimension /= size(mesh_(i) % lower_left)) then message = "Number of entries on must be the same as & - &the number of entries on ." + &the number of entries on ." call fatal_error() end if m % lower_left = mesh_(i) % lower_left ! Make sure either upper-right or width was specified if (associated(mesh_(i) % upper_right) .and. & - associated(mesh_(i) % width)) then + associated(mesh_(i) % width)) then message = "Cannot specify both and on a & - &tally mesh." + &tally mesh." call fatal_error() end if ! Make sure either upper-right or width was specified if (.not. associated(mesh_(i) % upper_right) .and. & - .not. associated(mesh_(i) % width)) then + .not. associated(mesh_(i) % width)) then message = "Must specify either and on a & - &tally mesh." + &tally mesh." call fatal_error() end if @@ -1336,7 +1337,7 @@ contains ! Check to ensure width has same dimensions if (size(mesh_(i) % width) /= size(mesh_(i) % lower_left)) then message = "Number of entries on must be the same as the & - &number of entries on ." + &number of entries on ." call fatal_error() end if @@ -1354,14 +1355,14 @@ contains ! Check to ensure width has same dimensions if (size(mesh_(i) % upper_right) /= size(mesh_(i) % lower_left)) then message = "Number of entries on must be the same as & - &the number of entries on ." + &the number of entries on ." call fatal_error() end if ! Check that upper-right is above lower-left if (any(mesh_(i) % upper_right < mesh_(i) % lower_left)) then message = "The coordinates must be greater than the & - & coordinates on a tally mesh." + & coordinates on a tally mesh." call fatal_error() end if @@ -1410,7 +1411,7 @@ contains if (size(tally_(i) % filters) > 0) then message = "Tally filters should be specified with multiple & - &elements. Did you forget to change your element?" + &elements. Did you forget to change your element?" call fatal_error() end if @@ -1442,7 +1443,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('cellborn') @@ -1456,7 +1457,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('material') @@ -1470,7 +1471,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('universe') @@ -1484,7 +1485,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('surface') @@ -1498,7 +1499,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('mesh') @@ -1520,7 +1521,7 @@ contains m => meshes(i_mesh) else message = "Could not find mesh " // trim(to_str(id)) // & - " specified on tally " // trim(to_str(t % id)) + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if @@ -1544,7 +1545,7 @@ contains allocate(t % filters(j) % real_bins(n_words)) do k = 1, n_words t % filters(j) % real_bins(k) = str_to_real(& - tally_(i) % filter(j) % bins(k)) + tally_(i) % filter(j) % bins(k)) end do case ('energyout') @@ -1558,7 +1559,7 @@ contains allocate(t % filters(j) % real_bins(n_words)) do k = 1, n_words t % filters(j) % real_bins(k) = str_to_real(& - tally_(i) % filter(j) % bins(k)) + tally_(i) % filter(j) % bins(k)) end do end select @@ -1572,9 +1573,9 @@ contains ! Check that both cell and surface weren't specified if (t % find_filter(FILTER_CELL) > 0 .and. & - t % find_filter(FILTER_SURFACE) > 0) then + t % find_filter(FILTER_SURFACE) > 0) then message = "Cannot specify both cell and surface filters for tally " & - // trim(to_str(t % id)) + // trim(to_str(t % id)) call fatal_error() end if @@ -1593,7 +1594,7 @@ contains ! Set bins to 1, 2, 3, ..., n_nuclides_total, -1 t % nuclide_bins(1:n_nuclides_total) = & - (/ (j, j=1, n_nuclides_total) /) + (/ (j, j=1, n_nuclides_total) /) t % nuclide_bins(n_nuclides_total + 1) = -1 ! Set number of nuclide bins @@ -1621,7 +1622,7 @@ contains key_list => dict_keys(nuclide_dict) do while (associated(key_list)) if (starts_with(key_list % data % key, & - tally_(i) % nuclides(j))) then + tally_(i) % nuclides(j))) then word = key_list % data % key exit end if @@ -1633,8 +1634,8 @@ contains ! Check if no nuclide was found if (.not. associated(key_list)) then message = "Could not find the nuclide " // trim(& - tally_(i) % nuclides(j)) // " specified in tally " & - // trim(to_str(t % id)) // " in any material." + tally_(i) % nuclides(j)) // " specified in tally " & + // trim(to_str(t % id)) // " in any material." call fatal_error() end if deallocate(key_list) @@ -1647,7 +1648,7 @@ contains ! Check to make sure nuclide specified is in problem if (.not. dict_has_key(nuclide_dict, word)) then message = "The nuclide " // trim(word) // " from tally " // & - trim(to_str(t % id)) // " is not present in any material." + trim(to_str(t % id)) // " is not present in any material." call fatal_error() end if @@ -1678,26 +1679,28 @@ contains n_words = size(tally_(i) % scores) n_new = 0 do j = 1, n_words - call lower_case(tally_(i) % scores(j)) - select case (tally_(i) % scores(j)) - case ('scatter-p0') - n_new = n_new + 0 - case ('scatter-p1') - n_new = n_new + 1 - case ('scatter-p2') - n_new = n_new + 2 - case ('scatter-p3') - n_new = n_new + 3 - case ('scatter-p4') - n_new = n_new + 4 - case ('scatter-p5') - n_new = n_new + 5 - end select - end do + call lower_case(tally_(i) % scores(j)) + select case (tally_(i) % scores(j)) + case ('scatter-p0') + n_new = n_new + 0 + case ('scatter-p1') + n_new = n_new + 1 + case ('scatter-p2') + n_new = n_new + 2 + case ('scatter-p3') + n_new = n_new + 3 + case ('scatter-p4') + n_new = n_new + 4 + case ('scatter-p5') + n_new = n_new + 5 + end select + end do n_scores = n_words + n_new ! Allocate accordingly allocate(t % score_bins(n_scores)) + allocate(t % scatt_order(n_scores)) + t % scatt_order = 0 j = 0 do l = 1, n_words j = j + 1 @@ -1705,7 +1708,7 @@ contains case ('flux') ! Prohibit user from tallying flux for an individual nuclide if (.not. (t % n_nuclide_bins == 1 .and. & - t % nuclide_bins(1) == -1)) then + t % nuclide_bins(1) == -1)) then message = "Cannot tally flux for an individual nuclide." call fatal_error() end if @@ -1719,7 +1722,7 @@ contains t % score_bins(j) = SCORE_TOTAL if (t % find_filter(FILTER_ENERGYOUT) > 0) then message = "Cannot tally total reaction rate with an & - &outgoing energy filter." + &outgoing energy filter." call fatal_error() end if case ('scatter') @@ -1729,92 +1732,118 @@ contains ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG - case ('scatter-0') !does the same as 'scatter', for convenience + case ('scatter-0') + ! Does the same as 'scatter', for convenience t % score_bins(j) = SCORE_SCATTER case ('scatter-1') - t % score_bins(j) = SCORE_SCATTER_1 + t % score_bins(j) = SCORE_SCATTER_N + t % scatt_order(j) = 1 ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('scatter-2') - t % score_bins(j) = SCORE_SCATTER_2 + t % score_bins(j) = SCORE_SCATTER_N + t % scatt_order(j) = 2 ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('scatter-3') - t % score_bins(j) = SCORE_SCATTER_3 + t % score_bins(j) = SCORE_SCATTER_N + t % scatt_order(j) = 3 ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('scatter-4') - t % score_bins(j) = SCORE_SCATTER_4 + t % score_bins(j) = SCORE_SCATTER_N + t % scatt_order(j) = 4 ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('scatter-5') - t % score_bins(j) = SCORE_SCATTER_5 + t % score_bins(j) = SCORE_SCATTER_N + t % scatt_order(j) = 5 ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('scatter-p0') - ! Same as a scatter-0 or scatter - t % score_bins(j) = SCORE_SCATTER - case ('scatter-p1') - ! Setup each of P0 and P1 - t % estimator = ESTIMATOR_ANALOG - t % score_bins(j) = SCORE_SCATTER - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_1 - - case ('scatter-p2') - ! Setup each of P0:P2 - t % estimator = ESTIMATOR_ANALOG - t % score_bins(j) = SCORE_SCATTER - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_1 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_2 - - case ('scatter-p3') - ! Setup each of P0:P3 - t % estimator = ESTIMATOR_ANALOG - t % score_bins(j) = SCORE_SCATTER - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_1 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_2 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_3 - - case ('scatter-p4') - ! Setup each of P0:P4 - t % estimator = ESTIMATOR_ANALOG - t % score_bins(j) = SCORE_SCATTER - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_1 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_2 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_3 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_4 - - case ('scatter-p5') - ! Setup each of P0:P5 - t % estimator = ESTIMATOR_ANALOG - t % score_bins(j) = SCORE_SCATTER - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_1 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_2 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_3 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_4 - j = j + 1 - t % score_bins(j) = SCORE_SCATTER_5 - + ! Same as a scatter-0 or scatter + t % score_bins(j) = SCORE_SCATTER + case ('scatter-p1') + ! Setup each of P0 and P1 + t % estimator = ESTIMATOR_ANALOG + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 1 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 1 + + case ('scatter-p2') + ! Setup each of P0:P2 + t % estimator = ESTIMATOR_ANALOG + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 2 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 2 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 2 + + case ('scatter-p3') + ! Setup each of P0:P3 + t % estimator = ESTIMATOR_ANALOG + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 3 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 3 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 3 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 3 + + case ('scatter-p4') + ! Setup each of P0:P4 + t % estimator = ESTIMATOR_ANALOG + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 4 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 4 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 4 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 4 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 4 + + case ('scatter-p5') + ! Setup each of P0:P5 + t % estimator = ESTIMATOR_ANALOG + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 5 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 5 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 5 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 5 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 5 + j = j + 1 + t % score_bins(j) = SCORE_SCATTER_PN + t % scatt_order(j) = 5 + case('transport') t % score_bins(j) = SCORE_TRANSPORT @@ -1849,14 +1878,14 @@ contains t % score_bins(j) = SCORE_ABSORPTION if (t % find_filter(FILTER_ENERGYOUT) > 0) then message = "Cannot tally absorption rate with an outgoing & - &energy filter." + &energy filter." call fatal_error() end if case ('fission') t % score_bins(j) = SCORE_FISSION if (t % find_filter(FILTER_ENERGYOUT) > 0) then message = "Cannot tally fission rate with an outgoing & - &energy filter." + &energy filter." call fatal_error() end if case ('nu-fission') @@ -1869,8 +1898,8 @@ contains ! for this tally if (n_words > 1) then message = "Cannot tally other scoring functions in the same & - &tally as surface currents. Separate other scoring & - &functions into a distinct tally." + &tally as surface currents. Separate other scoring & + &functions into a distinct tally." call fatal_error() end if @@ -1902,13 +1931,13 @@ contains t % filters(t % n_filters) % type = FILTER_SURFACE t % filters(t % n_filters) % n_bins = 2 * m % n_dimension allocate(t % filters(t % n_filters) % int_bins(& - 2 * m % n_dimension)) + 2 * m % n_dimension)) if (m % n_dimension == 2) then t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, & - OUT_RIGHT, IN_FRONT, OUT_FRONT /) + OUT_RIGHT, IN_FRONT, OUT_FRONT /) elseif (m % n_dimension == 3) then t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, & - OUT_RIGHT, IN_FRONT, OUT_FRONT, IN_TOP, OUT_TOP /) + OUT_RIGHT, IN_FRONT, OUT_FRONT, IN_TOP, OUT_TOP /) end if t % find_filter(FILTER_SURFACE) = t % n_filters @@ -1917,14 +1946,15 @@ contains case default message = "Unknown scoring function: " // & - trim(tally_(i) % scores(l)) + trim(tally_(i) % scores(l)) call fatal_error() end select end do t % n_score_bins = n_scores + t % n_nonPN_score_bins = n_words else message = "No specified on tally " // trim(to_str(t % id)) & - // "." + // "." call fatal_error() end if @@ -1942,7 +1972,7 @@ contains ! tally needs post-collision information if (t % estimator == ESTIMATOR_ANALOG) then message = "Cannot use track-length estimator for tally " & - // to_str(t % id) + // to_str(t % id) call fatal_error() end if @@ -1951,7 +1981,7 @@ contains case default message = "Invalid estimator '" // trim(tally_(i) % estimator) & - // "' on tally " // to_str(t % id) + // "' on tally " // to_str(t % id) call fatal_error() end select end if @@ -2055,7 +2085,7 @@ contains ! Set output file path pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // & - "_" // trim(plot_(i) % filename) // ".ppm" + "_" // trim(plot_(i) % filename) // ".ppm" ! Copy plot pixel size if (size(plot_(i) % pixels) == 2) then @@ -2071,7 +2101,7 @@ contains pl % not_found % rgb = plot_(i) % background else message = "Bad background RGB " & - // "in plot " // trim(to_str(pl % id)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if else @@ -2086,7 +2116,7 @@ contains ! pl % type = PLOT_TYPE_POINTS case default message = "Unsupported plot type '" // plot_(i) % type & - // "' in plot " // trim(to_str(pl % id)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -2100,7 +2130,7 @@ contains pl % basis = PLOT_BASIS_YZ case default message = "Unsupported plot basis '" // plot_(i) % basis & - // "' in plot " // trim(to_str(pl % id)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -2109,7 +2139,7 @@ contains pl % origin = plot_(i) % origin else message = "Origin must be length 3 " & - // "in plot " // trim(to_str(pl % id)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2121,7 +2151,7 @@ contains pl % width(2) = plot_(i) % width(2) else message = "Bad plot width " & - // "in plot " // trim(to_str(pl % id)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2149,7 +2179,7 @@ contains case default message = "Unsupported plot color type '" // plot_(i) % color & - // "' in plot " // trim(to_str(pl % id)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -2159,7 +2189,7 @@ contains do j = 1, n_cols if (size(plot_(i) % col_spec_(j) % rgb) /= 3) then message = "Bad RGB " & - // "in plot " // trim(to_str(pl % id)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2171,7 +2201,7 @@ contains pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb else message = "Could not find cell " // trim(to_str(col_id)) // & - " specified in plot " // trim(to_str(pl % id)) + " specified in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2181,7 +2211,7 @@ contains pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb else message = "Could not find material " // trim(to_str(col_id)) // & - " specified in plot " // trim(to_str(pl % id)) + " specified in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2193,7 +2223,7 @@ contains if (associated(plot_(i) % mask_)) then if (size(plot_(i) % mask_) > 1) then message = "Mutliple masks" // & - " specified in plot " // trim(to_str(pl % id)) + " specified in plot " // trim(to_str(pl % id)) call fatal_error() else if (.not. size(plot_(i) % mask_) == 0) then do j=1,size(pl % colors) @@ -2228,10 +2258,10 @@ contains ! Check if cross_sections.xml exists inquire(FILE=path_cross_sections, EXIST=file_exists) if (.not. file_exists) then - ! Could not find cross_sections.xml file - message = "Cross sections XML file '" // trim(path_cross_sections) // & + ! Could not find cross_sections.xml file + message = "Cross sections XML file '" // trim(path_cross_sections) // & "' does not exist!" - call fatal_error() + call fatal_error() end if message = "Reading cross sections XML file..." @@ -2247,25 +2277,25 @@ contains call read_xml_file_cross_sections_t(path_cross_sections) if (len_trim(directory_) > 0) then - ! Copy directory information if present - directory = trim(directory_) + ! Copy directory information if present + directory = trim(directory_) else - ! If no directory is listed in cross_sections.xml, by default select the - ! directory in which the cross_sections.xml file resides - i = index(path_cross_sections, "/", BACK=.true.) - directory = path_cross_sections(1:i) + ! If no directory is listed in cross_sections.xml, by default select the + ! directory in which the cross_sections.xml file resides + i = index(path_cross_sections, "/", BACK=.true.) + directory = path_cross_sections(1:i) end if ! determine whether binary/ascii if (filetype_ == 'ascii') then - filetype = ASCII + filetype = ASCII elseif (filetype_ == 'binary') then - filetype = BINARY + filetype = BINARY elseif (len_trim(filetype_) == 0) then - filetype = ASCII + filetype = ASCII else - message = "Unknown filetype in cross_sections.xml: " // trim(filetype_) - call fatal_error() + message = "Unknown filetype in cross_sections.xml: " // trim(filetype_) + call fatal_error() end if ! copy default record length and entries for binary files @@ -2274,57 +2304,57 @@ contains ! Allocate xs_listings array if (.not. associated(ace_tables_)) then - message = "No ACE table listings present in cross_sections.xml file!" - call fatal_error() + message = "No ACE table listings present in cross_sections.xml file!" + call fatal_error() else - n_listings = size(ace_tables_) - allocate(xs_listings(n_listings)) + n_listings = size(ace_tables_) + allocate(xs_listings(n_listings)) end if do i = 1, n_listings - listing => xs_listings(i) + listing => xs_listings(i) - ! copy a number of attributes - listing % name = trim(ace_tables_(i) % name) - listing % alias = trim(ace_tables_(i) % alias) - listing % zaid = ace_tables_(i) % zaid - listing % awr = ace_tables_(i) % awr - listing % kT = ace_tables_(i) % temperature - listing % location = ace_tables_(i) % location + ! copy a number of attributes + listing % name = trim(ace_tables_(i) % name) + listing % alias = trim(ace_tables_(i) % alias) + listing % zaid = ace_tables_(i) % zaid + listing % awr = ace_tables_(i) % awr + listing % kT = ace_tables_(i) % temperature + listing % location = ace_tables_(i) % location - ! determine type of cross section - if (ends_with(listing % name, 'c')) then + ! determine type of cross section + if (ends_with(listing % name, 'c')) then listing % type = ACE_NEUTRON - elseif (ends_with(listing % name, 't')) then + elseif (ends_with(listing % name, 't')) then listing % type = ACE_THERMAL - end if + end if - ! set filetype, record length, and number of entries - listing % filetype = filetype - listing % recl = recl - listing % entries = entries + ! set filetype, record length, and number of entries + listing % filetype = filetype + listing % recl = recl + listing % entries = entries - ! determine metastable state - if (ace_tables_(i) % metastable == 0) then + ! determine metastable state + if (ace_tables_(i) % metastable == 0) then listing % metastable = .false. - else + else listing % metastable = .true. - end if + end if - ! determine path of cross section table - if (starts_with(ace_tables_(i) % path, '/')) then + ! determine path of cross section table + if (starts_with(ace_tables_(i) % path, '/')) then listing % path = ace_tables_(i) % path - else + else if (ends_with(directory,'/')) then - listing % path = trim(directory) // trim(ace_tables_(i) % path) + listing % path = trim(directory) // trim(ace_tables_(i) % path) else - listing % path = trim(directory) // '/' // trim(ace_tables_(i) % path) + listing % path = trim(directory) // '/' // trim(ace_tables_(i) % path) end if - end if + end if - ! create dictionary entry for both name and alias - call dict_add_key(xs_listing_dict, listing % name, i) - call dict_add_key(xs_listing_dict, listing % alias, i) + ! create dictionary entry for both name and alias + call dict_add_key(xs_listing_dict, listing % name, i) + call dict_add_key(xs_listing_dict, listing % alias, i) end do end subroutine read_cross_sections_xml diff --git a/src/output.F90 b/src/output.F90 index 04cfcb8a9a..2771913a43 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -32,35 +32,35 @@ contains subroutine title() write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') & - ' .d88888b. 888b d888 .d8888b.', & - ' d88P" "Y88b 8888b d8888 d88P Y88b', & - ' 888 888 88888b.d88888 888 888', & - ' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', & - ' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', & - ' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', & - ' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', & - ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', & - '__________________888______________________________________________________', & - ' 888', & - ' 888' + ' .d88888b. 888b d888 .d8888b.', & + ' d88P" "Y88b 8888b d8888 d88P Y88b', & + ' 888 888 88888b.d88888 888 888', & + ' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', & + ' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', & + ' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', & + ' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', & + ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', & + '__________________888______________________________________________________', & + ' 888', & + ' 888' ! Write version information write(UNIT=OUTPUT_UNIT, FMT=*) & - ' Developed At: Massachusetts Institute of Technology' + ' Developed At: Massachusetts Institute of Technology' write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') & - VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE + VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE #ifdef GIT_SHA1 write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",6X,A)') GIT_SHA1 #endif ! Write the date and time write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A)') & - time_stamp() + time_stamp() #ifdef MPI ! Write number of processors write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",1X,A)') & - trim(to_str(n_procs)) + trim(to_str(n_procs)) #endif end subroutine title @@ -77,7 +77,7 @@ contains call date_and_time(DATE=date_, TIME=time_) current_time = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8) // & - " " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) + " " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) end function time_stamp @@ -125,13 +125,13 @@ contains select case (header_level) case (1) write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), & - repeat('=', n) // '> ' // trim(line) // ' <' // & - repeat('=', m), repeat('=', 75) + repeat('=', n) // '> ' // trim(line) // ' <' // & + repeat('=', m), repeat('=', 75) case (2) write(UNIT=unit_, FMT='(/2(1X,A/))') trim(line), repeat('-', 75) case (3) write(UNIT=unit_, FMT='(/1X,A/)') repeat('=', n) // '> ' // & - trim(line) // ' <' // repeat('=', m) + trim(line) // ' <' // repeat('=', m) end select end subroutine header @@ -145,11 +145,11 @@ contains if (master) then write(UNIT=OUTPUT_UNIT, FMT='(1X,A,1X,I1,".",I1,".",I1)') & - "OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE + "OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE write(UNIT=OUTPUT_UNIT, FMT=*) "Copyright (c) 2011-2012 & - &Massachusetts Institute of Technology" + &Massachusetts Institute of Technology" write(UNIT=OUTPUT_UNIT, FMT=*) "MIT/X license at & - &" + &" end if end subroutine print_version @@ -207,7 +207,7 @@ contains else ! Determine last space in current line i_end = i_start + index(message(i_start+1:i_start+line_wrap), & - ' ', BACK=.true.) + ' ', BACK=.true.) ! Write up to last space write(ou, fmt='(1X,A)') message(i_start+1:i_end-1) @@ -270,8 +270,8 @@ contains l => lattices(coord % lattice) write(ou,*) ' Lattice = ' // trim(to_str(l % id)) write(ou,*) ' Lattice position = (' // trim(to_str(& - p % coord % lattice_x)) // ',' // trim(to_str(& - p % coord % lattice_y)) // ')' + p % coord % lattice_x)) // ',' // trim(to_str(& + p % coord % lattice_y)) // ')' end if ! Print local coordinates @@ -390,7 +390,7 @@ contains case default index_surf = abs(c % surfaces(i)) string = trim(string) // ' ' // to_str(sign(& - surfaces(index_surf) % id, c % surfaces(i))) + surfaces(index_surf) % id, c % surfaces(i))) end select end do write(unit_,*) ' Surface Specification:' // trim(string) @@ -595,7 +595,7 @@ contains ! Write total atom density in atom/b-cm write(unit_,*) ' Atom Density = ' // trim(to_str(mat % density)) & - // ' atom/b-cm' + // ' atom/b-cm' ! Write atom density for each nuclide in material write(unit_,*) ' Nuclides:' @@ -603,7 +603,7 @@ contains nuc => nuclides(mat % nuclide(i)) density = mat % atom_density(i) string = ' ' // trim(nuc % name) // ' = ' // & - trim(to_str(density)) // ' atom/b-cm' + trim(to_str(density)) // ' atom/b-cm' write(unit_,*) trim(string) end do @@ -628,7 +628,9 @@ contains integer :: j ! index in filters array integer :: id ! user-specified id integer :: unit_ ! unit to write to + integer :: n ! scattering order to include in name character(MAX_LINE_LEN) :: string + character(MAX_WORD_LEN) :: pn_string type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() type(Universe), pointer :: u => null() @@ -740,7 +742,7 @@ contains string = "" do i = 1, t % filters(j) % n_bins + 1 string = trim(string) // ' ' // trim(to_str(& - t % filters(j) % real_bins(i))) + t % filters(j) % real_bins(i))) end do write(unit_,*) ' Incoming Energy Bins:' // trim(string) end if @@ -751,7 +753,7 @@ contains string = "" do i = 1, t % filters(j) % n_bins + 1 string = trim(string) // ' ' // trim(to_str(& - t % filters(j) % real_bins(i))) + t % filters(j) % real_bins(i))) end do write(unit_,*) ' Outgoing Energy Bins:' // trim(string) end if @@ -763,18 +765,20 @@ contains write(unit_,fmt='(A)',advance='no') ' total' else write(unit_,fmt='(A)',advance='no') ' ' // trim(adjustl(& - nuclides(t % nuclide_bins(i)) % name)) + nuclides(t % nuclide_bins(i)) % name)) end if if (mod(i,4) == 0 .and. i /= t % n_nuclide_bins) & - write(unit_,'(/18X)',advance='no') + write(unit_,'(/18X)',advance='no') end do write(unit_,*) ! Write score bins - string = "" - do i = 1, t % n_score_bins - select case (t % score_bins(i)) + string = "" + j = 0 + do i = 1, t % n_nonPN_score_bins + j = j + 1 + select case (t % score_bins(j)) case (SCORE_FLUX) string = trim(string) // ' flux' case (SCORE_TOTAL) @@ -783,16 +787,17 @@ contains string = trim(string) // ' scatter' case (SCORE_NU_SCATTER) string = trim(string) // ' nu-scatter' - case (SCORE_SCATTER_1) - string = trim(string) // ' scatter-1' - case (SCORE_SCATTER_2) - string = trim(string) // ' scatter-2' - case (SCORE_SCATTER_3) - string = trim(string) // ' scatter-3' - case (SCORE_SCATTER_4) - string = trim(string) // ' scatter-4' - case (SCORE_SCATTER_5) - string = trim(string) // ' scatter-5' + case (SCORE_SCATTER_N) + pn_string = ' scatter-' // trim(to_str(t % scatt_order(j))) + string = trim(string) // pn_string + case (SCORE_SCATTER_PN) + pn_string = ' scatter' + string = trim(string) // pn_string + do n = 1, t % scatt_order(j) + pn_string = ' scatter-' // trim(to_str(n)) + string = trim(string) // pn_string + end do + j = j + n - 1 case (SCORE_TRANSPORT) string = trim(string) // ' transport' case (SCORE_DIFFUSION) @@ -933,8 +938,8 @@ contains 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 + 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 @@ -971,11 +976,11 @@ contains write(unit_,*) ' Memory Requirements' write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes' write(unit_,*) ' Secondary angle distributions = ' // & - trim(to_str(size_angle_total)) // ' bytes' + trim(to_str(size_angle_total)) // ' bytes' write(unit_,*) ' Secondary energy distributions = ' // & - trim(to_str(size_energy_total)) // ' bytes' + trim(to_str(size_energy_total)) // ' bytes' write(unit_,*) ' Probability Tables = ' // & - trim(to_str(size_urr)) // ' bytes' + trim(to_str(size_urr)) // ' bytes' write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes' ! Blank line at end of nuclide @@ -1011,28 +1016,28 @@ contains ! Inelastic data write(unit_,*) ' # of Incoming Energies (Inelastic) = ' // & - trim(to_str(sab % n_inelastic_e_in)) + trim(to_str(sab % n_inelastic_e_in)) write(unit_,*) ' # of Outgoing Energies (Inelastic) = ' // & - trim(to_str(sab % n_inelastic_e_out)) + trim(to_str(sab % n_inelastic_e_out)) write(unit_,*) ' # of Outgoing Angles (Inelastic) = ' // & - trim(to_str(sab % n_inelastic_mu)) + trim(to_str(sab % n_inelastic_mu)) write(unit_,*) ' Threshold for Inelastic = ' // & - trim(to_str(sab % threshold_inelastic)) + trim(to_str(sab % threshold_inelastic)) ! Elastic data if (sab % n_elastic_e_in > 0) then write(unit_,*) ' # of Incoming Energies (Elastic) = ' // & - trim(to_str(sab % n_elastic_e_in)) + trim(to_str(sab % n_elastic_e_in)) write(unit_,*) ' # of Outgoing Angles (Elastic) = ' // & - trim(to_str(sab % n_elastic_mu)) + trim(to_str(sab % n_elastic_mu)) write(unit_,*) ' Threshold for Elastic = ' // & - trim(to_str(sab % threshold_elastic)) + trim(to_str(sab % threshold_elastic)) end if ! Determine memory used by S(a,b) table and write out size_sab = 8 * (sab % n_inelastic_e_in * (2 + sab % n_inelastic_e_out * & - (1 + sab % n_inelastic_mu)) + sab % n_elastic_e_in * & - (2 + sab % n_elastic_mu)) + (1 + sab % n_inelastic_mu)) + sab % n_elastic_e_in * & + (2 + sab % n_elastic_mu)) write(unit_,*) ' Memory Used = ' // trim(to_str(size_sab)) // ' bytes' ! Blank line at end @@ -1060,19 +1065,19 @@ contains call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1) write(UNIT=UNIT_SUMMARY, FMT=*) & - "Copyright: 2011-2012 Massachusetts Institute of Technology" + "Copyright: 2011-2012 Massachusetts Institute of Technology" write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') & - "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE + "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE #ifdef GIT_SHA1 write(UNIT=UNIT_SUMMARY, FMT='(1X,"Git SHA1:",6X,A)') GIT_SHA1 #endif write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A)') & - time_stamp() + time_stamp() ! Write information on number of processors #ifdef MPI write(UNIT=UNIT_SUMMARY, FMT='(1X,"MPI Processes:",1X,A)') & - trim(to_str(n_procs)) + trim(to_str(n_procs)) #endif ! Display problem summary @@ -1117,7 +1122,7 @@ contains call header("UNIONIZED ENERGY GRID", unit=UNIT_SUMMARY) write(UNIT_SUMMARY,*) "Points on energy grid: " // trim(to_str(n_grid)) write(UNIT_SUMMARY,*) "Extra storage required: " // trim(to_str(& - n_grid*n_nuclides_total*4)) // " bytes" + n_grid*n_nuclides_total*4)) // " bytes" ! print summary of variance reduction call header("VARIANCE REDUCTION", unit=UNIT_SUMMARY) @@ -1227,16 +1232,16 @@ contains ! write out information batch and option independent output write(UNIT=OUTPUT_UNIT, FMT='(2X,I5)', ADVANCE='NO') current_batch write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') & - k_batch(current_batch) + k_batch(current_batch) ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - entropy(current_batch) + entropy(current_batch) ! write out accumulated k-effective if after first active batch if (current_batch > n_inactive + 1) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') & - keff, keff_std + keff, keff_std else write(UNIT=OUTPUT_UNIT, FMT='(23X)', ADVANCE='NO') end if @@ -1244,11 +1249,11 @@ contains ! write out cmfd keff if it is active if (cmfd_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - cmfd % keff + cmfd % keff ! write out cmfd entopy if (cmfd_on .and. entropy_on) write(UNIT=OUTPUT_UNIT, & - FMT='(3X, F8.5)', ADVANCE='NO') cmfd % entropy + FMT='(3X, F8.5)', ADVANCE='NO') cmfd % entropy ! next line write(UNIT=OUTPUT_UNIT, FMT=*) @@ -1275,18 +1280,18 @@ contains ! Write plotting origin write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // & - " " // trim(to_str(pl % origin(2))) // " " // & - trim(to_str(pl % origin(3))) + " " // trim(to_str(pl % origin(2))) // " " // & + trim(to_str(pl % origin(3))) ! Write plotting width if (pl % type == PLOT_TYPE_SLICE) then write(ou,100) "Width:", trim(to_str(pl % width(1))) // & - " " // trim(to_str(pl % width(2))) + " " // trim(to_str(pl % width(2))) write(ou,100) "Coloring:", trim(to_str(pl % color_by)) write(ou,100) "Basis:", trim(to_str(pl % basis)) write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // & - trim(to_str(pl % pixels(2))) + trim(to_str(pl % pixels(2))) end if write(ou,*) @@ -1319,7 +1324,7 @@ contains write(ou,100) " Reading cross sections", time_read_xs % elapsed write(ou,100) " Unionizing energy grid", time_unionize % elapsed write(ou,100) "Total time in simulation", time_inactive % elapsed + & - time_active % elapsed + time_active % elapsed write(ou,100) " Time in transport only", time_transport % elapsed if(cmfd_run) write(ou,100) "Total CMFD time", time_cmfd % elapsed write(ou,100) " Time in inactive batches", time_inactive % elapsed @@ -1333,14 +1338,14 @@ contains if (restart_run) then total_particles = n_particles * (n_batches - & - restart_batch) * gen_per_batch + restart_batch) * gen_per_batch else total_particles = n_particles * n_batches * gen_per_batch end if ! display calculate rate speed = real(total_particles) / (time_inactive % elapsed + & - time_active % elapsed) + time_active % elapsed) string = to_str(speed) write(ou,101) "Calculation Rate", trim(string) @@ -1358,13 +1363,13 @@ contains ! write global tallies write(ou,102) "k-effective (Analog)", global_tallies(K_ANALOG) % sum, & - global_tallies(K_ANALOG) % sum_sq + global_tallies(K_ANALOG) % sum_sq write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) % sum, & - global_tallies(K_COLLISION) % sum_sq + global_tallies(K_COLLISION) % sum_sq write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum, & - global_tallies(K_TRACKLENGTH) % sum_sq + global_tallies(K_TRACKLENGTH) % sum_sq write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, & - global_tallies(LEAKAGE) % sum_sq + global_tallies(LEAKAGE) % sum_sq write(ou,*) ! format for write statements @@ -1385,17 +1390,21 @@ contains integer :: j ! level in tally hierarchy integer :: k ! loop index for scoring bins integer :: n ! loop index for nuclides + integer :: l ! loop index for nonPN scores integer :: type ! type of tally filter integer :: indent ! number of spaces to preceed output integer :: filter_index ! index in scores array for filters integer :: score_index ! scoring bin index integer :: i_nuclide ! index in nuclides array integer :: i_listing ! index in xs_listings array + integer :: nOrder ! loop index for scattering orders real(8) :: t_value ! t-values for confidence intervals real(8) :: alpha ! significance level for CI character(MAX_FILE_LEN) :: filename ! name of output file character(15) :: filter_name(N_FILTER_TYPES) ! names of tally filters - character(27) :: score_name(N_SCORE_TYPES) ! names of scoring function + character(27) :: score_names(N_SCORE_TYPES) ! names of scoring function + character(27) :: score_name ! names of scoring function + ! to be applied at write-time type(TallyObject), pointer :: t ! Skip if there are no tallies @@ -1412,25 +1421,22 @@ contains filter_name(FILTER_ENERGYOUT) = "Outgoing Energy" ! Initialize names for scores - score_name(abs(SCORE_FLUX)) = "Flux" - score_name(abs(SCORE_TOTAL)) = "Total Reaction Rate" - score_name(abs(SCORE_SCATTER)) = "Scattering Rate" - score_name(abs(SCORE_NU_SCATTER)) = "Scattering Production Rate" - score_name(abs(SCORE_SCATTER_1)) = "First Scattering Moment" - score_name(abs(SCORE_SCATTER_2)) = "Second Scattering Moment" - score_name(abs(SCORE_SCATTER_3)) = "Third Scattering Moment" - score_name(abs(SCORE_SCATTER_4)) = "Fourth Scattering Moment" - score_name(abs(SCORE_SCATTER_5)) = "Fifth Scattering Moment" - score_name(abs(SCORE_TRANSPORT)) = "Transport Rate" - score_name(abs(SCORE_DIFFUSION)) = "Diffusion Coefficient" - score_name(abs(SCORE_N_1N)) = "(n,1n) Rate" - score_name(abs(SCORE_N_2N)) = "(n,2n) Rate" - score_name(abs(SCORE_N_3N)) = "(n,3n) Rate" - score_name(abs(SCORE_N_4N)) = "(n,4n) Rate" - score_name(abs(SCORE_ABSORPTION)) = "Absorption Rate" - score_name(abs(SCORE_FISSION)) = "Fission Rate" - score_name(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate" - score_name(abs(SCORE_EVENTS)) = "Events" + score_names(abs(SCORE_FLUX)) = "Flux" + score_names(abs(SCORE_TOTAL)) = "Total Reaction Rate" + score_names(abs(SCORE_SCATTER)) = "Scattering Rate" + score_names(abs(SCORE_NU_SCATTER)) = "Scattering Production Rate" + score_names(abs(SCORE_SCATTER_N)) = "" + score_names(abs(SCORE_SCATTER_PN)) = "" + score_names(abs(SCORE_TRANSPORT)) = "Transport Rate" + score_names(abs(SCORE_DIFFUSION)) = "Diffusion Coefficient" + score_names(abs(SCORE_N_1N)) = "(n,1n) Rate" + score_names(abs(SCORE_N_2N)) = "(n,2n) Rate" + score_names(abs(SCORE_N_3N)) = "(n,3n) Rate" + score_names(abs(SCORE_N_4N)) = "(n,4n) Rate" + score_names(abs(SCORE_ABSORPTION)) = "Absorption Rate" + score_names(abs(SCORE_FISSION)) = "Fission Rate" + score_names(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate" + score_names(abs(SCORE_EVENTS)) = "Events" ! Create filename for tally output if (run_mode == MODE_TALLIES) then @@ -1468,10 +1474,10 @@ contains ! Write header block if (t % label == "") then call header("TALLY " // trim(to_str(t % id)), unit=UNIT_TALLY, & - level=3) + level=3) else call header("TALLY " // trim(to_str(t % id)) // ": " & - // trim(t % label), unit=UNIT_TALLY, level=3) + // trim(t % label), unit=UNIT_TALLY, level=3) endif ! Handle surface current tallies separately @@ -1520,7 +1526,7 @@ contains ! Print current filter information type = t % filters(j) % type write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - trim(filter_name(type)), trim(get_label(t, j)) + trim(filter_name(type)), trim(get_label(t, j)) indent = indent + 2 j = j + 1 end if @@ -1531,7 +1537,7 @@ contains if (t % n_filters > 0) then type = t % filters(j) % type write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - trim(filter_name(type)), trim(get_label(t, j)) + trim(filter_name(type)), trim(get_label(t, j)) end if ! Determine scoring index for this bin combination -- note that unlike @@ -1543,7 +1549,7 @@ contains else filter_index = 1 end if - + ! Write scores for this filter bin combination score_index = 0 if (t % n_filters > 0) indent = indent + 2 @@ -1552,20 +1558,52 @@ contains i_nuclide = t % nuclide_bins(n) if (i_nuclide == -1) then write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - "Total Material" + "Total Material" else i_listing = nuclides(i_nuclide) % listing write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - trim(xs_listings(i_listing) % alias) + trim(xs_listings(i_listing) % alias) end if indent = indent + 2 - do k = 1, t % n_score_bins + k = 0 + do l = 1, t % n_nonPN_score_bins + k = k + 1 score_index = score_index + 1 - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), score_name(abs(t % score_bins(k))), & - to_str(t % scores(score_index,filter_index) % sum), & - trim(to_str(t % scores(score_index,filter_index) % sum_sq)) + if (t % score_bins(k) == SCORE_SCATTER_N) then + if (t % scatt_order(k) == 0) then + score_name = "Scattering Rate" + else + score_name = 'Order-' // trim(to_str(t % scatt_order(k))) // & + ' Scattering Moment' + end if + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % scores(score_index,filter_index) % sum), & + trim(to_str(t % scores(score_index,filter_index) % sum_sq)) + else if (t % score_bins(k) == SCORE_SCATTER_PN) then + score_name = "Scattering Rate" + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % scores(score_index,filter_index) % sum), & + trim(to_str(t % scores(score_index,filter_index) % sum_sq)) + do nOrder = 1, t % scatt_order(k) + score_index = score_index + 1 + score_name = 'Order-' // trim(to_str(nOrder)) // & + ' Scattering Moment' + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % scores(score_index,filter_index) % sum), & + trim(to_str(t % scores(score_index,filter_index) % sum_sq)) + end do + k = k + nOrder - 1 + else + score_name = score_names(abs(t % score_bins(k))) + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % scores(score_index,filter_index) % sum), & + trim(to_str(t % scores(score_index,filter_index) % sum_sq)) + end if end do indent = indent - 2 @@ -1642,110 +1680,110 @@ contains ! Write incoming energy bin write(UNIT=UNIT_TALLY, FMT='(3X,A,1X,A)') & - "Incoming Energy", trim(get_label(t, i_filter_ein)) + "Incoming Energy", trim(get_label(t, i_filter_ein)) end if ! Left Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Left", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Left", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Left", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Left", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) ! Right Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Right", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Right", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Right", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Right", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) ! Back Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Back", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Back", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_FRONT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Back", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Back", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) ! Front Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Front", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Front", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_FRONT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Front", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Front", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) ! Bottom Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_TOP filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Bottom", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Bottom", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Bottom", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Bottom", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) ! Top Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_TOP filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Top", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Top", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Top", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Top", & + to_str(t % scores(1,filter_index) % sum), & + trim(to_str(t % scores(1,filter_index) % sum_sq)) end do end do @@ -1793,10 +1831,10 @@ contains call bin_to_mesh_indices(m, bin, ijk) if (m % n_dimension == 2) then label = "Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ")" + trim(to_str(ijk(2))) // ")" elseif (m % n_dimension == 3) then label = "Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" + trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" end if case (FILTER_ENERGYIN, FILTER_ENERGYOUT) E0 = t % filters(i_filter) % real_bins(bin) diff --git a/src/tally.F90 b/src/tally.F90 index daf3973e21..d7f05a5fb2 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -6,7 +6,7 @@ module tally use global use math, only: t_percentile, calc_pn use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices, & - mesh_indices_to_bin, mesh_intersects + mesh_indices_to_bin, mesh_intersects use mesh_header, only: StructuredMesh use output, only: header use search, only: binary_search @@ -34,6 +34,9 @@ contains integer :: j ! loop index for scoring bins integer :: k ! loop index for nuclide bins + integer :: n ! loop index for scattering order + integer :: l ! scoring bin loop index, allowing for changing + ! position during the loop integer :: filter_index ! single index for single bin integer :: score_bin ! scoring bin, e.g. SCORE_FLUX integer :: i_nuclide ! index in nuclides array @@ -119,7 +122,9 @@ contains end if ! Determine score for each bin - SCORE_LOOP: do j = 1, t % n_score_bins + j = 0 + SCORE_LOOP: do l = 1, t % n_nonPN_score_bins + j = j + 1 ! determine what type of score bin score_bin = t % score_bins(j) @@ -160,52 +165,36 @@ contains ! reaction with neutrons in the exit channel score = wgt - - case (SCORE_SCATTER_1) + + case (SCORE_SCATTER_N) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - ! The first scattering moment can be determined by using the - ! rate of scattering reactions multiplied by the cosine of the - ! change in neutron's angle due to the collision + ! Find the scattering order for a singly requested moment, and + ! store its moment contribution. - score = last_wgt * mu + if (t % scatt_order(j) == 1) then + score = last_wgt * mu ! avoid function call overhead + else + score = last_wgt * calc_pn(t % scatt_order(j), mu) + endif - case (SCORE_SCATTER_2) + case (SCORE_SCATTER_PN) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP + score_index = score_index - 1 + ! Find the scattering order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % scatt_order(j) + ! determine scoring bin index + score_index = score_index + 1 + ! get the score and tally it + score = last_wgt * calc_pn(n, mu) + if (n < t % scatt_order(j)) & + call add_to_score(t % scores(score_index, filter_index), score) + end do + j = j + t % scatt_order(j) - ! The second scattering moment can be determined in a similar - ! manner to the first scattering moment - - score = last_wgt * calc_pn(2, mu) - - case (SCORE_SCATTER_3) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! The third scattering moment can be determined in a similar - ! manner to the first scattering moment - - score = last_wgt * calc_pn(3, mu) - - case (SCORE_SCATTER_4) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! The fourth scattering moment can be determined in a similar - ! manner to the first scattering moment - - score = last_wgt * calc_pn(4, mu) - - case (SCORE_SCATTER_5) - ! Skip any event where the particle didn't scatter - if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP - - ! The fifth scattering moment can be determined in a similar - ! manner to the first scattering moment - - score = last_wgt * calc_pn(5, mu) case (SCORE_TRANSPORT) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP @@ -233,7 +222,7 @@ contains ! mu*Sigma_s)). score = last_wgt / (3.0_8 * score * (material_xs % total - & - mu * score)) + mu * score)) case (SCORE_N_1N) ! Skip any event where the particle didn't scatter @@ -326,6 +315,7 @@ contains case default message = "Invalid score type on tally " // to_str(t % id) // "." + write(*,*) score_index, score_bin, l,j, t % n_score_bins call fatal_error() end select @@ -393,7 +383,7 @@ contains ! change outgoing energy bin t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & - size(t % filters(i) % real_bins), E_out) + size(t % filters(i) % real_bins), E_out) ! determine scoring index i_filter = sum((t % matching_bins - 1) * t % stride) + 1 @@ -474,7 +464,7 @@ contains if (t % all_nuclides) then call score_all_nuclides(tracklength_tallies(curr_ptr % data), flux, & - filter_index) + filter_index) else NUCLIDE_BIN_LOOP: do k = 1, t % n_nuclide_bins @@ -511,20 +501,20 @@ contains select case(score_bin) case (SCORE_TOTAL) score = micro_xs(i_nuclide) % total * & - atom_density * flux + atom_density * flux case (SCORE_SCATTER) score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * & - atom_density * flux + micro_xs(i_nuclide) % absorption) * & + atom_density * flux case (SCORE_ABSORPTION) score = micro_xs(i_nuclide) % absorption * & - atom_density * flux + atom_density * flux case (SCORE_FISSION) score = micro_xs(i_nuclide) % fission * & - atom_density * flux + atom_density * flux case (SCORE_NU_FISSION) score = micro_xs(i_nuclide) % nu_fission * & - atom_density * flux + atom_density * flux case (SCORE_EVENTS) score = ONE case default @@ -634,7 +624,7 @@ contains score = micro_xs(i_nuclide) % total * atom_density * flux case (SCORE_SCATTER) score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * atom_density * flux + micro_xs(i_nuclide) % absorption) * atom_density * flux case (SCORE_ABSORPTION) score = micro_xs(i_nuclide) % absorption * atom_density * flux case (SCORE_FISSION) @@ -777,27 +767,27 @@ contains ! determine next universe bin ! TODO: Account for multiple universes when performing this filter t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, & - p % coord % universe, i_tally) + p % coord % universe, i_tally) case (FILTER_MATERIAL) t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, & - p % material, i_tally) + p % material, i_tally) case (FILTER_CELL) ! determine next cell bin ! TODO: Account for cells in multiple levels when performing this filter t % matching_bins(i) = get_next_bin(FILTER_CELL, & - p % coord % cell, i_tally) + p % coord % cell, i_tally) case (FILTER_CELLBORN) ! determine next cellborn bin t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, & - p % cell_born, i_tally) + p % cell_born, i_tally) case (FILTER_SURFACE) ! determine next surface bin t % matching_bins(i) = get_next_bin(FILTER_SURFACE, & - p % surface, i_tally) + p % surface, i_tally) case (FILTER_ENERGYIN) ! determine incoming energy bin @@ -805,12 +795,12 @@ contains ! check if energy of the particle is within energy bins if (p % E < t % filters(i) % real_bins(1) .or. & - p % E > t % filters(i) % real_bins(k + 1)) then + p % E > t % filters(i) % real_bins(k + 1)) then t % matching_bins(i) = NO_BIN_FOUND else ! search to find incoming energy bin t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & - k + 1, p % E) + k + 1, p % E) end if end select @@ -931,25 +921,25 @@ contains select case(score_bin) case (SCORE_TOTAL) score = micro_xs(i_nuclide) % total * & - atom_density * flux + atom_density * flux case (SCORE_SCATTER) score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * & - atom_density * flux + micro_xs(i_nuclide) % absorption) * & + atom_density * flux case (SCORE_ABSORPTION) score = micro_xs(i_nuclide) % absorption * & - atom_density * flux + atom_density * flux case (SCORE_FISSION) score = micro_xs(i_nuclide) % fission * & - atom_density * flux + atom_density * flux case (SCORE_NU_FISSION) score = micro_xs(i_nuclide) % nu_fission * & - atom_density * flux + atom_density * flux case (SCORE_EVENTS) score = ONE case default message = "Invalid score type on tally " // & - to_str(t % id) // "." + to_str(t % id) // "." call fatal_error() end select @@ -972,7 +962,7 @@ contains score = ONE case default message = "Invalid score type on tally " // & - to_str(t % id) // "." + to_str(t % id) // "." call fatal_error() end select end if @@ -1030,27 +1020,27 @@ contains ! determine next universe bin ! TODO: Account for multiple universes when performing this filter t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, & - p % coord % universe, i_tally) + p % coord % universe, i_tally) case (FILTER_MATERIAL) t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, & - p % material, i_tally) + p % material, i_tally) case (FILTER_CELL) ! determine next cell bin ! TODO: Account for cells in multiple levels when performing this filter t % matching_bins(i) = get_next_bin(FILTER_CELL, & - p % coord % cell, i_tally) + p % coord % cell, i_tally) case (FILTER_CELLBORN) ! determine next cellborn bin t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, & - p % cell_born, i_tally) + p % cell_born, i_tally) case (FILTER_SURFACE) ! determine next surface bin t % matching_bins(i) = get_next_bin(FILTER_SURFACE, & - p % surface, i_tally) + p % surface, i_tally) case (FILTER_ENERGYIN) ! determine incoming energy bin @@ -1065,12 +1055,12 @@ contains ! check if energy of the particle is within energy bins if (E < t % filters(i) % real_bins(1) .or. & - E > t % filters(i) % real_bins(n + 1)) then + E > t % filters(i) % real_bins(n + 1)) then t % matching_bins(i) = NO_BIN_FOUND else ! search to find incoming energy bin t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & - n + 1, E) + n + 1, E) end if case (FILTER_ENERGYOUT) @@ -1079,12 +1069,12 @@ contains ! check if energy of the particle is within energy bins if (E < t % filters(i) % real_bins(1) .or. & - E > t % filters(i) % real_bins(n + 1)) then + E > t % filters(i) % real_bins(n + 1)) then t % matching_bins(i) = NO_BIN_FOUND else ! search to find incoming energy bin t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & - n + 1, p % E) + n + 1, p % E) end if end select @@ -1175,14 +1165,14 @@ contains n = t % filters(j) % n_bins ! check if energy of the particle is within energy bins if (p % E < t % filters(j) % real_bins(1) .or. & - p % E > t % filters(j) % real_bins(n + 1)) then + p % E > t % filters(j) % real_bins(n + 1)) then curr_ptr => curr_ptr % next ! select next tally cycle end if ! search to find incoming energy bin t % matching_bins(j) = binary_search(t % filters(j) % real_bins, & - n + 1, p % E) + n + 1, p % E) end if ! ======================================================================= @@ -1200,7 +1190,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_TOP t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 call add_to_score(t % scores(1, filter_index), p % wgt) end if @@ -1211,7 +1201,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_TOP t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 call add_to_score(t % scores(1, filter_index), p % wgt) end if @@ -1227,7 +1217,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_FRONT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 call add_to_score(t % scores(1, filter_index), p % wgt) end if @@ -1238,7 +1228,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_FRONT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 call add_to_score(t % scores(1, filter_index), p % wgt) end if @@ -1254,7 +1244,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_RIGHT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 call add_to_score(t % scores(1, filter_index), p % wgt) end if @@ -1265,7 +1255,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_RIGHT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 call add_to_score(t % scores(1, filter_index), p % wgt) end if @@ -1318,7 +1308,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_RIGHT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if ijk0(1) = ijk0(1) + 1 xyz_cross(1) = xyz_cross(1) + m % width(1) @@ -1330,7 +1320,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_RIGHT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if end if elseif (distance == d(2)) then @@ -1340,7 +1330,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_FRONT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if ijk0(2) = ijk0(2) + 1 xyz_cross(2) = xyz_cross(2) + m % width(2) @@ -1352,7 +1342,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_FRONT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if end if else if (distance == d(3)) then @@ -1362,7 +1352,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_TOP t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if ijk0(3) = ijk0(3) + 1 xyz_cross(3) = xyz_cross(3) + m % width(3) @@ -1374,7 +1364,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_TOP t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if end if end if @@ -1385,7 +1375,7 @@ contains ! Check for errors if (filter_index <= 0 .or. filter_index > & - t % total_filter_bins) then + t % total_filter_bins) then message = "Score index outside range." call fatal_error() end if @@ -1442,7 +1432,7 @@ contains end if i_tally_check = tally_maps(filter_type) % items(filter_value) % & - elements(position(filter_type)) % index_tally + elements(position(filter_type)) % index_tally if (i_tally_check > i_tally) then ! Since the index being checked against is greater than the index we @@ -1454,7 +1444,7 @@ contains elseif (i_tally_check == i_tally) then ! Found a match bin = tally_maps(filter_type) % items(filter_value) % & - elements(position(filter_type)) % index_bin + elements(position(filter_type)) % index_bin return end if @@ -1539,14 +1529,14 @@ contains ! The MPI_IN_PLACE specifier allows the master to copy values into ! a receive buffer without having a temporary variable call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, & - MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Transfer values to value on master t % scores(:,:) % value = tally_temp else ! Receive buffer not significant at other processors call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, & - MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Reset value on other processors t % scores(:,:) % value = 0 @@ -1562,14 +1552,14 @@ contains if (master) then call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, & - MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Transfer values back to global_tallies on master global_tallies(:) % value = global_temp else ! Receive buffer not significant at other processors call MPI_REDUCE(global_temp, dummy, N_GLOBAL_TALLIES, & - MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Reset value on other processors global_tallies(:) % value = ZERO @@ -1580,11 +1570,11 @@ contains ! last realization if (master) then call MPI_REDUCE(MPI_IN_PLACE, total_weight, 1, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + 0, MPI_COMM_WORLD, mpi_err) else ! Receive buffer not significant at other processors call MPI_REDUCE(total_weight, dummy, 1, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + 0, MPI_COMM_WORLD, mpi_err) end if if (associated(curr_ptr)) nullify(curr_ptr) @@ -1690,7 +1680,7 @@ contains score % sum = score % sum/n score % sum_sq = sqrt((score % sum_sq/n - score % sum * & - score % sum) / (n - 1)) + score % sum) / (n - 1)) end subroutine statistics_score @@ -1911,7 +1901,7 @@ contains end if do i = n_cmfd_current_tallies + n_user_current_tallies, & - n_user_current_tallies + 1, -1 + n_user_current_tallies + 1, -1 ! allocate node allocate(curr_ptr) diff --git a/src/tally_header.F90 b/src/tally_header.F90 index f47f0ab1d9..47d992f4dd 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -97,8 +97,13 @@ module tally_header logical :: all_nuclides = .false. ! Values to score, e.g. flux, absorption, etc. + ! scat_order is the scattering order for each score. + ! It is to be 0 if the scattering order is 0, or if the score is not a + ! scattering response. integer :: n_score_bins = 0 integer, allocatable :: score_bins(:) + integer, allocatable :: scatt_order(:) + integer :: n_nonPN_score_bins = 0 ! Scores for each bin -- the first dimenion of the array is for scores ! (e.g. flux, total reaction rate, fission reaction rate, etc.) and the