crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Changed more variables to non-intrinsic names.

Browse source
This commit is contained in:
Paul Romano 2011-11-21 11:25:07 -08:00
parent b78676e7c8
commit 4ffd4bd160
3 changed files with 220 additions and 220 deletions
Download patch file Download diff file Expand all files Collapse all files

130
src/cross_section.f90
View file

@ -210,10 +210,10 @@ contains
! subroutines to parse the actual data.
!===============================================================================
subroutine read_ace_table(index_table, index)
subroutine read_ace_table(index_table, index_list)
integer, intent(in) :: index_table ! index in nuclides/sab_tables
integer, intent(in) :: index ! index in xs_listings
integer, intent(in) :: index_list ! index in xs_listings
integer :: i ! loop index for XSS records
integer :: j, j1, j2 ! indices in XSS
@ -230,7 +230,7 @@ contains
logical :: file_exists ! does ACE library exist?
character(7) :: readable ! is ACE library readable?
character(10) :: name ! name of ACE table
character(10) :: date ! date ACE library was processed
character(10) :: date_ ! date ACE library was processed
character(10) :: mat ! material identifier
character(70) :: comment ! comment for ACE table
character(MAX_FILE_LEN) :: filename ! path to ACE cross section library
@ -239,7 +239,7 @@ contains
type(XsListing), pointer :: listing => null()
! determine path, record length, and location of table
listing => xs_listings(index)
listing => xs_listings(index_list)
filename = listing % path
record_length = listing % recl
location = listing % location
@ -273,7 +273,7 @@ contains
end do
! Read first line of header
read(UNIT=in, FMT='(A10,2E12.0,1X,A10)') name, awr, kT, date
read(UNIT=in, FMT='(A10,2E12.0,1X,A10)') name, awr, kT, date_
! Read more header and NXS and JXS
read(UNIT=in, FMT=100) comment, mat, &
@ -296,7 +296,7 @@ contains
ACCESS='direct', RECL=record_length)
! Read all header information
read(UNIT=in, REC=location) name, awr, kT, date, &
read(UNIT=in, REC=location) name, awr, kT, date_, &
comment, mat, (zaids(i), awrs(i), i=1,16), NXS, JXS
! determine table length
@ -414,7 +414,7 @@ contains
integer :: LED ! location of energy distribution locators
integer :: LDIS ! location of all energy distributions
integer :: LOCC ! location of energy distributions for given MT
integer :: loc ! locator
integer :: lc ! locator
integer :: length ! length of data to allocate
type(DistEnergy), pointer :: edist => null()
@ -562,10 +562,10 @@ contains
! determine length of all precursor constants/yields/interp data
length = 0
loc = JXS(25)
lc = JXS(25)
do i = 1, NPCR
NR = int(XSS(loc + length + 1))
NE = int(XSS(loc + length + 2 + 2*NR))
NR = int(XSS(lc + length + 1))
NE = int(XSS(lc + length + 2 + 2*NR))
length = length + 3 + 2*NR + 2*NE
end do
@ -573,7 +573,7 @@ contains
allocate(nuc % nu_d_precursor_data(length))
! read delayed neutron precursor data
XSS_index = loc
XSS_index = lc
nuc % nu_d_precursor_data = get_real(length)
else
@ -839,7 +839,7 @@ contains
integer :: NR ! number of interpolation regions
integer :: NE ! number of incoming energies
integer :: IDAT ! location of first energy distribution for given MT
integer :: loc ! locator
integer :: lc ! locator
integer :: length ! length of data to allocate
integer :: length_interp_data ! length of interpolation data
@ -885,16 +885,16 @@ contains
edist % p_valid % y = get_real(NE)
! Set index to beginning of IDAT array
loc = LDIS + IDAT - 2
lc = LDIS + IDAT - 2
! determine length of energy distribution
length = length_energy_dist(loc, LAW, loc_law, length_interp_data)
length = length_energy_dist(lc, LAW, loc_law, length_interp_data)
! allocate secondary energy distribution array
allocate(edist % data(length))
! read secondary energy distribution
XSS_index = loc + 1
XSS_index = lc + 1
edist % data = get_real(length)
! read next energy distribution if present
@ -910,9 +910,9 @@ contains
! distribution array based on the secondary energy law and location in XSS
!===============================================================================
function length_energy_dist(loc, law, LOCC, lid) result(length)
function length_energy_dist(lc, law, LOCC, lid) result(length)
integer, intent(in) :: loc ! location in XSS array
integer, intent(in) :: lc ! location in XSS array
integer, intent(in) :: law ! energy distribution law
integer, intent(in) :: LOCC ! location of energy distribution
integer, intent(in) :: lid ! length of interpolation data
@ -934,9 +934,9 @@ contains
select case (law)
case (1)
! Tabular equiprobable energy bins
NR = int(XSS(loc + 1))
NE = int(XSS(loc + 2 + 2*NR))
NP = int(XSS(loc + 3 + 2*NR + NE))
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
NP = int(XSS(lc + 3 + 2*NR + NE))
length = 3 + 2*NR + NE + 3*NP*NE
case (2)
@ -949,72 +949,72 @@ contains
case (4)
! Continuous tabular distribution
NR = int(XSS(loc + 1))
NE = int(XSS(loc + 2 + 2*NR))
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
! determine length
NP = int(XSS(loc + length + 2))
NP = int(XSS(lc + length + 2))
length = length + 2 + 3*NP
! adjust location for this block
j = loc + 2 + 2*NR + NE + i
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
case (5)
! General evaporation spectrum
NR = int(XSS(loc + 1))
NE = int(XSS(loc + 2 + 2*NR))
NP = int(XSS(loc + 3 + 2*NR + 2*NE))
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
NP = int(XSS(lc + 3 + 2*NR + 2*NE))
length = 3 + 2*NR + 2*NE + NP
case (7)
! Maxwell fission spectrum
NR = int(XSS(loc + 1))
NE = int(XSS(loc + 2 + 2*NR))
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = 3 + 2*NR + 2*NE
case (9)
! Evaporation spectrum
NR = int(XSS(loc + 1))
NE = int(XSS(loc + 2 + 2*NR))
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = 3 + 2*NR + 2*NE
case (11)
! Watt spectrum
NRa = int(XSS(loc + 1))
NEa = int(XSS(loc + 2 + 2*NRa))
NRb = int(XSS(loc + 3 + 2*(NRa+NEa)))
NEb = int(XSS(loc + 4 + 2*(NRa+NEa+NRb)))
NRa = int(XSS(lc + 1))
NEa = int(XSS(lc + 2 + 2*NRa))
NRb = int(XSS(lc + 3 + 2*(NRa+NEa)))
NEb = int(XSS(lc + 4 + 2*(NRa+NEa+NRb)))
length = 5 + 2*(NRa + NEa + NRb + NEb)
case (44)
! Kalbach-Mann correlated scattering
NR = int(XSS(loc + 1))
NE = int(XSS(loc + 2 + 2*NR))
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
NP = int(XSS(loc + length + 2))
NP = int(XSS(lc + length + 2))
length = length + 2 + 5*NP
! adjust location for this block
j = loc + 2 + 2*NR + NE + i
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
case (61)
! Correlated energy and angle distribution
NR = int(XSS(loc + 1))
NE = int(XSS(loc + 2 + 2*NR))
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
! outgoing energy distribution
NP = int(XSS(loc + length + 2))
NP = int(XSS(lc + length + 2))
! adjust locators for angular distribution
do j = 1, NP
k = loc + length + 2 + 3*NP + j
k = lc + length + 2 + 3*NP + j
if (XSS(k) /= 0) XSS(k) = XSS(k) - LOCC - lid
end do
@ -1022,12 +1022,12 @@ contains
do j = 1, NP
! outgoing angle distribution -- NMU here is actually
! referred to as NP in the MCNP documentation
NMU = int(XSS(loc + length + 2))
NMU = int(XSS(lc + length + 2))
length = length + 2 + 3*NMU
end do
! adjust locators for energy distribution
j = loc + 2 + 2*NR + NE + i
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
@ -1037,9 +1037,9 @@ contains
case (67)
! Laboratory energy-angle law
NR = int(XSS(loc + 1))
NE = int(XSS(loc + 2 + 2*NR))
NMU = int(XSS(loc + 4 + 2*NR + 2*NE))
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
NMU = int(XSS(lc + 4 + 2*NR + 2*NE))
length = 4 + 2*(NR + NE + NMU)
end select
@ -1055,7 +1055,7 @@ contains
type(Nuclide), pointer :: nuc
integer :: JXS23 ! location of URR data
integer :: loc ! locator
integer :: lc ! locator
integer :: N ! # of incident energies
integer :: M ! # of probabilities
integer :: i ! index over incoming energies
@ -1070,19 +1070,19 @@ contains
nuc % urr_present = .true.
allocate(nuc % urr_data)
allocate(nuc % urr_data % params(6))
loc = JXS23
lc = JXS23
else
nuc % urr_present = .false.
return
end if
! read parameters
nuc % urr_data % params(1) = int(XSS(loc)) ! # of incident energies
nuc % urr_data % params(2) = int(XSS(loc + 1)) ! # of probabilities
nuc % urr_data % params(3) = int(XSS(loc + 2)) ! interpolation parameter
nuc % urr_data % params(4) = int(XSS(loc + 3)) ! inelastic competition flag
nuc % urr_data % params(5) = int(XSS(loc + 4)) ! other absorption flag
nuc % urr_data % params(6) = int(XSS(loc + 5)) ! factors flag
nuc % urr_data % params(1) = int(XSS(lc)) ! # of incident energies
nuc % urr_data % params(2) = int(XSS(lc + 1)) ! # of probabilities
nuc % urr_data % params(3) = int(XSS(lc + 2)) ! interpolation parameter
nuc % urr_data % params(4) = int(XSS(lc + 3)) ! inelastic competition flag
nuc % urr_data % params(5) = int(XSS(lc + 4)) ! other absorption flag
nuc % urr_data % params(6) = int(XSS(lc + 5)) ! factors flag
! allocate incident energies and probability tables
N = nuc % urr_data % params(1)
@ -1091,7 +1091,7 @@ contains
allocate(nuc % urr_data % prob(N,6,M))
! read incident energies
XSS_index = loc + 6
XSS_index = lc + 6
nuc % urr_data % energy = get_real(N)
! read probability tables
@ -1120,7 +1120,7 @@ contains
integer :: i ! index for incoming energies
integer :: j ! index for outgoing energies
integer :: k ! index for outoging angles
integer :: loc ! location in XSS array
integer :: lc ! location in XSS array
integer :: NE_in ! number of incoming energies
integer :: NE_out ! number of outgoing energies
integer :: NMU ! number of outgoing angles
@ -1153,19 +1153,19 @@ contains
allocate(table % inelastic_mu(NMU, NE_out, NE_in))
! read outgoing energy/angle distribution for inelastic scattering
loc = JXS(3) - 1
lc = JXS(3) - 1
do i = 1, NE_in
do j = 1, NE_out
! read outgoing energy
table % inelastic_e_out(j,i) = XSS(loc + 1)
table % inelastic_e_out(j,i) = XSS(lc + 1)
! read outgoing angles for this outgoing energy
do k = 1, NMU
table % inelastic_mu(k,j,i) = XSS(loc + 1 + k)
table % inelastic_mu(k,j,i) = XSS(lc + 1 + k)
end do
! advance pointer
loc = loc + 1 + NMU
lc = lc + 1 + NMU
end do
end do
@ -1202,12 +1202,12 @@ contains
! read equiprobable outgoing cosines for elastic scattering each
! incoming energy
if (JXS4 /= 0 .and. NMU /= 0) then
loc = JXS(6) - 1
lc = JXS(6) - 1
do i = 1, NE_in
do j = 1, NMU
table % elastic_mu(j,i) = XSS(loc + j)
table % elastic_mu(j,i) = XSS(lc + j)
end do
loc = loc + NMU
lc = lc + NMU
end do
end if

