anchor.diff:      ! Sort the sampled models for more efficient reading
       call heapsort(samp,nresamp)
   
  -    ! Reset file to read again
  +    ! Reset file to read again (do not overwrite nit set above)
       rewind(lu_log_file)
  -    call read_header(nit,nflt,npoles,log_format)
  +    call read_header(idum,nflt,npoles,log_format)
   
   	! Read search results, saving only sampled models
  -    call read_resampled_models(nit,nflt,npoles,log_format,locked,obj,slip,poles,temp,uncert,isAccepted)
  +    call read_resampled_models(nit,nflt,npoles,log_format,locked,obj,slip,poles,temp,uncert, &
  +                               isAccepted)
  +    nit = nresamp
   
   else
   
  
positive.code: integer, allocatable :: locked(:,:), plocked(:)
  double precision, allocatable :: temp(:)
  double precision, allocatable :: obj(:), obj_resamp(:)
  double precision, allocatable :: uncert(:)
  double precision :: min_obj, fraction_locked
  double precision, allocatable :: slip(:,:,:), mean_ss(:), mean_ds(:)
  double precision :: px, py, pz, plon, plat, pmag
  double precision, allocatable :: poles(:,:,:), mean_px(:), mean_py(:), mean_pz(:)
  double precision :: ddum
  logical :: ldum
  
  ios = 0
  
negative.code: return
  end subroutine
  
  
  !--------------------------------------------------------------------------------------------------!
  
  
  subroutine read_header(nit,nflt,npoles,log_format)
  
  