Fixed several errors in tallies. Added synchronize_tallies subroutine.

src/cross_section.f90

@@ -598,6 +598,9 @@ contains
             .or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) then
           nuc % fissionable = .true.
           nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma
+
+          ! Also need to add fission cross sections to absorption
+          nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma
        end if
 
        ! set defaults

src/main.f90

@@ -11,6 +11,7 @@ program main
   use tally,           only: calculate_keff
   use source,          only: get_source_particle
   use string,          only: int_to_str
+  use tally,           only: synchronize_tallies
   use timing,          only: timer_start, timer_stop
 
 #ifdef MPI
@@ -95,6 +96,9 @@ contains
        ! Start timer for inter-cycle synchronization
        call timer_start(time_intercycle)
 
+       ! Collect tallies
+       call synchronize_tallies()
+
        ! Distribute fission bank across processors evenly
        call synchronize_bank(i_cycle)
 

src/physics.f90

@@ -202,6 +202,7 @@ contains
     micro_xs(i) % interp_factor = f
 
     ! Initialize nuclide cross-sections to zero
+    micro_xs(i) % fission    = ZERO
     micro_xs(i) % nu_fission = ZERO
 
     ! Calculate microscopic nuclide total cross section
@@ -287,12 +288,16 @@ contains
     type(Nuclide),  pointer :: nuc
     type(Reaction), pointer :: rxn
 
+    ! Add to collision counter for particle
+    p % n_collision = p % n_collision + 1
+
+    ! Get pointer to current material
     mat => materials(p % material)
 
     ! Score collision estimator tallies for any macro tallies -- we can do this
     ! before sampling the nuclide since 
     if (tallies_on) then
-       ! call score_tally(p)
+       call score_tally(p)
     end if
 
     ! sample nuclide

src/tally.f90

@@ -15,6 +15,8 @@ module tally
 
   implicit none
 
+  integer, allocatable :: position(:)
+
 contains
 
 !===============================================================================
@@ -109,6 +111,10 @@ contains
     allocate(tally_maps(MAP_MESH)     % items(100)) ! TODO: Change this
     allocate(tally_maps(MAP_BORNIN)   % items(n_cells))
 
+    ! Allocate and initialize tally map positioning for finding bins
+    allocate(position(TALLY_MAP_TYPES))
+    position = 0
+
     do i = 1, n_tallies
        t => tallies(i)
 
@@ -189,6 +195,7 @@ contains
        end if
 
        ! Allocate scores for tally
+       t % n_total_bins = filter_bins
        allocate(t % scores(filter_bins, score_bins))
 
     end do
@@ -348,7 +355,8 @@ contains
           case (MACRO_TOTAL)
              score = p % wgt
           case (MACRO_SCATTER)
-             score = p % wgt * material_xs % scatter / material_xs % total
+             score = p % wgt * (material_xs % total - material_xs % absorption) &
+                  / material_xs % total
           case (MACRO_ABSORPTION)
              score = p % wgt * material_xs % absorption / material_xs % total
           case (MACRO_FISSION)
@@ -364,6 +372,9 @@ contains
 
     end do
 
+    ! Reset tally map positioning
+    position = 0
+
   end subroutine score_tally
 
 !===============================================================================
@@ -379,16 +390,7 @@ contains
 
     integer :: index_tally
     integer :: index_bin
-    integer, allocatable, save :: position(:)
-    logical, save :: first_entry = .true.
-    integer       :: n
-
-    ! initialize position to zero
-    if (first_entry) then
-       allocate(position(TALLY_MAP_TYPES))
-       position = 0
-       first_entry = .false.
-    end if
+    integer :: n
 
     n = size(tally_maps(i_map) % items(i_cell) % elements)
 
@@ -500,4 +502,41 @@ contains
     
   end subroutine add_to_score
 
+!===============================================================================
+! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history
+! within the cycle to a new random variable
+!===============================================================================
+
+  subroutine synchronize_tallies()
+
+    integer :: i
+    integer :: j
+    integer :: k
+    real(8) :: val
+    type(TallyObject), pointer :: t
+
+
+    do i = 1, n_tallies
+       t => tallies(i)
+
+       do j = 1, t % n_total_bins
+          do k = 1, t % n_macro_bins
+             ! Add the sum and square of the sum of contributions from each
+             ! history within a cycle to the variables val and val_sq. This will
+             ! later allow us to calculate a variance on the tallies
+
+             val = t % scores(j,k) % val_history / n_particles
+             t % scores(j,k) % val    = t % scores(j,k) % val    + val
+             t % scores(j,k) % val_sq = t % scores(j,k) % val_sq + val*val
+
+             ! Reset the within-cycle accumulation variable
+
+             t % scores(j,k) % val_history = ZERO
+          end do
+       end do
+
+    end do
+
+  end subroutine synchronize_tallies
+
 end module tally

src/tally_header.f90

@@ -86,6 +86,7 @@ module tally_header
      integer :: n_bornin_bins   = 0
      integer :: n_energy_in     = 0
      integer :: n_energy_out    = 0
+     integer :: n_total_bins    = 0
 
      ! Macroscopic properties to score