58
src/initialize.f90
View file

@ -197,9 +197,9 @@ contains
subroutine prepare_universes()
integer :: i ! index in cells array
integer :: index ! index in universes array
integer :: count ! number of cells in a universe
integer :: i ! index in cells array
integer :: i_univ ! index in universes array
integer :: n_cells_in_univ ! number of cells in a universe
integer, allocatable :: index_cell_in_univ(:) ! the index in the univ%cells
! array for each universe
type(ListKeyValueII), pointer :: key_list => null()
@ -216,19 +216,19 @@ contains
key_list => dict_keys(universe_dict)
do while (associated(key_list))
! find index of universe in universes array
index = key_list%data%value
univ => universes(index)
i_univ = key_list%data%value
univ => universes(i_univ)
univ % id = key_list%data%key
! check for lowest level universe
if (univ % id == 0) BASE_UNIVERSE = index
if (univ % id == 0) BASE_UNIVERSE = i_univ
! find cell count for this universe
count = dict_get_key(cells_in_univ_dict, key_list%data%key)
n_cells_in_univ = dict_get_key(cells_in_univ_dict, key_list%data%key)
! allocate cell list for universe
allocate(univ % cells(count))
univ % n_cells = count
allocate(univ % cells(n_cells_in_univ))
univ % n_cells = n_cells_in_univ
! move to next universe
key_list => key_list%next
@ -244,13 +244,13 @@ contains
c => cells(i)
! get pointer to corresponding universe
index = dict_get_key(universe_dict, c % universe)
univ => universes(index)
i_univ = dict_get_key(universe_dict, c % universe)
univ => universes(i_univ)
! increment the index for the cells array within the Universe object and
! then store the index of the Cell object in that array
index_cell_in_univ(index) = index_cell_in_univ(index) + 1
univ % cells(index_cell_in_univ(index)) = i
index_cell_in_univ(i_univ) = index_cell_in_univ(i_univ) + 1
univ % cells(index_cell_in_univ(i_univ)) = i
end do
end subroutine prepare_universes
@ -268,7 +268,7 @@ contains
integer :: i ! index in cells array
integer :: j ! index over surface list
integer :: k
integer :: index ! index in surfaces/materials array
integer :: i_array ! index in surfaces/materials array
integer :: id ! user-specified id
type(Cell), pointer :: c => null()
type(Lattice), pointer :: l => null()
@ -283,8 +283,8 @@ contains
id = c % surfaces(j)
if (id < OP_DIFFERENCE) then
if (dict_has_key(surface_dict, abs(id))) then
index = dict_get_key(surface_dict, abs(id))
c % surfaces(j) = sign(index, id)
i_array = dict_get_key(surface_dict, abs(id))
c % surfaces(j) = sign(i_array, id)
else
message = "Could not find surface " // trim(int_to_str(abs(id))) // &
& " specified on cell " // trim(int_to_str(c % id))
@ -464,9 +464,9 @@ contains
integer, intent(in) :: parent ! cell containing universe
integer, intent(in) :: level ! level of universe
integer :: i ! index for cells in universe
integer :: x,y ! indices for lattice positions
integer :: index ! index in cells array
integer :: i ! index for cells in universe
integer :: x,y ! indices for lattice positions
integer :: i_cell ! index in cells array
integer :: universe_num
type(Cell), pointer :: c => null()
type(Universe), pointer :: subuniverse => null()
@ -478,15 +478,15 @@ contains
! loop over all cells in the universe
do i = 1, univ % n_cells
index = univ % cells(i)
c => cells(index)
i_cell = univ % cells(i)
c => cells(i_cell)
c%parent = parent
! if this cell is filled with another universe, recursively
! call this subroutine
if (c % type == CELL_FILL) then
subuniverse => universes(c % fill)
call build_universe(subuniverse, index, level + 1)
call build_universe(subuniverse, i_cell, level + 1)
end if
! if this cell is filled by a lattice, need to build the
@ -503,7 +503,7 @@ contains
call dict_add_key(dict, universe_num, 0)
subuniverse => universes(universe_num)
call build_universe(subuniverse, index, level + 1)
call build_universe(subuniverse, i_cell, level + 1)
end if
end do
end do
@ -519,7 +519,7 @@ contains
subroutine normalize_ao()
integer :: index ! index used for several purposes
integer :: index_list ! index in xs_listings array
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: n ! length of string
@ -562,8 +562,8 @@ contains
end if
if (dict_has_key(xs_listing_dict, key)) then
index = dict_get_key(xs_listing_dict, key)
mat % xs_listing(j) = index
index_list = dict_get_key(xs_listing_dict, key)
mat % xs_listing(j) = index_list
else
message = "Cannot find cross-section " // trim(key) // &
" in specified cross_sections.xml file."
@ -571,7 +571,7 @@ contains
end if
! determine atomic weight ratio
awr = xs_listings(index) % awr
awr = xs_listings(index_list) % awr
! if given weight percent, convert all values so that they are divided
! by awr. thus, when a sum is done over the values, it's actually
@ -592,8 +592,8 @@ contains
if (.not. density_in_atom) then
sum_percent = ZERO
do j = 1, mat % n_nuclides
index = mat % xs_listing(j)
awr = xs_listings(index) % awr
index_list = mat % xs_listing(j)
awr = xs_listings(index_list) % awr
x = mat % atom_percent(j)
sum_percent = sum_percent + x*awr
end do

