diff --git a/src/elements.f90 b/src/elements.f90
index e9fc928..0dffecc 100644
--- a/src/elements.f90
+++ b/src/elements.f90
@@ -5,6 +5,7 @@ module elements
     use parameters
     use functions
     use subroutines
+    use sorts
     use box
     implicit none
 
@@ -395,54 +396,69 @@ module elements
         return
     end subroutine rhombshape
 
-    subroutine delete_atoms(num, index)
+    subroutine delete_atoms(num, in_index)
         !This subroutine deletes atoms from the atom arrays 
         integer, intent(in) :: num
-        integer, intent(inout), dimension(num) :: index
+        integer, intent(in), dimension(num) :: in_index
+        
+        real(kind=dp), dimension(num) :: for_sort
+        integer, dimension(num) :: sorted_index
 
         integer :: i
 
-        call heapsort(index)
+        for_sort = in_index 
+
+        call dpquicksort(for_sort)
+        do i = 1, num
+            sorted_index(i) = nint(for_sort(i))
+        end do
 
         !We go from largest index to smallest index just to make sure that we don't miss 
         !accidentally overwrite values which need to be deleted
         do i = num, 1, -1
-            if(index(i) == atom_num) then 
-                r_atom(:,index(i)) = 0.0_dp
-                type_atom(index(i)) = 0
+            if(sorted_index(i) == atom_num) then 
+                r_atom(:,sorted_index(i)) = 0.0_dp
+                type_atom(sorted_index(i)) = 0
             else
-                r_atom(:,index(i)) = r_atom(:, atom_num)
-                type_atom(index(i)) = type_atom(atom_num)
+                r_atom(:,sorted_index(i)) = r_atom(:, atom_num)
+                type_atom(sorted_index(i)) = type_atom(atom_num)
             end if            
             atom_num = atom_num - 1
         end do 
     end subroutine
 
-    subroutine delete_elements(num, index)
+    subroutine delete_elements(num, in_index)
         !This subroutine deletes elements from the element array 
         integer, intent(in) :: num
-        integer, intent(inout), dimension(num) :: index
+        integer, intent(in), dimension(num) :: in_index
 
         integer :: i
+        real(kind=dp), dimension(num) :: for_sort
+        integer, dimension(num) :: sorted_index
 
-        call heapsort(index)
+        
+        for_sort = in_index
+        call dpquicksort(for_sort)
+        do i = 1, num
+            sorted_index(i) = nint(for_sort(i))
+        end do
 
         !We go from largest index to smallest index just to make sure that we don't miss 
         !accidentally overwrite values which need to be deleted
         do i = num, 1, -1
-            if(index(i) == ele_num) then 
-                r_node(:,:,:,index(i)) = 0.0_dp
-                type_ele(index(i)) =''
-                size_ele(index(i)) = 0
-                lat_ele(index(i)) = 0
-                sbox_ele(index(i)) = 0
+            if(sorted_index(i) == ele_num) then 
+                r_node(:,:,:,sorted_index(i)) = 0.0_dp
+                type_ele(sorted_index(i)) =''
+                size_ele(sorted_index(i)) = 0
+                lat_ele(sorted_index(i)) = 0
+                sbox_ele(sorted_index(i)) = 0
             else
-                node_num = node_num - ng_node(lat_ele(index(i)))
-                r_node(:,:,:,index(i)) = r_node(:,:,:,ele_num)
-                type_ele(index(i)) = type_ele(ele_num)
-                size_ele(index(i)) = size_ele(ele_num)
-                lat_ele(index(i)) = lat_ele(ele_num)
-                sbox_ele(index(i)) = sbox_ele(ele_num)
+                node_num = node_num - ng_node(lat_ele(sorted_index(i)))
+                r_node(:,:,:,sorted_index(i)) = r_node(:,:,:,ele_num)
+                type_ele(sorted_index(i)) = type_ele(ele_num)
+                size_ele(sorted_index(i)) = size_ele(ele_num)
+                lat_ele(sorted_index(i)) = lat_ele(ele_num)
+                sbox_ele(sorted_index(i)) = sbox_ele(ele_num)
             end if            
             ele_num = ele_num - 1
         end do 
@@ -632,4 +648,4 @@ module elements
     end select
     end subroutine
 
