anchor.diff:  module tally
  +
  +  use ISO_FORTRAN_ENV
   
     use constants
     use error,         only: fatal_error
  
positive.code: contains
  
  !===============================================================================
  ! CALCULATE_KEFF
  !===============================================================================
  
    subroutine calculate_keff(i_cycle)
  
      integer, intent(in) :: i_cycle ! index of current cycle
  
      integer(8)              :: total_bank ! total number of source sites
      integer                 :: n          ! active cycle number
      real(8)                 :: kcoll      ! keff collision estimator         
      real(8), save           :: k1 = 0.    ! accumulated keff
      real(8), save           :: k2 = 0.    ! accumulated keff**2
      real(8)                 :: std        ! stdev of keff over active cycles
  
      message = "Calculate cycle keff..."
      call write_message(8)
  
      ! set k1 and k2 at beginning of run
      if (i_cycle == 1) then
         k1 = ZERO
         k2 = ZERO
      end if
  
  #ifdef MPI
  
negative.code: 
      ! Collect statistics and print output
      if (master) then
         ! Since the creation of bank sites was originally weighted by the last
         ! cycle keff, we need to multiply by that keff to get the current cycle's
         ! value
  
         k_cycle = real(total_bank)/real(n_particles)*keff
  
         if (i_cycle > n_inactive) then
            ! Active cycle number
            n = i_cycle - n_inactive
  
            ! Accumulate cycle estimate of k
            k_sum =    k_sum    + k_cycle
            k_sum_sq = k_sum_sq + k_cycle*k_cycle
  
            ! Determine mean and standard deviation of mean
            keff = k_sum/n
            keff_std = sqrt((k_sum_sq/n - keff*keff)/n)
  
            ! Display output for this cycle
            if (i_cycle > n_inactive+1) then
               write(UNIT=OUTPUT_UNIT, FMT=101) i_cycle, k_cycle, keff, keff_std
            else
               write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle
            end if
         else
            ! Display output for inactive cycle
            write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle
            keff = k_cycle
         end if
      end if
  
  #ifdef MPI
      ! Broadcast new keff value to all processors
      call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err)
  #endif
  
  100 format (2X,I5,2X,F8.5)
  101 format (2X,I5,2X,F8.5,9X,F8.5,1X,F8.5)
  
    end subroutine calculate_keff
  
  