diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 1b7177db0d..7b18695de5 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -186,7 +186,6 @@ search.o: constants.o search.o: error.o search.o: global.o -source.o: ace_header.o source.o: bank_header.o source.o: constants.o source.o: error.o diff --git a/src/ace.F90 b/src/ace.F90 index 85f91a8495..13e4caff7b 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -937,7 +937,9 @@ contains integer, intent(in) :: lid ! length of interpolation data integer :: length ! length of energy distribution (LDAT) - integer :: i,j,k + integer :: i ! loop index for incoming energies + integer :: j ! loop index for outgoing energies + integer :: k ! dummy index in XSS integer :: NR ! number of interpolation regions integer :: NE ! number of incoming energies integer :: NP ! number of points in outgoing energy distribution diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 3caa5eb563..f3e305c391 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -61,8 +61,8 @@ module ace_header integer :: inelastic_flag ! inelastic competition flag integer :: absorption_flag ! other absorption flag logical :: multiply_smooth ! multiply by smooth cross section? - real(8), allocatable :: energy(:) - real(8), allocatable :: prob(:,:,:) + real(8), allocatable :: energy(:) ! incident energies + real(8), allocatable :: prob(:,:,:) ! actual probabibility tables end type UrrData !=============================================================================== diff --git a/src/constants.F90 b/src/constants.F90 index 43ca861901..9168727771 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -77,15 +77,15 @@ module constants ! Lattice types integer, parameter :: & - LATTICE_RECT = 1, & - LATTICE_HEX = 2 + LATTICE_RECT = 1, & ! Rectangular lattice + LATTICE_HEX = 2 ! Hexagonal lattice ! Lattice boundary crossings integer, parameter :: & - LATTICE_LEFT = 1, & - LATTICE_RIGHT = 2, & - LATTICE_BOTTOM = 3, & - LATTICE_TOP = 4 + LATTICE_LEFT = 1, & ! Flag for crossing left lattice boundary + LATTICE_RIGHT = 2, & ! Flag for crossing right lattice boundary + LATTICE_BOTTOM = 3, & ! Flag for crossing bottom lattice boundary + LATTICE_TOP = 4 ! Flag for crossing top lattice boundary ! Surface types integer, parameter :: & @@ -320,4 +320,3 @@ module constants integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file end module constants - diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 39f5fcd31f..3fb1d940d2 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -20,11 +20,11 @@ contains subroutine unionized_grid() - integer :: i ! index in materials array - integer :: j ! index over nuclides in material - integer :: index_list ! index in xs_listings array - character(12) :: name ! name of isotope, e.g. 92235.03c - character(12) :: alias ! alias of nuclide, e.g. U-235.03c + integer :: i ! index in materials array + integer :: j ! index over nuclides in material + integer :: index_list ! index in xs_listings array + character(12) :: name ! name of isotope, e.g. 92235.03c + character(12) :: alias ! alias of nuclide, e.g. U-235.03c type(ListReal), pointer :: list => null() type(ListReal), pointer :: current => null() type(Material), pointer :: mat => null() @@ -74,6 +74,9 @@ contains call list_delete(list) call dict_delete(already_added) + ! Set pointers to unionized energy grid for each nuclide + call grid_pointers() + end subroutine unionized_grid !=============================================================================== @@ -86,9 +89,9 @@ contains type(ListReal), pointer :: list real(8), intent(in) :: energy(:) - integer :: i - integer :: n - real(8) :: E + integer :: i ! index in energy array + integer :: n ! size of energy array + real(8) :: E ! actual energy value type(ListReal), pointer :: current => null() type(ListReal), pointer :: previous => null() type(ListReal), pointer :: head => null() @@ -169,22 +172,21 @@ contains end subroutine add_grid_points !=============================================================================== -! ORIGINAL_INDICES +! GRID_POINTERS creates an array of pointers (ints) for each nuclide to link +! each point on the nuclide energy grid to one on the unionized energy grid !=============================================================================== - subroutine original_indices() + subroutine grid_pointers() - integer :: i - integer :: j - integer :: index_e - integer :: n_grid_nuclide - real(8) :: union_energy - real(8) :: energy - type(Nuclide), pointer :: nuc + integer :: i ! loop index for nuclides + integer :: j ! loop index for nuclide energy grid + integer :: index_e ! index on union energy grid + real(8) :: union_energy ! energy on union grid + real(8) :: energy ! energy on nuclide grid + type(Nuclide), pointer :: nuc => null() do i = 1, n_nuclides_total nuc => nuclides(i) - n_grid_nuclide = size(nuc % energy) allocate(nuc % grid_index(n_grid)) index_e = 1 @@ -192,7 +194,7 @@ contains do j = 1, n_grid union_energy = e_grid(j) - if (union_energy >= energy .and. index_e < n_grid_nuclide) then + if (union_energy >= energy .and. index_e < nuc % n_grid) then index_e = index_e + 1 energy = nuc % energy(index_e) end if @@ -200,6 +202,6 @@ contains end do end do - end subroutine original_indices + end subroutine grid_pointers end module energy_grid diff --git a/src/initialize.F90 b/src/initialize.F90 index e8dda1d3d1..b3196ae4ad 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -6,7 +6,7 @@ module initialize use datatypes, only: dict_create, dict_add_key, dict_get_key, & dict_has_key, dict_keys use datatypes_header, only: ListKeyValueII, DictionaryII - use energy_grid, only: unionized_grid, original_indices + use energy_grid, only: unionized_grid use error, only: fatal_error use geometry, only: neighbor_lists use geometry_header, only: Cell, Surface, Universe, Lattice, BASE_UNIVERSE @@ -106,7 +106,6 @@ contains ! Construct unionized energy grid from cross-sections call timer_start(time_unionize) call unionized_grid() - call original_indices() call timer_stop(time_unionize) ! Create tally map @@ -143,11 +142,9 @@ contains subroutine setup_mpi() #ifdef MPI - integer :: i integer :: bank_blocks(4) ! Count for each datatype integer :: bank_types(4) ! Datatypes integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements - integer(MPI_ADDRESS_KIND) :: base type(Bank) :: b mpi_enabled = .true. @@ -187,10 +184,7 @@ contains call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err) ! Adjust displacements - base = bank_disp(1) - do i = 1, 4 - bank_disp(i) = bank_disp(i) - base - end do + bank_disp = bank_disp - bank_disp(1) ! Define MPI_BANK for fission sites bank_blocks = (/ 1, 3, 3, 1 /) @@ -365,11 +359,11 @@ contains subroutine adjust_indices() - integer :: i ! index in cells array - integer :: j ! index over surface list - integer :: k - integer :: i_array ! index in surfaces/materials array - integer :: id ! user-specified id + integer :: i ! index for various purposes + integer :: j ! index for various purposes + integer :: k ! loop index for lattices + integer :: i_array ! index in surfaces/materials array + integer :: id ! user-specified id type(Cell), pointer :: c => null() type(Lattice), pointer :: l => null() type(TallyObject), pointer :: t => null() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 49fac274ea..9d83c4b1df 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -485,14 +485,15 @@ contains use xml_data_materials_t - integer :: i, j - integer :: n - real(8) :: val - logical :: file_exists - character(3) :: default_xs - character(MAX_WORD_LEN) :: units - character(12) :: name - character(MAX_LINE_LEN) :: filename + integer :: i ! loop index for materials + integer :: j ! loop index for nuclides + integer :: n ! number of nuclides + real(8) :: val ! value entered for density + logical :: file_exists ! does materials.xml exist? + character(3) :: default_xs ! default xs identifier (e.g. 70c) + character(12) :: name ! name of nuclide + character(MAX_WORD_LEN) :: units ! units on density + character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml type(Material), pointer :: m => null() type(nuclide_xml), pointer :: nuc => null() type(sab_xml), pointer :: sab => null() @@ -1093,8 +1094,8 @@ contains use xml_data_plot_t - logical :: file_exists ! does tallies.xml file exist? - character(MAX_LINE_LEN) :: filename + logical :: file_exists ! does tallies.xml file exist? + character(MAX_LINE_LEN) :: filename ! absolute path to tallies.xml ! Check if plot.xml exists filename = trim(path_input) // "plot.xml" @@ -1150,7 +1151,7 @@ contains integer :: recl ! default record length integer :: entries ! default number of entries logical :: file_exists ! does cross_sections.xml exist? - character(MAX_WORD_LEN) :: directory + character(MAX_WORD_LEN) :: directory ! directory with cross sections type(XsListing), pointer :: listing => null() ! Check if cross_sections.xml exists diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 05f58313ff..3280d65a93 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -280,12 +280,12 @@ contains #ifdef MPI subroutine reduce_tallies() - integer :: i - integer :: n - integer :: m - integer :: n_bins - real(8), allocatable :: tally_temp(:,:) - type(TallyObject), pointer :: t + integer :: i ! loop index for tallies + integer :: n ! number of filter bins + integer :: m ! number of score bins + integer :: n_bins ! total number of bins + real(8), allocatable :: tally_temp(:,:) ! contiguous array of scores + type(TallyObject), pointer :: t => null() do i = 1, n_tallies t => tallies(i) diff --git a/src/mesh.F90 b/src/mesh.F90 index 69f7b47eaf..22c2f28c86 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -14,13 +14,13 @@ contains subroutine get_mesh_bin(m, xyz, bin) - type(StructuredMesh), pointer :: m - real(8), intent(in) :: xyz(:) - integer, intent(out) :: bin + type(StructuredMesh), pointer :: m ! mesh pointer + real(8), intent(in) :: xyz(:) ! coordinates + integer, intent(out) :: bin ! tally bin - integer :: n - integer :: ijk(3) - logical :: in_mesh + integer :: n ! size of mesh (2 or 3) + integer :: ijk(3) ! indices in mesh + logical :: in_mesh ! was given coordinate in mesh at all? ! Get number of dimensions n = m % n_dimension @@ -68,9 +68,9 @@ contains subroutine get_mesh_indices(m, xyz, ijk, in_mesh) type(StructuredMesh), pointer :: m - real(8), intent(in) :: xyz(:) - integer, intent(out) :: ijk(:) - logical, intent(out) :: in_mesh + real(8), intent(in) :: xyz(:) ! coordinates to check + integer, intent(out) :: ijk(:) ! indices in mesh + logical, intent(out) :: in_mesh ! were given coords in mesh? ! Find particle in mesh ijk = ceiling((xyz - m % lower_left)/m % width) @@ -95,11 +95,9 @@ contains integer, intent(in) :: ijk(:) integer :: bin - integer :: n_x - integer :: n_y - integer :: n_z + integer :: n_y ! number of mesh cells in y direction + integer :: n_z ! number of mesh cells in z direction - n_x = m % dimension(1) n_y = m % dimension(2) if (m % n_dimension == 2) then @@ -122,18 +120,16 @@ contains integer, intent(in) :: bin integer, intent(out) :: ijk(:) - integer :: n_y - integer :: n_z + integer :: n_y ! number of mesh cells in y direction + integer :: n_z ! number of mesh cells in z direction + + n_y = m % dimension(2) if (m % n_dimension == 2) then - n_y = m % dimension(2) - ijk(1) = (bin - 1)/n_y + 1 ijk(2) = mod(bin - 1, n_y) + 1 else if (m % n_dimension == 3) then - n_y = m % dimension(2) n_z = m % dimension(3) - ijk(1) = (bin - 1)/(n_y*n_z) + 1 ijk(2) = mod(bin - 1, n_y*n_z)/n_z + 1 ijk(3) = mod(bin - 1, n_z) + 1 diff --git a/src/output.F90 b/src/output.F90 index 72256f063a..e1226e1ea3 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -81,14 +81,14 @@ contains subroutine header(msg, unit, level) - character(*), intent(in) :: msg - integer, optional :: unit - integer, optional :: level + character(*), intent(in) :: msg ! header message + integer, optional :: unit ! unit to write to + integer, optional :: level ! specified header level integer :: n integer :: m - integer :: unit_ - integer :: header_level + integer :: unit_ ! unit to write to + integer :: header_level ! actual header level character(MAX_LINE_LEN) :: line ! set default level @@ -738,7 +738,7 @@ contains integer :: size_angle integer :: size_energy type(Reaction), pointer :: rxn => null() - type(UrrData), pointer :: urr => null() + type(UrrData), pointer :: urr => null() ! set default unit for writing information if (present(unit)) then diff --git a/src/physics.F90 b/src/physics.F90 index e203c0ca3c..4eb9852834 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1200,9 +1200,14 @@ contains real(8), intent(inout) :: w real(8), intent(in) :: mu ! cosine of angle in lab - real(8) :: phi, sinphi, cosphi - real(8) :: a,b - real(8) :: u0, v0, w0 + real(8) :: phi ! azimuthal angle + real(8) :: sinphi ! sine of azimuthal angle + real(8) :: cosphi ! cosine of azimuthal angle + real(8) :: a ! sqrt(1 - mu^2) + real(8) :: b ! sqrt(1 - w^2) + real(8) :: u0 ! original cosine in x direction + real(8) :: v0 ! original cosine in y direction + real(8) :: w0 ! original cosine in z direction ! Copy original directional cosines u0 = u @@ -1234,17 +1239,18 @@ contains end subroutine rotate_angle !=============================================================================== -! SAMPLE_ENERGY +! SAMPLE_ENERGY samples an outgoing energy distribution, either for a secondary +! neutron from a collision or for a prompt/delayed fission neutron !=============================================================================== recursive subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q) type(DistEnergy), pointer :: edist - real(8), intent(in) :: E_in - real(8), intent(out) :: E_out - real(8), intent(inout), optional :: mu_out - real(8), intent(in), optional :: A - real(8), intent(in), optional :: Q + real(8), intent(in) :: E_in ! incoming energy of neutron + real(8), intent(out) :: E_out ! outgoing energy + real(8), intent(inout), optional :: mu_out ! outgoing cosine of angle + real(8), intent(in), optional :: A ! mass number of nuclide + real(8), intent(in), optional :: Q ! Q-value of reaction integer :: i ! index on incoming energy grid integer :: k ! sampled index on outgoing grid @@ -1282,7 +1288,6 @@ contains real(8) :: p_k ! angular pdf in bin k real(8) :: p_k1 ! angular pdf in bin k+1 - real(8) :: E_cm real(8) :: r ! interpolation factor on incoming energy real(8) :: frac ! interpolation factor on outgoing energy real(8) :: U ! restriction energy @@ -1381,9 +1386,7 @@ contains ! ======================================================================= ! INELASTIC LEVEL SCATTERING - E_cm = edist%data(2) * (E_in - edist%data(1)) - - E_out = E_cm + E_out = edist%data(2) * (E_in - edist%data(1)) case (4) ! ======================================================================= diff --git a/src/source.F90 b/src/source.F90 index 8878179059..fd9cbf1c68 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -1,8 +1,7 @@ module source - use ace_header, only: Nuclide use bank_header, only: Bank - use constants, only: ONE, MAX_LINE_LEN + use constants, only: ONE use error, only: fatal_error use geometry_header, only: BASE_UNIVERSE use global @@ -163,12 +162,9 @@ contains subroutine initialize_particle() - ! TODO: if information on the cell, lattice, universe, and material is - ! passed through the fission bank to the source bank, no lookup would be - ! needed at the beginning of a cycle - - p % type = NEUTRON - p % alive = .true. + ! Set particle to neutron that's alive + p % type = NEUTRON + p % alive = .true. ! clear attributes p % surface = NONE diff --git a/src/string.F90 b/src/string.F90 index 2800d84c12..c40f828fca 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -159,7 +159,7 @@ contains integer :: i integer :: ic - do i = 1,len(word) + do i = 1, len(word) ic = ichar(word(i:i)) if (ic >= 65 .and. ic < 90) word(i:i) = char(ic+32) end do @@ -177,7 +177,7 @@ contains integer :: i integer :: ic - do i = 1,len(word) + do i = 1, len(word) ic = ichar(word(i:i)) if (ic >= 97 .and. ic < 122) word(i:i) = char(ic-32) end do @@ -214,7 +214,8 @@ contains character(*) :: str ! string to check character(*) :: seq ! sequence of characters - integer :: i, i_start + integer :: i + integer :: i_start integer :: str_len integer :: seq_len diff --git a/src/tally.F90 b/src/tally.F90 index fd108ff2bf..908e7879c7 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -239,16 +239,16 @@ contains subroutine score_analog_tally() - integer :: i - integer :: j - integer :: bins(N_FILTER_TYPES) - integer :: score_index - integer :: score_bin - real(8) :: score - real(8) :: last_wgt - real(8) :: wgt - real(8) :: mu - logical :: found_bin + integer :: i ! loop index for analog tallies + integer :: j ! loop index for scoring bins + integer :: bins(N_FILTER_TYPES) ! scoring bin combination + integer :: score_index ! single index for single bin + integer :: score_bin ! scoring bin, e.g. SCORE_FLUX + real(8) :: score ! analog tally score + real(8) :: last_wgt ! pre-collision particle weight + real(8) :: wgt ! post-collision particle weight + real(8) :: mu ! cosine of angle of collision + logical :: found_bin ! scoring bin found? type(TallyObject), pointer :: t => null() ! Copy particle's pre- and post-collision weight and angle @@ -465,13 +465,13 @@ contains type(TallyObject), pointer :: t integer, intent(inout) :: bins(N_FILTER_TYPES) - integer, intent(in) :: j + integer, intent(in) :: j ! index for score - integer :: k - integer :: bin_energyout - integer :: score_index - real(8) :: score - real(8) :: E_out + integer :: k ! loop index for bank sites + integer :: bin_energyout ! original outgoing energy bin + integer :: score_index ! index for scoring bin combination + real(8) :: score ! actualy score + real(8) :: E_out ! energy of fission bank site ! save original outgoing energy bin and score index bin_energyout = bins(FILTER_ENERGYOUT) @@ -515,14 +515,14 @@ contains real(8), intent(in) :: distance - integer :: i - integer :: j - integer :: bins(N_FILTER_TYPES) - integer :: score_index - integer :: score_bin - real(8) :: flux - real(8) :: score - logical :: found_bin + integer :: i ! loop index for tracklength tallies + integer :: j ! loop index for scoring bins + integer :: bins(N_FILTER_TYPES) ! scoring bin combination + integer :: score_index ! single index for single bin + integer :: score_bin ! scoring bin, e.g. SCORE_FLUX + real(8) :: flux ! tracklength estimate of flux + real(8) :: score ! actual score (e.g., flux*xs) + logical :: found_bin ! scoring bin found? type(TallyObject), pointer :: t => null() ! Determine track-length estimate of flux @@ -604,9 +604,9 @@ contains integer, intent(out) :: bins(N_FILTER_TYPES) logical, intent(out) :: found_bin - integer :: i - integer :: n - integer :: mesh_bin + integer :: i ! loop index for filters + integer :: n ! number of bins for single filter + integer :: mesh_bin ! index for mesh bin type(TallyObject), pointer :: t => null() type(StructuredMesh), pointer :: m => null() @@ -1424,7 +1424,7 @@ contains type(TallyObject), pointer :: t ! tally object integer, intent(in) :: filter_type ! type of filter integer, intent(in) :: bin ! bin in filter array - character(30) :: label ! user-specified identifier + character(30) :: label ! user-specified identifier integer :: i ! index in cells/surfaces/etc array integer, allocatable :: ijk(:) ! indices in mesh @@ -1486,7 +1486,7 @@ contains real(8) :: val2 ! sum(x*x) real(8) :: mean ! mean value real(8) :: std ! standard deviation of the mean - type(TallyObject), pointer :: t + type(TallyObject), pointer :: t => null() do i = 1, n_tallies t => tallies(i)