diff --git a/docs/source/usersguide/mgxs_library.rst b/docs/source/usersguide/mgxs_library.rst
index d0dc1e54f..478b06002 100644
--- a/docs/source/usersguide/mgxs_library.rst
+++ b/docs/source/usersguide/mgxs_library.rst
@@ -114,20 +114,25 @@ attributes/sub-elements required to describe the meta-data:
     *Default*: "isotropic"
 
   :num_azimuthal:
-    This element provides the number of equi-width bins that the azimuthal
-    angular domain is subdivided in the case of angle-dependent cross sections
-    (i.e., "angle" is passed to the ``representation`` element).
+    This element provides the number of equal width angular bins that the
+    azimuthal angular domain is subdivided in the case of angle-dependent
+    cross sections (i.e., "angle" is passed to the ``representation`` element).
+    Note that these bins are equal in azimuthal angle widths, not equal in the
+    cosine of the azimuthal angle widths.
 
-    *Default*: If ``representation`` is "angle", this must be provided.  If
-    not, this parameter is not used.
+    *Default*: If ``representation`` is "angle", this must be provided.  This
+    parameter is not used for other ``representation`` types.
 
   :num_polar:
-    This element provides the number of equi-width bins that the polar angular
-    domain is subdivided in the case of angle-dependent cross sections
-    (i.e., "angle" is passed to the ``representation`` element).
+    This element provides the number of equal width angular bins that the
+    polar angular domain is subdivided in the case of angle-dependent
+    cross sections (i.e., "angle" is passed to the ``representation`` element).
+    Note that these bins are equal in polar angle widths, not equal in the
+    cosine of the polar angle widths.
 
-    *Default*: If ``representation`` is "angle", this must be provided.  If
-    not, this parameter is not used.
+
+    *Default*: If ``representation`` is "angle", this must be provided.  This
+    parameter is not used for other ``representation`` types.
 
   :scatt_type:
     This element provides the representation of the angular distribution
diff --git a/examples/python/pincell_multigroup/build-xml.py b/examples/python/pincell_multigroup/build-xml.py
index ba15370c3..17fbae994 100644
--- a/examples/python/pincell_multigroup/build-xml.py
+++ b/examples/python/pincell_multigroup/build-xml.py
@@ -72,7 +72,7 @@ mg_cross_sections_file.export_to_xml()
 uo2_data = openmc.Macroscopic('UO2', '300K')
 h2o_data = openmc.Macroscopic('LWTR', '300K')
 
-# Instantiate some Materials and register the appropriate Nuclides
+# Instantiate some Materials and register the appropriate Macroscopic objects
 uo2 = openmc.Material(material_id=1, name='UO2 fuel')
 uo2.set_density('macro', 1.0)
 uo2.add_macroscopic(uo2_data)
diff --git a/openmc/macroscopic.py b/openmc/macroscopic.py
index 094fa5042..9f67a50ba 100644
--- a/openmc/macroscopic.py
+++ b/openmc/macroscopic.py
@@ -8,7 +8,7 @@ if sys.version_info[0] >= 3:
 
 
 class Macroscopic(object):
-    """A nuclide that can be used in a material.
+    """A Macroscopic object that can be used in a material.
 
     Parameters
     ----------
diff --git a/src/global.F90 b/src/global.F90
index 9b68eab6a..dda4e0aea 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -490,10 +490,6 @@ contains
     end if
 
     if (allocated(macro_xs)) then
-    ! First call the clear routines
-      do i = 1, size(macro_xs)
-        call macro_xs(i) % obj % clear()
-      end do
       deallocate(macro_xs)
     end if
 
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index e726b6f40..2a21bac21 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -2071,7 +2071,7 @@ contains
         call get_list_item(node_macro_list, 1, node_nuc)
 
         ! Check for empty name on nuclide
