diff --git a/src/cross_section.f90 b/src/cross_section.f90
index 99c3564226..271e412d95 100644
--- a/src/cross_section.f90
+++ b/src/cross_section.f90
@@ -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
     
diff --git a/src/initialize.f90 b/src/initialize.f90
index bd181da2ec..0a0250f874 100644
--- a/src/initialize.f90
+++ b/src/initialize.f90
@@ -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
diff --git a/src/physics.f90 b/src/physics.f90
index 65bcc3a0a6..7e9fed018e 100644
--- a/src/physics.f90
+++ b/src/physics.f90
@@ -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