-end module elements
\ No newline at end of file
+end module elements
diff --git a/src/sorts.f90 b/src/sorts.f90
new file mode 100644
index 0000000..b65abc4
--- /dev/null
+++ b/src/sorts.f90
@@ -0,0 +1,854 @@
+! These various sort routines are provided as
+! examples, rather than being fully optimised
+! for the best possible performance. MJR 9/2019
+!
+! So far the collection is:
+!
+!  bubble_sort           (order N**2)
+!  selection_sort        (order N**2)
+!  inline_selection_sort (order N**2)
+!  insertion_sort        (order N**2)
+!  binary_insertion_sort (order N**2)
+!  pair_insertion_sort   (order N**2)
+!  double_insertion_sort (order N**2)
+!  shell_sort            (order N**(4/3)) (possibly)
+!  merge_sort            (order N log N)
+!  alt_merge_sort        (order N log N)
+!  heap_sort             (order N log N)
+!  heap_sort_nr          (order N log N)
+!  smoothsort            (order N log N)
+!  quicksort             (order N log N) (usually)
+!  quicksort_nr          (order N log N) (usually)
+!  dpquicksort           (order N log N) (usually)
+
+! 14/10/19 and 18/10/19 indexing errors in heap sorts corrected
+
+! 4/4/20 smooth sort and non-recursive quicksort added, and some
+!        minor code tidying
+
+
+module sorts
+  use parameters
+  implicit none
+  integer, parameter :: prec=dp, i64=selected_int_kind(15)
+
+  ! array L used by smoothsort
+
+  integer, parameter :: L(0:43) = (/ 1, 1, 3, 5, 9, 15, 25, 41, 67, 109, &
+       177, 287, 465, 753, 1219, 1973, 3193, 5167, 8361, &
+       13529, 21891, 35421, 57313, 92735, 150049, &
+       242785, 392835, 635621, 1028457, 1664079, 2692537, &
+       4356617, 7049155, 11405773, 18454929, 29860703, &
+       48315633, 78176337, 126491971, 204668309, 331160281, &
+       535828591, 866988873, 1402817465 /)
+
+
+  private :: binary_search, merge_arrays, heapify, heapify_nr, &
+       smooth_heapify, smooth_extract, sift_in, interheap_sift, L
+
+contains
+
+  recursive subroutine quicksort(array)
+    real(prec), intent(inout)::array(:)
+    real(prec) :: temp,pivot
+    integer :: i,j,last,left,right
+
+    last=size(array)
+
+    if (last.lt.50) then ! use insertion sort on small arrays
+       do i=2,last
+          temp=array(i)
+          do j=i-1,1,-1
+             if (array(j).le.temp) exit
+             array(j+1)=array(j)
+          enddo
+          array(j+1)=temp
+       enddo
+       return
+    endif
+    ! find median of three pivot
+    ! and place sentinels at first and last elements
+    temp=array(last/2)
+    array(last/2)=array(2)
+    if (temp.gt.array(last)) then
+       array(2)=array(last)
+       array(last)=temp
+    else
+       array(2)=temp
+    endif
+    if (array(1).gt.array(last)) then
+       temp=array(1)
+       array(1)=array(last)
+       array(last)=temp
+    endif
+    if (array(1).gt.array(2)) then
+       temp=array(1)
+       array(1)=array(2)
+       array(2)=temp
+    endif
+    pivot=array(2)
+
+    left=3
+    right=last-1
+    do
+       do while(array(left).lt.pivot)
+          left=left+1
+       enddo
+       do while(array(right).gt.pivot)
+          right=right-1
+       enddo
+       if (left.ge.right) exit
+       temp=array(left)
+       array(left)=array(right)
+       array(right)=temp
+       left=left+1
+       right=right-1
+    enddo
+    if (left.eq.right) left=left+1
+    call quicksort(array(1:left-1))
+    call quicksort(array(left:))
+
+  end subroutine quicksort
+
+  ! This version maintains its own stack, to avoid needing to call
+  ! itself recursively. By always pushing the larger "half" to the
+  ! stack, and moving directly to calculate the smaller "half",
+  ! it can guarantee that the stack needs no more than log_2(N)
+  ! entries
+  subroutine quicksort_nr(array)
+    real(prec), intent(inout)::array(:)
+    real(prec) :: temp,pivot
+    integer :: i,j,left,right,low,high
+    ! If your compiler lacks storage_size(), replace
+    ! storage_size(i) by 64
+    integer :: stack(2,storage_size(i)),stack_ptr
+
+    low=1
+    high=size(array)
+    stack_ptr=1
+
+    do
+
+       if (high-low.lt.50) then ! use insertion sort on small arrays
+          do i=low+1,high
+             temp=array(i)
+             do j=i-1,low,-1
+                if (array(j).le.temp) exit
+                array(j+1)=array(j)
+             enddo
+             array(j+1)=temp
+          enddo
+          ! now pop from stack
+          if (stack_ptr.eq.1) return
+          stack_ptr=stack_ptr-1
+          low=stack(1,stack_ptr)
+          high=stack(2,stack_ptr)
+          cycle
+       endif
+
+       ! find median of three pivot
+       ! and place sentinels at first and last elements
+       temp=array((low+high)/2)
+       array((low+high)/2)=array(low+1)
+       if (temp.gt.array(high)) then
+          array(low+1)=array(high)
+          array(high)=temp
+       else
+          array(low+1)=temp
+       endif
+       if (array(low).gt.array(high)) then
+          temp=array(low)
+          array(low)=array(high)
+          array(high)=temp
+       endif
+       if (array(low).gt.array(low+1)) then
+          temp=array(low)
+          array(low)=array(low+1)
+          array(low+1)=temp
+       endif
+       pivot=array(low+1)
+
+       left=low+2
+       right=high-1
+       do
+          do while(array(left).lt.pivot)
+             left=left+1
+          enddo
+          do while(array(right).gt.pivot)
+             right=right-1
+          enddo
+          if (left.ge.right) exit
+          temp=array(left)
+          array(left)=array(right)
+          array(right)=temp
+          left=left+1
+          right=right-1
+       enddo
+       if (left.eq.right) left=left+1
+       !          call quicksort(array(1:left-1))
+       !          call quicksort(array(left:))
+       if (left.lt.(low+high)/2) then
+          stack(1,stack_ptr)=left
+          stack(2,stack_ptr)=high
+          stack_ptr=stack_ptr+1
+          high=left-1
+       else
+          stack(1,stack_ptr)=low
+          stack(2,stack_ptr)=left-1
+          stack_ptr=stack_ptr+1
+          low=left
+       endif
+
+    enddo
+  end subroutine quicksort_nr
+
+
+  ! dual pivot quicksort
+  recursive subroutine dpquicksort(array)
+    real(prec), intent(inout)::array(:)
+    real(prec) :: temp,p1,p2
+    integer :: i,j,last,l,k,g
+
+    last=size(array)
+
+    if (last.lt.40) then ! use insertion sort on small arrays
+       do i=2,last
+          temp=array(i)
+          do j=i-1,1,-1
+             if (array(j).le.temp) exit
+             array(j+1)=array(j)
+          enddo
+          array(j+1)=temp
+       enddo
+       return
+    endif
+    p1=array(last/3)
+    p2=array(2*last/3)
+    if (p2.lt.p1) then
+       temp=p1
+       p1=p2
+       p2=temp
+    endif
+    array(last/3)=array(1)
+    array(1)=p1
+    array(2*last/3)=array(last)
+    array(last)=p2
+
+    g=last
+    l=2
+    do while (array(l).lt.p1)
+       l=l+1
+    enddo
+    k=l
+
+    do while(k.lt.g)
+       temp=array(k)
+       if (temp.lt.p1) then
+          array(k)=array(l)
+          array(l)=temp
+          l=l+1
+       else if (temp.gt.p2) then
+          do while(array(g-1).gt.p2)
+             g=g-1
+          enddo
+          if (k.ge.g) exit
+          g=g-1
+          if (array(g).lt.p1) then
+             array(k)=array(l)
+             array(l)=array(g)
+             array(g)=temp
+             l=l+1
+          else
+             array(k)=array(g)
+             array(g)=temp
+          endif
+       endif
+       k=k+1
+    enddo
+    if (l.gt.2) then
+       array(1)=array(l-1)
+       array(l-1)=p1
+       call dpquicksort(array(1:l-2))
+    endif
+    call dpquicksort(array(l:g-1))
+    if (g.lt.last) then
+       array(last)=array(g)
+       array(g)=p2
+       call dpquicksort(array(g+1:last))
+    endif
+
+  end subroutine dpquicksort
+
+  function binary_search(array, x, start, end)
+    integer binary_search
+    integer, intent(in) :: start,end
+    real(prec), intent(in) :: x,array(:)
+    integer :: a,b,mid
+
+    a=start
+    b=end
+
+    do
+       if (a.eq.b) then
+          if (array(a).gt.x) then
+             binary_search=a
+             return 
+          else
+             binary_search=a+1
+             return
+          endif
+       endif
+       if (a.gt.b) then
+          binary_search=a
+          return
+       end if
+
+
+       mid=(a+b)/2
+       if (array(mid).lt.x) then
+          a=mid+1
+       else if (array(mid).gt.x) then
+          b=mid-1
+       else
+          binary_search=mid
+          return
+       endif
+
+    end do
+  end function binary_search
+
+  subroutine binary_insertion_sort(array)
+    real(prec), intent(inout) :: array(:)
+    integer :: i,j,pos
+    real(prec) :: x
+
+    do i=2,size(array)
+       x=array(i)
+       pos=binary_search(array,x,1,i-1)
+       do j=i,pos+1,-1
+          array(j)=array(j-1)
+       enddo
+       array(pos)=x
+    enddo
+  end subroutine binary_insertion_sort
+
+  subroutine insertion_sort(array)
+    real(prec), intent(inout) :: array(:)
+    integer :: i,j
+    real(prec) :: temp
+
+    do i=2,size(array)
+       temp=array(i)
+       do j=i-1,1,-1
+          if (array(j).le.temp) exit
+          array(j+1)=array(j)
+       enddo
+       array(j+1)=temp
+    enddo
+  end subroutine insertion_sort
+
+  subroutine double_insertion_sort(array)
+    real(prec), intent(inout) :: array(:)
+    integer :: i,j,left,right,last
+    real(prec) :: temp,p,next
+
+    last=size(array)
+
+    p=0.5*(array(1)+array(last))
+    if (array(1).gt.array(last)) then
+       temp=array(last)
+       array(last)=array(1)
+       array(1)=temp
+    endif
+
+    left=1
+    right=last
+    temp=array(2)
+
+    do i=2,last-1
+       if (temp.lt.p) then
+          do j=left,1,-1
+             if (array(j).le.temp) exit
+             array(j+1)=array(j)
+          enddo
+          array(j+1)=temp
+          temp=array(left+2)
+          left=left+1
+       else
+          next=array(right-1)
+          do j=right,last
+             if (array(j).ge.temp) exit
+             array(j-1)=array(j)
+          enddo
+          array(j-1)=temp
+          temp=next
+          right=right-1
+       endif
+    enddo
+
+  end subroutine double_insertion_sort
+
+
+  subroutine pair_insertion_sort(array)
+    real(prec), intent(inout) :: array(:)
+    integer :: i,j,last
+    real(prec) :: t1,t2
+
+    last=size(array)
+    do i=2,last-1,2
+       t1=min(array(i),array(i+1))
+       t2=max(array(i),array(i+1))
+       j=i-1
+       do while((j.ge.1).and.(array(j).gt.t2))
+          array(j+2)=array(j)
+          j=j-1
+       enddo
+       array(j+2)=t2
+       do while((j.ge.1).and.(array(j).gt.t1))
+          array(j+1)=array(j)
+          j=j-1
+       enddo
+       array(j+1)=t1
+    end do
+
+    if(mod(last,2).eq.0)then
+       t1=array(last)
+       do j=last-1,1,-1
+          if (array(j).le.t1) exit
+          array(j+1)=array(j)
+       end do
+       array(j+1)=t1
+    endif
+
+  end subroutine pair_insertion_sort
+
+
+  subroutine merge_arrays(a,b,out)
+    ! routine used by merge_sort()
+    real(prec), intent(in) :: a(:),b(:)
+    real(prec), intent(out) :: out(:)
+    integer :: ai,bi,oi,i
+
+    ai=1
+    bi=1
+    oi=1
+
+    do
+       if (ai.gt.size(a)) then
+          do i=bi,size(b)
+             out(oi)=b(i)
+             oi=oi+1
+          enddo
+          return
+       endif
+       if (bi.gt.size(b)) then
+          do i=ai,size(a)
+             out(oi)=a(i)
+             oi=oi+1
+          enddo
+          return
+       endif
+
+       if (a(ai).lt.b(bi)) then
+          out(oi)=a(ai)
+          ai=ai+1
+       else
+          out(oi)=b(bi)
+          bi=bi+1
+       endif
+       oi=oi+1
+    enddo
+  end subroutine merge_arrays
+
+  recursive subroutine merge_sort(array)
+    ! With a little skill, the extra copy of "array=temp"
+    ! can be avoided
+    real(prec), intent(inout) :: array(:)
+    real(prec), allocatable :: temp(:)
+    integer :: last
+
+    last=size(array)
+
+    if (last.lt.40) then
+       call insertion_sort(array)
+       return
+    endif
+
+    call merge_sort(array(1:last/2))
+    call merge_sort(array(last/2+1:last))
+    allocate(temp(last))
+    call merge_arrays(array(1:last/2),array(last/2+1:last),temp)
+    array=temp
+    deallocate(temp)
+  end subroutine merge_sort
+
+  recursive subroutine alt_merge_sort(array)
+    ! Alternate merges between array and temp
+    ! to avoid copy
+    real(prec), intent(inout) :: array(:)
+    real(prec), allocatable :: temp(:)
+    integer :: last
+
+    last=size(array)
+
+    if (last.lt.40) then
+       call insertion_sort(array)
+       return
+    endif
+
+    call alt_merge_sort(array(1:last/4))
+    call alt_merge_sort(array(last/4+1:last/2))
+    call alt_merge_sort(array(last/2+1:last/2+last/4))
+    call alt_merge_sort(array(last/2+last/4+1:last))
+
+    allocate(temp(last))
+    call merge_arrays(array(1:last/4),array(last/4+1:last/2),temp(1:last/2))
+    call merge_arrays(array(last/2+1:last/2+last/4),  &
+         array(last/2+last/4+1:last),temp(last/2+1:last))
+    call merge_arrays(temp(1:last/2),temp(last/2+1:last),array)
+    deallocate(temp)
+  end subroutine alt_merge_sort
+
+  subroutine selection_sort(array)
+    ! This version uses Fortran's minloc intrinsic
+    real(prec), intent(inout) :: array(:)
+    real(prec) :: temp
+    integer :: i,j
+
+    do i=1,size(array)-1
+       j=minloc(array(i:),1)+i-1
+       temp=array(i)
+       array(i)=array(j)
+       array(j)=temp
+    enddo
+
+  end subroutine selection_sort
+
+  subroutine inline_selection_sort(array)
+    ! This version does not use Fortran's minloc intrinsic
+    real(prec), intent(inout) :: array(:)
+    real(prec) :: temp
+    integer :: i,j,k
+
+    do i=1,size(array)-1
+       temp=array(i)
+       j=i
+       do k=i+1,size(array)
+          if (array(k).lt.temp) then
+             temp=array(k)
+             j=k
+          endif
+       enddo
+       array(j)=array(i)
+       array(i)=temp
+    enddo
+
+  end subroutine inline_selection_sort
+
+  subroutine shell_sort(array)
+    ! There are many choices of intervals
+    ! The sequence (3**k-1)/2 used below is given by Knuth
+    ! The sort relies on the insertion sort of nearly sorted
+    ! data being fast (i.e. not order N**2)
+    real(prec), intent(inout) :: array(:)
+    integer :: i,interval,last
+
+    interval=1
+    last=size(array)
+    do while(interval<last/3)
+       interval=3*interval+1
+    enddo
+
+    do while(interval>0)
+       do i=1,interval
+          ! The following syntax passes a subarray with elements
+          ! array(i),array(i+interval),array(i+2*interval),...
+          call insertion_sort(array(i::interval))
+       enddo
+       interval=(interval-1)/3
+    enddo
+
+  end subroutine shell_sort
+
+  subroutine bubble_sort(array)
+
+    real(prec), intent(inout) :: array(:)
+    real(prec) :: temp
+    integer :: i,j,last
+
+    last=size(array)
+    do i=last-1,1,-1
+       do j=1,i
+          if (array(j+1).lt.array(j)) then
+             temp=array(j+1)
+             array(j+1)=array(j)
+             array(j)=temp
+          endif
+       enddo
+    enddo
+
+  end subroutine bubble_sort
+
+  subroutine heap_sort(array)
+    real(prec), intent(inout) :: array(:)
+    real(prec) :: temp
+    integer :: i,last
+
+    last=size(array)
+
+    ! Build heap
+    do i=last/2,1,-1
+       call heapify(array,i)
+    enddo
+
+    ! Unpick heap
+    do i=last,2,-1
+       temp=array(1)
+       array(1)=array(i)
+       array(i)=temp
+       call heapify(array(1:i-1),1)
+    end do
+
+  end subroutine heap_sort
+
+  recursive subroutine heapify(array,root)
+    real(prec), intent(inout) :: array(:)
+    real(prec) :: temp
+    integer :: left,right,root,last,largest
+
+    last=size(array)
+    left=2*root
+    right=left+1
+    largest=root
+
+    if (left.le.last) then
+       if(array(left).gt.array(largest)) largest=left
+    endif
+
+    if (right.le.last) then
+       if(array(right).gt.array(largest)) largest=right
+    endif
+
+    if (largest.ne.root) then
+       temp=array(root)
+       array(root)=array(largest)
+       array(largest)=temp
+       call heapify(array,largest)
+    end if
+
+  end subroutine heapify
+
+
+  subroutine heap_sort_nr(array)
+    ! The same but with a non-recursive heapify
+    real(prec), intent(inout) :: array(:)
+    real(prec) :: temp
+    integer :: i,last
+
+    last=size(array)
+
+    ! Build heap
+    do i=last/2,1,-1
+       call heapify_nr(array,i)
+    enddo
+
+    ! Unpick heap
+    do i=last,2,-1
+       temp=array(1)
+       array(1)=array(i)
+       array(i)=temp
+       call heapify_nr(array(1:i-1),1)
+    end do
+
+  end subroutine heap_sort_nr
+
+  subroutine heapify_nr(array,i)
+    ! A non-recursive heapify
+    real(prec), intent(inout) :: array(:)
+    real(prec) :: root_val
+    integer :: i,left,right,root,last,largest
+
+    last=size(array)
+    root=i
+    left=2*root
+    root_val=array(root)
+    largest=root
+
+    do while(left.le.last)
+       right=left+1
+
+       if (left.le.last) then
+          if(array(left).gt.array(largest)) largest=left
+       endif
+
+       if (right.le.last) then
+          if(array(right).gt.array(largest)) largest=right
+       endif
+
+       if (largest.eq.root) exit
+
+       array(root)=array(largest)
+       array(largest)=root_val
+       root=largest
+       left=2*root
+    enddo
+
+    !    array(largest)=root_val
+
+  end subroutine heapify_nr
+
+  ! Smoothsort is introduced by Edsger Dijkstra
+  ! https://doi.org/10.1016/0167-6423(82)90016-8
+  !
+  ! This Fortran version owes much to https://www.keithschwarz.com/smoothsort/
+  ! and the C version given by Martin Knoblauch Revuelta at
+  ! http://code.google.com/archive/p/combsortcs2p-and-other-sorting-algorithms
+
+  subroutine smoothsort(array)
+    real(prec), intent(inout) :: array(:)
+    integer(i64) :: mask
+    integer :: offset
+
+    if (size(array).le.1) return
+
+    call smooth_heapify(array,mask,offset)
+
+    call smooth_extract(array,mask,offset)
+  end subroutine smoothsort
+
+  subroutine smooth_heapify(array,mask,offset)
+    real(prec), intent(inout) :: array(0:)
+    integer(i64), intent(out) :: mask
+    integer, intent(out) :: offset
+    integer :: i,last
+    logical :: wbf
+
+    mask=1
+    offset=1
+    last=size(array)
+
+    do i=1,last-1
+       if (iand(mask,2_i64).eq.2) then
+          mask=ior(ishft(mask,-2),1_i64)
+          offset=offset+2
+       else if (offset.eq.1) then
+          mask=ior(ishft(mask,1),1_i64)
+          offset=0
+       else
+          mask=ior(ishft(mask,offset-1),1_i64)
+          offset=1
+       endif
+       wbf=(((iand(mask,2_i64).eq.2).and.(i+1.lt.last)).or. &
+            ((offset.gt.0).and.(1+i+L(offset-1).lt.last)))
+       if (wbf) then
+          call sift_in(array,i,offset)
+       else
+          call interheap_sift(array,i,mask,offset)
+       endif
+    enddo
+
+  end subroutine smooth_heapify
+  subroutine smooth_extract(array,mask,offset)
+    real(prec), intent(inout) :: array(0:)
+    integer :: i,j,last,ch(2)
+    integer(i64),value :: mask
+    integer, value :: offset
+
+    last=size(array)
+
+    do i=last-1,2,-1
+       if (offset.lt.2) then
+          do
+             mask=ishft(mask,-1)
+             offset=offset+1
+             if (iand(mask,1_i64).eq.1) exit
+          enddo
+       else
+          ch(2)=i-1
+          ch(1)=ch(2)-L(offset-2)
+          mask=iand(mask,not(1_i64))
+          do j=1,2
+             mask=ior(ishft(mask,1),1_i64)
+             offset=offset-1
+             call interheap_sift(array,ch(j),mask,offset)
+          end do
+       end if
+    end do
+
+
+  end subroutine smooth_extract
+
+
+  subroutine sift_in(array,root,sz)
+    real(prec), intent(inout) :: array(0:)
+    real(prec) :: tmp
+    integer :: left,right,next
+    integer,value :: root,sz
+
+    if (sz.lt.2) return
+
+    tmp=array(root)
+
+    do
+       right=root-1
+       left=right-L(sz-2)
+
+       if (array(right).lt.array(left)) then
+          next=left
+          sz=sz-1
+       else
+          next=right
+          sz=sz-2
+       end if
+
+       if (array(next).le.tmp) exit
+
+       array(root)=array(next)
+
+       root=next
+       if (sz.le.1) exit
+    end do
+
+    array(root)=tmp
+
+  end subroutine sift_in
+
+  subroutine interheap_sift(array,root,mask,offset)
+    real(prec), intent(inout) :: array(0:)
+    real(prec) :: tmp,mx
+    integer :: left,right,next
+    integer(i64),value :: mask
+    integer, value :: root,offset
+
+    tmp=array(root)
+
+    do while (mask.ne.1)
+       mx=tmp
+
+       if (offset.gt.1) then
+          right=root-1
+          left=right-L(offset-2)
+          mx=max(mx,array(left),array(right))
+       end if
+
+       next=root-L(offset)
+
+       if (array(next).le.mx) exit
+
+       array(root)=array(next)
+       root=next
+
+       do
+          mask=ishft(mask,-1)
+          offset=offset+1
+          if (iand(mask,1_i64).eq.1) exit
+       end do
+    end do
+
+
+    array(root)=tmp
+    call sift_in(array,root,offset)
+  end subroutine interheap_sift
+
+
+end module sorts
diff --git a/src/subroutines.f90 b/src/subroutines.f90
index b40b344..ece2294 100644
--- a/src/subroutines.f90
+++ b/src/subroutines.f90
@@ -173,52 +173,6 @@ module subroutines
     end subroutine parse_ori_vec
     
 