-        if (.not.check_for_node(node_nuc, "name")) then
+        if (.not. check_for_node(node_nuc, "name")) then
           call fatal_error("No name specified on macroscopic data in material " &
                // trim(to_str(mat % id)))
         end if
diff --git a/src/macroxs.F90 b/src/macroxs.F90
index c26e9e888..f44f87f66 100644
--- a/src/macroxs.F90
+++ b/src/macroxs.F90
@@ -140,7 +140,7 @@ contains
 !===============================================================================
 
   subroutine macroxs_sample_scatter(scatt, gin, gout, mu, wgt)
-    Class(ScattData_Base), intent(in) :: scatt ! Scattering Object to Use
+    class(ScattData_Base), intent(in) :: scatt ! Scattering Object to Use
     integer,             intent(in)   :: gin   ! Incoming neutron group
     integer,            intent(out)   :: gout  ! Sampled outgoin group
     real(8),            intent(out)   :: mu    ! Sampled change in angle
diff --git a/src/macroxs_header.F90 b/src/macroxs_header.F90
index 63e6fecae..17ee8a1e0 100644
--- a/src/macroxs_header.F90
+++ b/src/macroxs_header.F90
@@ -19,10 +19,9 @@ module macroxs_header
     integer :: order
 
     ! Type-Bound procedures
-    contains
-      procedure(macroxs_init_),   deferred, pass :: init     ! initializes object
-      procedure(macroxs_clear_),  deferred, pass :: clear    ! Deallocates object
-      procedure(macroxs_get_xs_), deferred, pass :: get_xs   ! Return xs
+  contains
+    procedure(macroxs_init_),   deferred, pass :: init     ! initializes object
+    procedure(macroxs_get_xs_), deferred, pass :: get_xs   ! Return xs
   end type MacroXS_Base
 
   abstract interface
@@ -79,11 +78,10 @@ module macroxs_header
     real(8), allocatable :: scattxs(:)    ! scattering xs
     real(8), allocatable :: chi(:,:)      ! fission spectra
 
-    ! Type-Bound procedures
-    contains
-      procedure, pass :: init        => macroxs_iso_init        ! inits object
-      procedure, pass :: clear       => macroxs_iso_clear       ! Deallocates object
-      procedure, pass :: get_xs      => macroxs_iso_get_xs      ! Returns xs
+  ! Type-Bound procedures
+  contains
+    procedure, pass :: init        => macroxs_iso_init        ! inits object
+    procedure, pass :: get_xs      => macroxs_iso_get_xs      ! Returns xs
   end type MacroXS_Iso
 
   type, extends(MacroXS_Base) :: MacroXS_Angle
@@ -99,11 +97,10 @@ module macroxs_header
     real(8), allocatable :: polar(:)          ! polar angles
     real(8), allocatable :: azimuthal(:)      ! azimuthal angles
 
-    ! Type-Bound procedures
-    contains
-      procedure, pass :: init     => macroxs_angle_init   ! inits object
-      procedure, pass :: clear    => macroxs_angle_clear  ! Deallocates object
-      procedure, pass :: get_xs   => macroxs_angle_get_xs ! Returns xs
+  ! Type-Bound procedures
+  contains
+    procedure, pass :: init     => macroxs_angle_init   ! inits object
+    procedure, pass :: get_xs   => macroxs_angle_get_xs ! Returns xs
   end type MacroXS_Angle
 
 !===============================================================================
@@ -657,68 +654,6 @@ contains
 
   end subroutine macroxs_angle_init
 
-!===============================================================================
-! MACROXS*_CLEAR resets and deallocates data in MacroXS.
-!===============================================================================
-
-  subroutine macroxs_iso_clear(this)
-
-    class(MacroXS_Iso), intent(inout) :: this ! The MacroXS to clear
-
-    if (allocated(this % total)) then
-      deallocate(this % total, this % absorption, &
-                 this % nu_fission)
-    end if
-
-    if (allocated(this % fission)) then
-      deallocate(this % fission)
-    end if
-
-    if (allocated(this % k_fission)) then
-      deallocate(this % k_fission)
-    end if
-
-    call this % scatter % clear()
-
-    if (allocated(this % chi)) then
-      deallocate(this % chi)
-    end if
-
-  end subroutine macroxs_iso_clear
-
-  subroutine macroxs_angle_clear(this)
-
-    class(MacroXS_Angle), intent(inout) :: this ! The MacroXS to clear
-    integer :: i, j
-
-    if (allocated(this % total)) then
-      deallocate(this % total, this % absorption, &
-                 this % nu_fission)
-    end if
-
-    if (allocated(this % fission)) then
-      deallocate(this % fission)
-    end if
-
-    if (allocated(this % k_fission)) then
-      deallocate(this % k_fission)
-    end if
-
-    do i = 1, size(this % scatter,dim=2)
-      do j = 1, size(this % scatter,dim=1)
-        call this % scatter(j,i) % obj % clear()
-      end do
-    end do
-    if (allocated(this % scatter)) then
-      deallocate(this % scatter)
-    end if
-
-    if (allocated(this % chi)) then
-      deallocate(this % chi)
-    end if
-
-  end subroutine macroxs_angle_clear
-
 !===============================================================================
 ! MACROXS_*_GET_XS returns the requested data type
 !===============================================================================
@@ -790,130 +725,4 @@ contains
 
   end function macroxs_angle_get_xs
 
-!===============================================================================
-! THIN_GRID thins an (x,y) set while also thinning an associated y2
-!===============================================================================
-
-  subroutine thin_grid(xout, yout, yout2, tol, compression, maxerr)
-    real(8), allocatable, intent(inout) :: xout(:)     ! Resultant x grid
-    real(8), allocatable, intent(inout) :: yout(:)     ! Resultant y values
-    real(8), allocatable, intent(inout) :: yout2(:)     ! Secondary y values
-    real(8), intent(in)                 :: tol         ! Desired fractional error to maintain
-    real(8), intent(out)                :: compression ! Data reduction fraction
-    real(8), intent(inout)              :: maxerr      ! Maximum error due to compression
-
-    real(8), allocatable :: xin(:)      ! Incoming x grid
-    real(8), allocatable :: yin(:)      ! Incoming y values
-    real(8), allocatable :: yin2(:)     ! Secondary Incoming y values
-    integer :: k, klo, khi
-    integer :: all_ok
-    real(8) :: x1, y1, x2, y2, x, y, testval
-    integer :: num_keep, remove_it
-    real(8) :: initial_size
-    real(8) :: error
-    real(8) :: x_frac
-
-    initial_size = real(size(xout), 8)
-
-    allocate(xin(size(xout)))
-    xin = xout
-    allocate(yin(size(yout)))
-    yin = yout
-    allocate(yin2(size(yout2)))
-    yin2 = yout2
-
-    all_ok = size(yin)
-    maxerr = 0.0_8
-
-    ! This loop will step through each entry in dim==3 and check to see if
-    ! all of the values in other 2 dims can be replaced with linear interp.
-    ! If not, the value will be saved to a new array, if so, it will be
-    ! skipped.
-
-    xout = 0.0_8
-    yout = 0.0_8
-
-    ! Keep first point's data
-    xout(1) = xin(1)
-    yout(1) = yin(1)
-    yout2(1) = yin2(1)
-
-    ! Initialize data
-    num_keep = 1
-    klo = 1
-    khi = 3
-    k = 2
-    do while (khi <= size(xin))
-      remove_it = 0
-      x1 = xin(klo)
-      x2 = xin(khi)
-      x  = xin(k)
-      x_frac = 1.0_8 / (x2 - x1) * (x - x1)  ! Linear interp.
-
-      ! Check for removal. Otherwise, it stays. This is accomplished by leaving
-      ! remove_it as 0, entering the else portion of if(remove_it==all_ok)
-      y1 = yin(klo)
-      y2 = yin(khi)
-      y  = yin(k)
-
-      testval = y1 + (y2 - y1) * x_frac
-      error = abs(testval - y)
-      if (y /= 0.0_8) then
-        error = error / y
-      end if
-      if (error <= tol) then
-        remove_it = remove_it + 1
-        if (error > maxerr) then
-          maxerr = abs(testval - y)
-        end if
-      end if
-      ! Now place the point in to the proper bin and advance iterators.
-      if (remove_it /= 0) then
-        ! Then don't put it in the new grid but advance iterators
-        k = k + 1
-        khi = khi + 1
-      else
-        ! Put it in new grid and advance iterators accordingly
-        num_keep = num_keep + 1
-        xout(num_keep) = xin(k)
-        yout(num_keep) = yin(k)
-        yout2(num_keep) = yin2(k)
-        klo = k
-        k = k + 1
-        khi = khi + 1
-      end if
-    end do
-    ! Save the last point's data
-    num_keep = num_keep + 1
-    xout(num_keep) = xin(size(xin))
-    yout(num_keep) = yin(size(xin))
-    yout2(num_keep) = yin2(size(xin))
-
-    ! Finally, xout and yout were sized to match xin and yin since we knew
-    ! they would be no larger than those.  Now we must resize these arrays
-    ! and copy only the useful data in. Will use xin/yin for temp arrays.
-    xin = xout(1:num_keep)
-    yin = yout(1:num_keep)
-    yin2 = yout2(1:num_keep)
-
-    deallocate(xout)
-    deallocate(yout)
-    deallocate(yout2)
-    allocate(xout(num_keep))
-    allocate(yout(size(yin)))
-    allocate(yout2(size(yin2)))
-
-    xout = xin(1:num_keep)
-    yout = yin(1:num_keep)
-    yout2 = yin2(1:num_keep)
-
-    ! Clean up
-    deallocate(xin)
-    deallocate(yin)
-    deallocate(yin2)
-
-    compression = (initial_size - real(size(xout),8)) / initial_size
-
-  end subroutine thin_grid
-
 end module macroxs_header
diff --git a/src/math.F90 b/src/math.F90
index aedf18235..36c65a240 100644
--- a/src/math.F90
+++ b/src/math.F90
@@ -558,7 +558,7 @@ contains
   end function calc_rn
 
 !===============================================================================
-! EXPAND_HARMONIC expands a given series of harmonics
+! EXPAND_HARMONIC expands a given series of real spherical harmonics
 !===============================================================================
   pure function expand_harmonic(data, order, uvw) result(val)
     real(8), intent(in) :: data(:)
diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90
index 880764eda..f73f3f5a2 100644
--- a/src/mgxs_data.F90
+++ b/src/mgxs_data.F90
@@ -76,8 +76,7 @@ contains
           get_fiss = .true.
         end if
       end do
-      if (get_kfiss .and. get_fiss) &
-           exit
+      if (get_kfiss .and. get_fiss) exit
     end do
 
     ! ==========================================================================
diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90
index 398e91eb2..95f9833b4 100644
--- a/src/nuclide_header.F90
+++ b/src/nuclide_header.F90
@@ -155,7 +155,7 @@ module nuclide_header
     end function nuclide_calc_f_
   end interface
 
-  !===============================================================================
+!===============================================================================
 ! NUCLIDE_ISO contains the base MGXS data for a nuclide specifically for
 ! isotropically weighted MGXS
 !===============================================================================
@@ -688,11 +688,12 @@ module nuclide_header
     end function nuclide_angle_get_xs
 
 !===============================================================================
-! NUCLIDE_*_CALC_F Finds the value of f(mu), the scattering probability, given mu
+! NUCLIDE_*_CALC_F Finds the value of f(mu), the scattering angle probability,
+! given mu
 !===============================================================================
 
-  pure function nuclide_mg_iso_calc_f(this, gin, gout, mu, uvw, i_azi, i_pol) &
-       result(f)
+    pure function nuclide_mg_iso_calc_f(this, gin, gout, mu, uvw, i_azi, i_pol) &
+         result(f)
       class(Nuclide_Iso), intent(in) :: this
       integer, intent(in)            :: gin  ! Incoming Energy Group
       integer, intent(in)            :: gout ! Outgoing Energy Group
@@ -791,6 +792,7 @@ module nuclide_header
       end if
 
     end function nuclide_mg_angle_calc_f
+
 !===============================================================================
 ! find_angle finds the closest angle on the data grid and returns that index
 !===============================================================================
diff --git a/src/particle_header.F90 b/src/particle_header.F90
index 0e9fc0263..0426acd92 100644
--- a/src/particle_header.F90
+++ b/src/particle_header.F90
@@ -132,7 +132,6 @@ contains
     this % fission           = .false.
     this % delayed_group     = 0
     this % n_delayed_bank(:) = 0
-    ! Initialize this % g so there is always at least some initialized value
     this % g = 1
 
     ! Set up base level coordinates
diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90
index 29ab616e0..091a77741 100644
--- a/src/scattdata_header.F90
+++ b/src/scattdata_header.F90
@@ -16,11 +16,10 @@ module scattdata_header
     real(8), allocatable :: mult(:,:)   ! (Gout x Gin)
     real(8), allocatable :: data(:,:,:) ! (Order/Nmu x Gout x Gin)
 
-    ! Type-Bound procedures
-    contains
-      procedure(init_), deferred, pass   :: init   ! Initializes ScattData
-      procedure(calc_f_), deferred, pass :: calc_f ! Calculates f, given mu
-      procedure(clear_), deferred, pass  :: clear  ! Deallocates ScattData
+  ! Type-Bound procedures
+  contains
+    procedure(init_), deferred, pass   :: init   ! Initializes ScattData
+    procedure(calc_f_), deferred, pass :: calc_f ! Calculates f, given mu
   end type ScattData_Base
 
   abstract interface
@@ -42,37 +41,29 @@ module scattdata_header
       real(8)                           :: f    ! Return value of f(mu)
 
     end function calc_f_
-
-    subroutine clear_(this)
-      import ScattData_Base
-      class(ScattData_Base), intent(inout) :: this ! The ScattData to clear
-    end subroutine clear_
   end interface
 
   type, extends(ScattData_Base) :: ScattData_Legendre
-    contains
-      procedure, pass :: init   => scattdata_legendre_init
-      procedure, pass :: calc_f => scattdata_legendre_calc_f
-      procedure, pass :: clear  => scattdata_legendre_clear
+  contains
+    procedure, pass :: init   => scattdata_legendre_init
+    procedure, pass :: calc_f => scattdata_legendre_calc_f
   end type ScattData_Legendre
 
   type, extends(ScattData_Base) :: ScattData_Histogram
     real(8), allocatable :: mu(:) ! Mu bins
     real(8)              :: dmu   ! Mu spacing
-    contains
-      procedure, pass :: init   => scattdata_histogram_init
-      procedure, pass :: calc_f => scattdata_histogram_calc_f
-      procedure, pass :: clear  => scattdata_histogram_clear
+  contains
+    procedure, pass :: init   => scattdata_histogram_init
+    procedure, pass :: calc_f => scattdata_histogram_calc_f
   end type ScattData_Histogram
 
   type, extends(ScattData_Base) :: ScattData_Tabular
     real(8), allocatable :: mu(:)      ! Mu bins
     real(8)              :: dmu        ! Mu spacing
     real(8), allocatable :: fmu(:,:,:) ! PDF of f(mu)
-    contains
-      procedure, pass :: init   => scattdata_tabular_init
-      procedure, pass :: calc_f => scattdata_tabular_calc_f
-      procedure, pass :: clear  => scattdata_tabular_clear
+  contains
+    procedure, pass :: init   => scattdata_tabular_init
+    procedure, pass :: calc_f => scattdata_tabular_calc_f
   end type ScattData_Tabular
 
 !===============================================================================
@@ -239,56 +230,6 @@ contains
 
     end subroutine scattdata_tabular_init
 
-!===============================================================================
-! SCATTDATA_CLEAR resets and deallocates data in ScattData.
-!===============================================================================
-
-    subroutine scattdata_base_clear(this)
-      class(ScattData_Base), intent(inout) :: this
-
-      if (allocated(this % energy)) then
-        deallocate(this % energy)
-      end if
-
-      if (allocated(this % mult)) then
-        deallocate(this % mult)
-      end if
-
-      if (allocated(this % data)) then
-        deallocate(this % data)
-      end if
-
-    end subroutine scattdata_base_clear
-
-    subroutine scattdata_legendre_clear(this)
-      class(ScattData_Legendre), intent(inout) :: this
-
-      call scattdata_base_clear(this)
-
-    end subroutine scattdata_legendre_clear
-
-    subroutine scattdata_histogram_clear(this)
-      class(ScattData_Histogram), intent(inout) :: this
-
-      call scattdata_base_clear(this)
-
-      if (allocated(this % mu)) then
-        deallocate(this % mu)
-      end if
-
-    end subroutine scattdata_histogram_clear
-
-    subroutine scattdata_tabular_clear(this)
-      class(ScattData_Tabular), intent(inout) :: this
-
-      call scattdata_base_clear(this)
-
-      if (allocated(this % mu)) then
-        deallocate(this % mu)
-      end if
-
-    end subroutine scattdata_tabular_clear
-
 !===============================================================================
 ! SCATTDATA_*_CALC_F Calculates the value of f given mu (and gin,gout pair)
 !===============================================================================
@@ -350,6 +291,4 @@ contains
 
     end function scattdata_tabular_calc_f
 
-
-
 end module scattdata_header
\ No newline at end of file
diff --git a/src/tally.F90 b/src/tally.F90
index bbe099c3e..e127cf037 100644
--- a/src/tally.F90
+++ b/src/tally.F90
@@ -26,6 +26,8 @@ module tally
 
   integer :: position(N_FILTER_TYPES - 3) = 0 ! Tally map positioning array
 
+!$omp threadprivate(position)
+
   procedure(score_general_intfc),    pointer :: score_general => null()
   procedure(get_scoring_bins_intfc), pointer :: get_scoring_bins => null()
 
@@ -52,8 +54,6 @@ module tally
 
   end interface
 
-!$omp threadprivate(position)
-
 contains
 
 !===============================================================================
@@ -1255,97 +1255,95 @@ contains
     real(8) :: uvw(3)
 
     select case(score_bin)
-
-
-      case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
-        ! Find the scattering order for a singly requested moment, and
-        ! store its moment contribution.
-        if (t % moment_order(i) == 1) then
-          score = score * p % mu ! avoid function call overhead
-        else
-          score = score * calc_pn(t % moment_order(i), p % mu)
-        endif
+    case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
+      ! Find the scattering order for a singly requested moment, and
+      ! store its moment contribution.
+      if (t % moment_order(i) == 1) then
+        score = score * p % mu ! avoid function call overhead
+      else
+        score = score * calc_pn(t % moment_order(i), p % mu)
+      endif
 !$omp atomic
-        t % results(score_index, filter_index) % value = &
-             t % results(score_index, filter_index) % value + score
+      t % results(score_index, filter_index) % value = &
+           t % results(score_index, filter_index) % value + score
 
 
-      case(SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN)
-        score_index = score_index - 1
-        num_nm = 1
-        ! Find the order for a collection of requested moments
-        ! and store the moment contribution of each
-        do n = 0, t % moment_order(i)
-          ! determine scoring bin index
-          score_index = score_index + num_nm
-          ! Update number of total n,m bins for this n (m = [-n: n])
-          num_nm = 2 * n + 1
+    case(SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN)
+      score_index = score_index - 1
+      num_nm = 1
+      ! Find the order for a collection of requested moments
+      ! and store the moment contribution of each
+      do n = 0, t % moment_order(i)
+        ! determine scoring bin index
+        score_index = score_index + num_nm
+        ! Update number of total n,m bins for this n (m = [-n: n])
+        num_nm = 2 * n + 1
 
-          ! multiply score by the angular flux moments and store
+        ! multiply score by the angular flux moments and store
 !$omp critical (score_general_scatt_yn)
-          t % results(score_index: score_index + num_nm - 1, filter_index) &
-               % value = t &
-               % results(score_index: score_index + num_nm - 1, filter_index)&
-               % value &
-               + score * calc_pn(n, p % mu) * calc_rn(n, p % last_uvw)
+        t % results(score_index: score_index + num_nm - 1, filter_index) &
+             % value = t &
+             % results(score_index: score_index + num_nm - 1, filter_index)&
+             % value &
+             + score * calc_pn(n, p % mu) * calc_rn(n, p % last_uvw)
 !$omp end critical (score_general_scatt_yn)
-        end do
-        i = i + (t % moment_order(i) + 1)**2 - 1
+      end do
+      i = i + (t % moment_order(i) + 1)**2 - 1
 
 
-      case(SCORE_FLUX_YN, SCORE_TOTAL_YN)
-        score_index = score_index - 1
-        num_nm = 1
-        if (t % estimator == ESTIMATOR_ANALOG .or. &
-             t % estimator == ESTIMATOR_COLLISION) then
-          uvw = p % last_uvw
-        else if (t % estimator == ESTIMATOR_TRACKLENGTH) then
-          uvw = p % coord(1) % uvw
-        end if
-        ! Find the order for a collection of requested moments
-        ! and store the moment contribution of each
-        do n = 0, t % moment_order(i)
-          ! determine scoring bin index
-          score_index = score_index + num_nm
-          ! Update number of total n,m bins for this n (m = [-n: n])
-          num_nm = 2 * n + 1
+    case(SCORE_FLUX_YN, SCORE_TOTAL_YN)
+      score_index = score_index - 1
+      num_nm = 1
+      if (t % estimator == ESTIMATOR_ANALOG .or. &
+           t % estimator == ESTIMATOR_COLLISION) then
+        uvw = p % last_uvw
+      else if (t % estimator == ESTIMATOR_TRACKLENGTH) then
+        uvw = p % coord(1) % uvw
+      end if
+      ! Find the order for a collection of requested moments
+      ! and store the moment contribution of each
+      do n = 0, t % moment_order(i)
+        ! determine scoring bin index
+        score_index = score_index + num_nm
+        ! Update number of total n,m bins for this n (m = [-n: n])
+        num_nm = 2 * n + 1
 
-          ! multiply score by the angular flux moments and store
+        ! multiply score by the angular flux moments and store
 !$omp critical (score_general_flux_tot_yn)
-          t % results(score_index: score_index + num_nm - 1, filter_index) &
-               % value = t &
-               % results(score_index: score_index + num_nm - 1, filter_index)&
-               % value &
-               + score * calc_rn(n, uvw)
+        t % results(score_index: score_index + num_nm - 1, filter_index) &
+             % value = t &
+             % results(score_index: score_index + num_nm - 1, filter_index)&
+             % value &
+             + score * calc_rn(n, uvw)
 !$omp end critical (score_general_flux_tot_yn)
-        end do
-        i = i + (t % moment_order(i) + 1)**2 - 1
+      end do
+      i = i + (t % moment_order(i) + 1)**2 - 1
 
 
-      case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
-        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 % moment_order(i)
-          ! determine scoring bin index
-          score_index = score_index + 1
+    case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
+      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 % moment_order(i)
+        ! determine scoring bin index
+        score_index = score_index + 1
 
-          ! get the score and tally it
-!$omp atomic
-          t % results(score_index, filter_index) % value = &
-               t % results(score_index, filter_index) % value &
-               + score * calc_pn(n, p % mu)
-        end do
-        i = i + t % moment_order(i)
-
-
-      case default
+        ! get the score and tally it
 !$omp atomic
         t % results(score_index, filter_index) % value = &
-             t % results(score_index, filter_index) % value + score
+             t % results(score_index, filter_index) % value &
+             + score * calc_pn(n, p % mu)
+      end do
+      i = i + t % moment_order(i)
 
 
-      end select
+    case default
+!$omp atomic
+      t % results(score_index, filter_index) % value = &
+           t % results(score_index, filter_index) % value + score
+
+
+    end select
 
   end subroutine expand_and_score