anchor.diff:          end do
   
       end subroutine reduce_slots
  -     
  -    subroutine insertion_rank_int(arr, rank, tolerance)
  -        
  -        ! ranking algorithm modified for integer valued arrays
  -
  -        ! Rank a real(p) array in increasing order using the insertion sort
  -        ! algorithm.
  -        !
  -        ! Resultant ranking is *stable* and insertion sort is really pretty
  -        ! decent, especially when the array is small.  Naturally one should
  -        ! investigate quicksort and the like for large arrays, but insertion
  -        ! sort is a good compromise and easy to code.
  -        !
  -        ! Based upon the F90 insertion sort implementation in Rosetta Stone
  -        ! http://rosettacode.org/wiki/Sorting_algorithms/Insertion_sort#Fortran.
  -        !
  -        ! ***WARNING***
  -        ! We assume that the arrays are 1-indexed due to a feature of how array
  -        ! bounds are handled in procedure arguments in Fortran.
  -        !
  +      
  +    subroutine find_processors(nparticles_proc, up_thresh, low_thresh, rec_dummy, don_dummy, donor_slots, proc_map)
  +        
  +        ! Find donor/receiver processors.
  +        ! Put these into varying size array receivers/donors.
  +
           ! In:
  -        !   arr: array of real values.
  -        !   tolerance (optional, default 0.0): tolerance for comparing values.
  -        !   If present, then two values which differ by less than tolerance are
  -        !   treated as equal.  This allows one to avoid changing the order of
  -        !   items that have identical values to within a desired precision.
  -        !   Naturally, one should have this as small as possible.
  -        ! In/Out:
  -        !    rank: on output rank(i) contains the ranked index (in increasing
  -        !    order) of the value in arr(i), that is arr(rank(i)) returns the
  -        !    i-th element in the sorted list of values of arr.
  -        !    NOTE: rank must have at least the dimensions of arr on input.
  -
  -        integer(lint), intent(in) :: arr(:)
  -        integer, intent(inout) :: rank(:) ! inout to avoid automatic deallocation
  -                                          ! of an allocatable array on entry
  -        integer :: i, j, tmp, tolerance,tol
  -
  -        forall (i=1:size(arr)) rank(i) = i
  -        tol=tolerance
  -        do i = 2, size(arr)
  -            j = i - 1
  -            tmp = rank(i)
  -            do while ( j >= 1 )
  -                if (arr(rank(j)) - arr(tmp) < tol) exit
  -                rank(j+1) = rank(j)
  -                j = j - 1
  -            end do
  -            rank(j+1) = tmp
  -        end do
  -
  -    end subroutine insertion_rank_int
  -    
  -    subroutine find_processors(rec_dummy, don_dummy, nparticles_proc, donor_slots)
  -        
  -        ! Find donor/receiver processors
  -        ! put these into varying size array receivers/donors
  -        ! In:
  -        !   nparticles_proc: num particles on each processor
  -        ! In/Out:
  -        !   rec_dummy/don_dummy: arrays which contain donor/receivers, entries are then kept in 
  -        !       smaller arrays 
  +        !   nparticles_proc: number particles on each processor.
  +        !   upper/lower_thresh: upper/lower thresholds for load imblance.
  +        ! In/Out:
  +        !   rec_dummy/don_dummy: arrays which contain donor/receivers processors.
           !   donor_slots: number of slots which we can donate, this varies as more entries in proc_map are 
           !       modified.
  +        !   proc_map: array which maps determinants to processors.
  +        !       proc_map(modulo(hash(d),num_slots*nprocs)=processor.
   
           use parallel, only: nprocs
  -        
  +        use ranking, only: insertion_rank_int
  +
  +        integer(lint), intent(in) :: nparticles_proc(0:nprocs-1)
  +        integer, intent(in) :: proc_map(0:p_map_size-1)
  +        integer, intent(in) :: up_thresh, low_thresh
  +        integer, intent(inout) :: donor_slots
  +
           integer ::  i, j, k, upper, lower
  -        integer, allocatable, dimension(:) ::  tmp_rec, tmp_don
  +        integer, allocatable, dimension(:) ::  tmp_rec, tmp_don, rec_sort
           integer, allocatable :: rec_dummy(:), don_dummy(:)
  -        integer(lint), intent(in) :: nparticles_proc(nprocs)
  -        integer, intent(inout) :: donor_slots
  +        integer :: rank_nparticles(nprocs)
  +
           allocate(tmp_rec(nprocs))
           allocate(tmp_don(nprocs))
  -
  +        k=1
           j=1
  -        k=1
  -        ! find donors/receivers
  -        do i=1, nprocs
  +
  +        ! Find donors/receivers processors.
  +        
  +        do i=0, size(nparticles_proc)-1
               if(nparticles_proc(i) .lt. low_thresh) then
  -                ! rank = 0 is a problem
  -                tmp_rec(j)=i-1
  +                tmp_rec(j)=i
                   j=j+1
               else if (nparticles_proc(i) .gt. up_thresh) then
  -                tmp_don(k)=i-1
  +                tmp_don(k)=i
                   k=k+1
               end if 
           end do
  -
  -        ! put index of processors into smaller array
  +        
  +        ! Put processor ID into smaller array.
           
           allocate(rec_dummy(j-1))
  +        allocate(rec_sort(j-1))
           allocate(don_dummy(k-1))
   
  -        do i=1, j-1
  -            rec_dummy(i)=tmp_rec(i)
  -        end do
  -        
  -        do i=1, k-1
  -            don_dummy(i)=tmp_don(i)
  -        end do
  -
  -        ! calculate number of donor slots which we can move
  +        don_dummy=tmp_don(:k-1)
  +        rec_dummy=tmp_rec(:j-1)
  +
  +        ! Sort receiver processers.
  +        
  +        call insertion_rank_int(nparticles_proc, rank_nparticles, 0) 
  +        do i=1, size(rec_dummy)
  +            rec_sort(i)=rank_nparticles(i)
  +        end do
  +        rec_dummy=rec_sort
  +
  +        ! Calculate number of donor slots which we can move.
           
           donor_slots=0
  -        do i=1, p_map_size
  +        do i=0, size(proc_map)-1
               do j=1, size(don_dummy)
                   if(proc_map(i)==don_dummy(j)) then
                       donor_slots=donor_slots+1
  
positive.code: module load_balancing
  
  ! Module for performing load balancing
  ! Determinants are assigned to processors using proc_map(num_slots*nprocs) which 
  ! is initialised so that its entries cyclically contain processors 0,..,nprocs-1
  ! Population for each entry in proc_map is calculated 
  ! Attempt to donate slots from processors with above average populations to those
  ! with below average populations
  
      use loadbal_data
      implicit none
      ! Upper and lower population thresholds
      integer :: up_thresh, low_thresh
      ! How much of a load imbalance is acceptable
      real(dp) :: perc_diff=0.01
      
  contains 
      
      subroutine redistribute_processors(proc_map)
          
          ! Main subroutine in module, carries out load balancing is relatively logical order
          ! In/Out:
          ! proc_map: array which maps determinants to processors
          !       proc_map(modulo(hash(d),num_slots*nprocs)=processor
  
negative.code: 
      end subroutine find_processors
      
      subroutine particles_per_proc(proc_map, slot_list, nparticles_proc)
  
          ! Find number of particles per processor and store in array nparticles_proc.
  
          ! In:
          !   proc_map: array which maps determinants to processors.
          !       proc_map(modulo(hash(d),num_slots*nprocs)=processor.
          !   slot_list: array containing populations of slots in proc_map across all
          !       processors.
          ! Out:
          !   nparticles_proc(nprocs): array containing population on each processor.
  
          use parallel, only : nprocs
  
          integer, intent(in) :: proc_map(0:p_map_size-1)
          integer(lint), intent(in) :: slot_list(0:p_map_size-1) 
          integer(lint), intent(out) :: nparticles_proc(0:nprocs-1)
          
          integer :: i
  
          nparticles_proc=0  
  
          do i=0, size(proc_map)-1
              nparticles_proc(proc_map(i))=nparticles_proc(proc_map(i))+slot_list(i)
          end do
  
      end subroutine particles_per_proc
  
      subroutine initialise_slot_pop(proc_map, num_slots, slot_pop)
  
          ! In: 
          !   proc_map(p_map_size): array which maps determinants to processors.
          !       proc_map(modulo(hash(d),num_slots*nprocs)=processor
          !   num_slots: number of slots which we divide slot_pop (and similar arrays) into.
          ! In/Out:
          !   slot_pop(p_map_size): array containing population of slots in proc_map
          !       processor dependendent.
  
          use parallel, only: nprocs, iproc
          use hashing 
          use basis, only: basis_length
          use determinants, only: det_info, alloc_det_info, dealloc_det_info
          use fciqmc_data, only: tot_walkers, walker_dets, walker_population
  
          integer, intent(in):: num_slots
          integer, intent(in) :: proc_map(0:p_map_size-1)
          integer(lint), intent(out) :: slot_pop(0:size(proc_map)-1)        
   
          integer :: i, det_pos
          type(det_info) :: cdet
  
          slot_pop=0
          do i=1, tot_walkers
              cdet%f => walker_dets(:,i)
              det_pos=modulo(murmurhash_bit_string(cdet%f, basis_length, 7),num_slots*nprocs)
              slot_pop(det_pos)=slot_pop(det_pos)+abs(walker_population(1,i))
          end do
     
  