-    subroutine heapsort(a)
-     
-       integer, intent(in out) :: a(0:)
-       integer :: start, n, bottom
-       real :: temp
-     
-       n = size(a)
-       do start = (n - 2) / 2, 0, -1
-         call siftdown(a, start, n);
-       end do
-     
-       do bottom = n - 1, 1, -1
-         temp = a(0)
-         a(0) = a(bottom)
-         a(bottom) = temp;
-         call siftdown(a, 0, bottom)
-       end do
-     
-    end subroutine heapsort
-     
-    subroutine siftdown(a, start, bottom)
-     
-      integer, intent(in out) :: a(0:)
-      integer, intent(in) :: start, bottom
-      integer :: child, root
-      real :: temp
-     
-      root = start
-      do while(root*2 + 1 < bottom)
-        child = root * 2 + 1
-     
-        if (child + 1 < bottom) then
-          if (a(child) < a(child+1)) child = child + 1
-        end if
-     
-        if (a(root) < a(child)) then
-          temp = a(child)
-          a(child) = a (root)
-          a(root) = temp
-          root = child
-        else
-          return
-        end if  
-      end do      
-     
-    end subroutine siftdown
 
 
     subroutine apply_periodic(r)
@@ -344,4 +298,4 @@ module subroutines
         return
     end subroutine check_right_ortho
 
-end module subroutines
\ No newline at end of file
+end module subroutines