diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES
index 587fb06cca..bcc3b533ff 100644
--- a/src/DEPENDENCIES
+++ b/src/DEPENDENCIES
@@ -40,6 +40,7 @@ energy_grid.o: output.o
 finalize.o: global.o
 finalize.o: hdf5_interface.o
 finalize.o: output.o
+finalize.o: source.o
 finalize.o: tally.o
 finalize.o: timing.o
 
diff --git a/src/constants.F90 b/src/constants.F90
index 9cb5c99d55..4edebf1291 100644
--- a/src/constants.F90
+++ b/src/constants.F90
@@ -333,5 +333,6 @@ module constants
   integer, parameter :: UNIT_TALLY   = 12 ! unit # for writing tally file
   integer, parameter :: UNIT_PLOT    = 13 ! unit # for writing plot file
   integer, parameter :: UNIT_XS      = 14 ! unit # for writing xs summary file
+  integer, parameter :: UNIT_SOURCE  = 15 ! unit # for writing source file
 
 end module constants
diff --git a/src/finalize.F90 b/src/finalize.F90
index c4dcaf0623..e239c6db9f 100644
--- a/src/finalize.F90
+++ b/src/finalize.F90
@@ -2,6 +2,7 @@ module finalize
 
   use global
   use output,         only: print_runtime
+  use source,         only: write_source_binary
   use tally,          only: write_tallies, tally_statistics
   use timing,         only: timer_start, timer_stop
 
@@ -23,10 +24,13 @@ contains
     ! Start finalization timer
     call timer_start(time_finalize)
 
-    ! Calculate statistics for tallies and write to tallies.out
     if (.not. plotting) then
+       ! Calculate statistics for tallies and write to tallies.out
        call tally_statistics()
        if (master) call write_tallies()
+
+       ! Write out binary source
+       if (write_source) call write_source_binary()
     end if
 
     ! stop timers and show timing statistics
diff --git a/src/global.F90 b/src/global.F90
index 752e0db3a7..6cbb99cfa4 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -156,6 +156,9 @@ module global
   real(8), allocatable :: entropy_p(:,:,:) ! % of source sites in each cell
   type(StructuredMesh), pointer :: entropy_mesh
 
+  ! Write source at end of simulation
+  logical :: write_source = .false.
+
   ! ============================================================================
   ! PARALLEL PROCESSING VARIABLES
 
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index bb9841a14e..f90a54a247 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -68,6 +68,7 @@ contains
     verbosity_ = 0
     energy_grid_ = "union"
     seed_ = 0_8
+    write_source = ""
 
     ! Parse settings.xml file
     call read_xml_file_settings_t(filename)
@@ -231,6 +232,9 @@ contains
        entropy_on = .true.
     end if
 
+    ! Check if the user has specified to write binary source file
+    if (trim(write_source_) == 'on') write_source = .true.
+
   end subroutine read_settings_xml
 
 !===============================================================================
diff --git a/src/source.F90 b/src/source.F90
index e9088e3ddb..0652d22ac3 100644
--- a/src/source.F90
+++ b/src/source.F90
@@ -189,4 +189,31 @@ contains
 
   end subroutine initialize_particle
 
+!===============================================================================
+! WRITE_SOURCE writes out the final source distribution to a binary file that
+! can be used as a starting source in a new simulation
+!===============================================================================
+
+  subroutine write_source_binary()
+    
+    integer(8) :: i
+
+    open(UNIT=UNIT_SOURCE, FILE='source.binary', STATUS='replace', &
+         ACCESS='stream')
+
+    ! write the number of particles
+    write(UNIT=UNIT_SOURCE) n_particles
+
+    ! write information from the source bank
+    do i = 1_8, work
+       write(UNIT=UNIT_SOURCE) source_bank(i) % xyz
+       write(UNIT=UNIT_SOURCE) source_bank(i) % uvw
+       write(UNIT=UNIT_SOURCE) source_bank(i) % E
+    end do
+
+    ! close binary source file
+    close(UNIT=UNIT_SOURCE)
+
+  end subroutine write_source_binary
+
 end module source
diff --git a/src/utils/read_source.py b/src/utils/read_source.py
new file mode 100755
index 0000000000..6da8ff5008
--- /dev/null
+++ b/src/utils/read_source.py
@@ -0,0 +1,32 @@
+#!/usr/bin/env python
+
+import struct
+
+class SourceFile(object):
+    def __init__(self, filename):
+        # Open source file for reading
+        self.f = open(filename, 'r')
+
+        # Read number of source sites
+        self.n_sites = struct.unpack('l', self.f.read(8))[0]
+
+        # Create list to store source sites
+        self.sites = []
+
+        for i in range(self.n_sites):
+            # Read position, angle, and energy
+            xyz = struct.unpack('3d', self.f.read(24))
+            uvw = struct.unpack('3d', self.f.read(24))
+            E = struct.unpack('d', self.f.read(8))
+
+            # Create source site and append to list
+            self.sites.append(SourceSite(xyz, uvw, E))
+
+class SourceSite(object):
+    def __init__(self, xyz, uvw, E):
+        self.xyz = xyz
+        self.uvw = uvw
+        self.E = E
+
+    def __repr__(self):
+        return "<SourceSite: xyz={0} at E={1}>".format(self.xyz, self.E)
diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml
index dbe8bd6226..f2b387316d 100644
--- a/src/xml-fortran/templates/settings_t.xml
+++ b/src/xml-fortran/templates/settings_t.xml
@@ -37,5 +37,6 @@
   <variable name="survival_" tag="survival_biasing" type="word" length="3" />
   <variable name="verbosity_" tag="verbosity" type="integer" />
   <variable name="trace_" tag="trace" type="integer-array" />
+  <variable name="write_source_" tag="write_source" type="word" length="3" />
 
 </template>