252
src/physics.f90
View file

@ -1009,7 +1009,7 @@ contains
integer :: i ! loop index
integer :: j ! index on nu energy grid / precursor group
integer :: loc ! index before start of energies/nu values
integer :: lc ! index before start of energies/nu values
integer :: NR ! number of interpolation regions
integer :: NE ! number of energies tabulated
integer :: nu ! actual number of neutrons produced
@ -1091,23 +1091,23 @@ contains
! sampled delayed precursor group
xi = prn()
loc = 1
lc = 1
prob = ZERO
do j = 1, nuc % n_precursor
! determine number of interpolation regions and energies
NR = int(nuc % nu_d_precursor_data(loc + 1))
NE = int(nuc % nu_d_precursor_data(loc + 2 + 2*NR))
NR = int(nuc % nu_d_precursor_data(lc + 1))
NE = int(nuc % nu_d_precursor_data(lc + 2 + 2*NR))
! determine delayed neutron precursor yield for group j
yield = interpolate_tab1(nuc % nu_d_precursor_data( &
loc+1:loc+2+2*NR+2*NE), E)
lc+1:lc+2+2*NR+2*NE), E)
! Check if this group is sampled
prob = prob + yield
if (xi < prob) exit
! advance pointer
loc = loc + 2 + 2*NR + 2*NE + 1
lc = lc + 2 + 2*NR + 2*NE + 1
end do
! sample from energy distribution for group j
@ -1250,7 +1250,7 @@ contains
integer :: type ! angular distribution type
integer :: i ! incoming energy bin
integer :: n ! number of incoming energy bins
integer :: loc ! location in data array
integer :: lc ! location in data array
integer :: NP ! number of points in cos distribution
integer :: k ! index on cosine grid
real(8) :: r ! interpolation factor on incoming energy
@ -1290,7 +1290,7 @@ contains
if (r > prn()) i = i + 1
! check whether this is a 32-equiprobable bin or a tabular distribution
loc = rxn % adist % location(i)
lc = rxn % adist % location(i)
type = rxn % adist % type(i)
if (type == ANGLE_ISOTROPIC) then
mu = TWO * prn() - ONE
@ -1300,26 +1300,26 @@ contains
k = 1 + int(32.0_8*xi)
! calculate cosine
mu0 = rxn % adist % data(loc + k)
mu1 = rxn % adist % data(loc + k+1)
mu0 = rxn % adist % data(lc + k)
mu1 = rxn % adist % data(lc + k+1)
mu = mu0 + (32.0_8 * xi - k) * (mu1 - mu0)
elseif (type == ANGLE_TABULAR) then
interp = int(rxn % adist % data(loc + 1))
NP = int(rxn % adist % data(loc + 2))
interp = int(rxn % adist % data(lc + 1))
NP = int(rxn % adist % data(lc + 2))
! determine outgoing cosine bin
xi = prn()
loc = loc + 2
c_k = rxn % adist % data(loc + 2*NP + 1)
lc = lc + 2
c_k = rxn % adist % data(lc + 2*NP + 1)
do k = 1, NP-1
c_k1 = rxn % adist % data(loc + 2*NP + k+1)
c_k1 = rxn % adist % data(lc + 2*NP + k+1)
if (xi < c_k1) exit
c_k = c_k1
end do
p0 = rxn % adist % data(loc + NP + k)
mu0 = rxn % adist % data(loc + k)
p0 = rxn % adist % data(lc + NP + k)
mu0 = rxn % adist % data(lc + k)
if (interp == HISTOGRAM) then
! Histogram interpolation
mu = mu0 + (xi - c_k)/p0
@ -1327,8 +1327,8 @@ contains
elseif (interp == LINEAR_LINEAR) then
! Linear-linear interpolation -- not sure how you come about the
! formula given in the MCNP manual
p1 = rxn % adist % data(loc + NP + k+1)
mu1 = rxn % adist % data(loc + k+1)
p1 = rxn % adist % data(lc + NP + k+1)
mu1 = rxn % adist % data(lc + k+1)
frac = (p1 - p0)/(mu1 - mu0)
if (frac == ZERO) then
@ -1418,7 +1418,7 @@ contains
integer :: i ! index on incoming energy grid
integer :: k ! sampled index on outgoing grid
integer :: l ! sampled index on incoming grid
integer :: loc ! dummy index
integer :: lc ! dummy index
integer :: NR ! number of interpolation regions
integer :: NE ! number of energies
integer :: NET ! number of outgoing energies
@ -1500,23 +1500,23 @@ contains
end if
! determine index on incoming energy grid and interpolation factor
loc = 2 + 2*NR
i = binary_search(edist % data(loc+1), NE, E_in)
r = (E_in - edist%data(loc+i)) / &
& (edist%data(loc+i+1) - edist%data(loc+i))
lc = 2 + 2*NR
i = binary_search(edist % data(lc+1), NE, E_in)
r = (E_in - edist%data(lc+i)) / &
& (edist%data(lc+i+1) - edist%data(lc+i))
! Sample outgoing energy bin
r1 = prn()
k = 1 + int(NET * r1)
! Determine E_1 and E_K
loc = 3 + 3*NR + NE + (i-1)*NET
E_i_1 = edist % data(loc + 1)
E_i_K = edist % data(loc + NET)
lc = 3 + 3*NR + NE + (i-1)*NET
E_i_1 = edist % data(lc + 1)
E_i_K = edist % data(lc + NET)
loc = 3 + 3*NR + NE + i*NET
E_i1_1 = edist % data(loc + 1)
E_i1_K = edist % data(loc + NET)
lc = 3 + 3*NR + NE + i*NET
E_i1_1 = edist % data(lc + 1)
E_i1_K = edist % data(lc + NET)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
@ -1530,9 +1530,9 @@ contains
end if
! Determine E_l_k and E_l_k+1
loc = 3 + 2*NR + NE + (l-1)*NET
E_l_k = edist % data(loc+k)
E_l_k1 = edist % data(loc+k+1)
lc = 3 + 2*NR + NE + (l-1)*NET
E_l_k = edist % data(lc+k)
E_l_k1 = edist % data(lc+k+1)
! Determine E' (denoted here as E_out)
r2 = prn()
@ -1569,17 +1569,17 @@ contains
! find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last
! bins
loc = 2 + 2*NR
if (E_in < edist % data(loc+1)) then
lc = 2 + 2*NR
if (E_in < edist % data(lc+1)) then
i = 1
r = ZERO
elseif (E_in > edist % data(loc+NE)) then
elseif (E_in > edist % data(lc+NE)) then
i = NE - 1
r = ONE
else
i = binary_search(edist % data(loc+1), NE, E_in)
r = (E_in - edist%data(loc+i)) / &
& (edist%data(loc+i+1) - edist%data(loc+i))
i = binary_search(edist % data(lc+1), NE, E_in)
r = (E_in - edist%data(lc+i)) / &
& (edist%data(lc+i+1) - edist%data(lc+i))
end if
! Sample between the ith and (i+1)th bin
@ -1591,25 +1591,25 @@ contains
end if
! interpolation for energy E1 and EK
loc = int(edist%data(2 + 2*NR + NE + i))
NP = int(edist%data(loc + 2))
E_i_1 = edist%data(loc + 2 + 1)
E_i_K = edist%data(loc + 2 + NP)
lc = int(edist%data(2 + 2*NR + NE + i))
NP = int(edist%data(lc + 2))
E_i_1 = edist%data(lc + 2 + 1)
E_i_K = edist%data(lc + 2 + NP)
loc = int(edist%data(2 + 2*NR + NE + i + 1))
NP = int(edist%data(loc + 2))
E_i1_1 = edist%data(loc + 2 + 1)
E_i1_K = edist%data(loc + 2 + NP)
lc = int(edist%data(2 + 2*NR + NE + i + 1))
NP = int(edist%data(lc + 2))
E_i1_1 = edist%data(lc + 2 + 1)
E_i1_K = edist%data(lc + 2 + NP)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! determine location of outgoing energies, pdf, cdf for E(l)
loc = int(edist % data(2 + 2*NR + NE + l))
lc = int(edist % data(2 + 2*NR + NE + l))
! determine type of interpolation and number of discrete lines
INTTp = int(edist % data(loc + 1))
NP = int(edist % data(loc + 2))
INTTp = int(edist % data(lc + 1))
NP = int(edist % data(lc + 2))
if (INTTp > 10) then
INTT = mod(INTTp,10)
ND = (INTTp - INTT)/10
@ -1627,16 +1627,16 @@ contains
! determine outgoing energy bin
r1 = prn()
loc = loc + 2 ! start of EOUT
c_k = edist % data(loc + 2*NP + 1)
lc = lc + 2 ! start of EOUT
c_k = edist % data(lc + 2*NP + 1)
do k = 1, NP-1
c_k1 = edist % data(loc + 2*NP + k+1)
c_k1 = edist % data(lc + 2*NP + k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
E_l_k = edist % data(loc+k)
p_l_k = edist % data(loc+NP+k)
E_l_k = edist % data(lc+k)
p_l_k = edist % data(lc+NP+k)
if (INTT == HISTOGRAM) then
! Histogram interpolation
E_out = E_l_k + (r1 - c_k)/p_l_k
@ -1644,8 +1644,8 @@ contains
elseif (INTT == LINEAR_LINEAR) then
! Linear-linear interpolation -- not sure how you come about the
! formula given in the MCNP manual
E_l_k1 = edist % data(loc+k+1)
p_l_k1 = edist % data(loc+NP+k+1)
E_l_k1 = edist % data(lc+k+1)
p_l_k1 = edist % data(lc+NP+k+1)
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
if (frac == ZERO) then
@ -1694,8 +1694,8 @@ contains
T = interpolate_tab1(edist % data, E_in)
! determine restriction energy
loc = 2 + 2*NR + 2*NE
U = edist % data(loc + 1)
lc = 2 + 2*NR + 2*NE
U = edist % data(lc + 1)
! sample outgoing energy based on evaporation spectrum probability
! density function
@ -1719,8 +1719,8 @@ contains
Watt_a = interpolate_tab1(edist % data, E_in)
! determine Watt parameter 'b' from tabulated function
loc = 3 + 2*(NR + NE)
Watt_b = interpolate_tab1(edist % data, E_in, loc)
lc = 3 + 2*(NR + NE)
Watt_b = interpolate_tab1(edist % data, E_in, lc)
! Sample energy-dependent Watt fission spectrum
E_out = watt_spectrum(Watt_a, Watt_b)
@ -1748,17 +1748,17 @@ contains
! find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last
! bins
loc = 2 + 2*NR
if (E_in < edist % data(loc+1)) then
lc = 2 + 2*NR
if (E_in < edist % data(lc+1)) then
i = 1
r = ZERO
elseif (E_in > edist % data(loc+NE)) then
elseif (E_in > edist % data(lc+NE)) then
i = NE - 1
r = ONE
else
i = binary_search(edist % data(loc+1), NE, E_in)
r = (E_in - edist%data(loc+i)) / &
& (edist%data(loc+i+1) - edist%data(loc+i))
i = binary_search(edist % data(lc+1), NE, E_in)
r = (E_in - edist%data(lc+i)) / &
& (edist%data(lc+i+1) - edist%data(lc+i))
end if
! Sample between the ith and (i+1)th bin
@ -1770,26 +1770,26 @@ contains
end if
! determine endpoints on grid i
loc = int(edist%data(2+2*NR+NE + i)) ! start of LDAT for i
NP = int(edist%data(loc + 2))
E_i_1 = edist%data(loc + 2 + 1)
E_i_K = edist%data(loc + 2 + NP)
lc = int(edist%data(2+2*NR+NE + i)) ! start of LDAT for i
NP = int(edist%data(lc + 2))
E_i_1 = edist%data(lc + 2 + 1)
E_i_K = edist%data(lc + 2 + NP)
! determine endpoints on grid i+1
loc = int(edist%data(2+2*NR+NE + i+1)) ! start of LDAT for i+1
NP = int(edist%data(loc + 2))
E_i1_1 = edist%data(loc + 2 + 1)
E_i1_K = edist%data(loc + 2 + NP)
lc = int(edist%data(2+2*NR+NE + i+1)) ! start of LDAT for i+1
NP = int(edist%data(lc + 2))
E_i1_1 = edist%data(lc + 2 + 1)
E_i1_K = edist%data(lc + 2 + NP)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! determine location of outgoing energies, pdf, cdf for E(l)
loc = int(edist % data(2 + 2*NR + NE + l))
lc = int(edist % data(2 + 2*NR + NE + l))
! determine type of interpolation and number of discrete lines
INTTp = int(edist % data(loc + 1))
NP = int(edist % data(loc + 2))
INTTp = int(edist % data(lc + 1))
NP = int(edist % data(lc + 2))
if (INTTp > 10) then
INTT = mod(INTTp,10)
ND = (INTTp - INTT)/10
@ -1807,29 +1807,29 @@ contains
! determine outgoing energy bin
r1 = prn()
loc = loc + 2 ! start of EOUT
c_k = edist % data(loc + 2*NP + 1)
lc = lc + 2 ! start of EOUT
c_k = edist % data(lc + 2*NP + 1)
do k = 1, NP-1
c_k1 = edist % data(loc + 2*NP + k+1)
c_k1 = edist % data(lc + 2*NP + k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
E_l_k = edist % data(loc+k)
p_l_k = edist % data(loc+NP+k)
E_l_k = edist % data(lc+k)
p_l_k = edist % data(lc+NP+k)
if (INTT == HISTOGRAM) then
! Histogram interpolation
E_out = E_l_k + (r1 - c_k)/p_l_k
! Determine Kalbach-Mann parameters
KM_R = edist % data(loc + 3*NP + k)
KM_A = edist % data(loc + 4*NP + k)
KM_R = edist % data(lc + 3*NP + k)
KM_A = edist % data(lc + 4*NP + k)
elseif (INTT == LINEAR_LINEAR) then
! Linear-linear interpolation -- not sure how you come about the
! formula given in the MCNP manual
E_l_k1 = edist % data(loc+k+1)
p_l_k1 = edist % data(loc+NP+k+1)
E_l_k1 = edist % data(lc+k+1)
p_l_k1 = edist % data(lc+NP+k+1)
! Find E prime
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
@ -1841,10 +1841,10 @@ contains
end if
! Determine Kalbach-Mann parameters
R_k = edist % data(loc + 3*NP + k)
R_k1 = edist % data(loc + 3*NP + k+1)
A_k = edist % data(loc + 4*NP + k)
A_k1 = edist % data(loc + 4*NP + k+1)
R_k = edist % data(lc + 3*NP + k)
R_k1 = edist % data(lc + 3*NP + k+1)
A_k = edist % data(lc + 4*NP + k)
A_k1 = edist % data(lc + 4*NP + k+1)
KM_R = R_k + (R_k1 - R_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k)
KM_A = A_k + (A_k1 - A_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k)
@ -1891,17 +1891,17 @@ contains
! find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last
! bins
loc = 2 + 2*NR
if (E_in < edist % data(loc+1)) then
lc = 2 + 2*NR
if (E_in < edist % data(lc+1)) then
i = 1
r = ZERO
elseif (E_in > edist % data(loc+NE)) then
elseif (E_in > edist % data(lc+NE)) then
i = NE - 1
r = ONE
else
i = binary_search(edist % data(loc+1), NE, E_in)
r = (E_in - edist%data(loc+i)) / &
& (edist%data(loc+i+1) - edist%data(loc+i))
i = binary_search(edist % data(lc+1), NE, E_in)
r = (E_in - edist%data(lc+i)) / &
& (edist%data(lc+i+1) - edist%data(lc+i))
end if
! Sample between the ith and (i+1)th bin
@ -1913,26 +1913,26 @@ contains
end if
! determine endpoints on grid i
loc = int(edist%data(2+2*NR+NE + i)) ! start of LDAT for i
NP = int(edist%data(loc + 2))
E_i_1 = edist%data(loc + 2 + 1)
E_i_K = edist%data(loc + 2 + NP)
lc = int(edist%data(2+2*NR+NE + i)) ! start of LDAT for i
NP = int(edist%data(lc + 2))
E_i_1 = edist%data(lc + 2 + 1)
E_i_K = edist%data(lc + 2 + NP)
! determine endpoints on grid i+1
loc = int(edist%data(2+2*NR+NE + i+1)) ! start of LDAT for i+1
NP = int(edist%data(loc + 2))
E_i1_1 = edist%data(loc + 2 + 1)
E_i1_K = edist%data(loc + 2 + NP)
lc = int(edist%data(2+2*NR+NE + i+1)) ! start of LDAT for i+1
NP = int(edist%data(lc + 2))
E_i1_1 = edist%data(lc + 2 + 1)
E_i1_K = edist%data(lc + 2 + NP)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! determine location of outgoing energies, pdf, cdf for E(l)
loc = int(edist % data(2 + 2*NR + NE + l))
lc = int(edist % data(2 + 2*NR + NE + l))
! determine type of interpolation and number of discrete lines
INTTp = int(edist % data(loc + 1))
NP = int(edist % data(loc + 2))
INTTp = int(edist % data(lc + 1))
NP = int(edist % data(lc + 2))
if (INTTp > 10) then
INTT = mod(INTTp,10)
ND = (INTTp - INTT)/10
@ -1950,16 +1950,16 @@ contains
! determine outgoing energy bin
r1 = prn()
loc = loc + 2 ! start of EOUT
c_k = edist % data(loc + 2*NP + 1)
lc = lc + 2 ! start of EOUT
c_k = edist % data(lc + 2*NP + 1)
do k = 1, NP-1
c_k1 = edist % data(loc + 2*NP + k+1)
c_k1 = edist % data(lc + 2*NP + k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
E_l_k = edist % data(loc+k)
p_l_k = edist % data(loc+NP+k)
E_l_k = edist % data(lc+k)
p_l_k = edist % data(lc+NP+k)
if (INTT == HISTOGRAM) then
! Histogram interpolation
E_out = E_l_k + (r1 - c_k)/p_l_k
@ -1967,8 +1967,8 @@ contains
elseif (INTT == LINEAR_LINEAR) then
! Linear-linear interpolation -- not sure how you come about the
! formula given in the MCNP manual
E_l_k1 = edist % data(loc+k+1)
p_l_k1 = edist % data(loc+NP+k+1)
E_l_k1 = edist % data(lc+k+1)
p_l_k1 = edist % data(lc+NP+k+1)
! Find E prime
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
@ -1991,30 +1991,30 @@ contains
end if
! Find location of correlated angular distribution
loc = int(edist % data(loc+3*NP+k))
lc = int(edist % data(lc+3*NP+k))
! Check if angular distribution is isotropic
if (loc == 0) then
if (lc == 0) then
mu_out = TWO * prn() - ONE
return
end if
! interpolation type and number of points in angular distribution
JJ = int(edist % data(loc + 1))
NP = int(edist % data(loc + 2))
JJ = int(edist % data(lc + 1))
NP = int(edist % data(lc + 2))
! determine outgoing cosine bin
r3 = prn()
loc = loc + 2
c_k = edist % data(loc + 2*NP + 1)
lc = lc + 2
c_k = edist % data(lc + 2*NP + 1)
do k = 1, NP-1
c_k1 = edist % data(loc + 2*NP + k+1)
c_k1 = edist % data(lc + 2*NP + k+1)
if (r3 < c_k1) exit
c_k = c_k1
end do
p_k = edist % data(loc + NP + k)
mu_k = edist % data(loc + k)
p_k = edist % data(lc + NP + k)
mu_k = edist % data(lc + k)
if (JJ == HISTOGRAM) then
! Histogram interpolation
mu_out = mu_k + (r3 - c_k)/p_k
@ -2022,8 +2022,8 @@ contains
elseif (JJ == LINEAR_LINEAR) then
! Linear-linear interpolation -- not sure how you come about the
! formula given in the MCNP manual
p_k1 = edist % data(loc + NP + k+1)
mu_k1 = edist % data(loc + k+1)
p_k1 = edist % data(lc + NP + k+1)
mu_k1 = edist % data(lc + k+1)
frac = (p_k1 - p_k)/(mu_k1 - mu_k)
if (frac == ZERO) then