r_quicksort Subroutine

  pure subroutine r_quicksort (arr_len,key)
    implicit none
    !> Size of arrays
    integer, intent(in) :: arr_len

    !> Arrays to be sorted
    real, intent(inout), dimension(1:arr_len) :: key
    integer, parameter :: insertion_limit=16 !< Lists smaller than this use insertion sort

    !Internals
    integer, dimension(:), allocatable :: stack_inds
    integer :: i, j, right_lim, left_lim, stack_count
    real :: ak !Used to hold selected elements
    integer :: k !Temp used for middle of range
    logical :: partitioned!, done_exit

    !Intialise vars
    allocate(stack_inds(max(arr_len/2,200)))
    stack_count=0 !How many sections/partitions
    left_lim=1 !Left index of current section
    right_lim=arr_len !Right index of current section

    !Infinite loop --> explicit return forced when list sorted
    do while(.true.)
       if(right_lim-left_lim.lt.insertion_limit)then
          !Insertion sort : Loop over current section
          call insertsort(1+right_lim-left_lim,key(left_lim:right_lim))

          !Check if work remaining

          !If stack is empty then done sorting so return
          if(stack_count.eq.0)then
             deallocate(stack_inds)
             return
          endif

          !If not get next section limits from stack
          right_lim=stack_inds(stack_count)
          left_lim=stack_inds(stack_count-1)

          !Reduce stack count
          stack_count=stack_count-2
       else
          !What index are we looking at. Written as left
          !plus half_range to avoid potential overflow if
          !written as (left+right)/2
          k = left_lim + (right_lim - left_lim) / 2

          !Swap values | KEY
          call swap_elem(key,k,left_lim+1)

          !Order so key(left_lim)<=key(left_lim+1)<=key(right_lim)
          if(key(left_lim).gt.key(right_lim)) then
             call swap_elem(key,left_lim,right_lim)
          endif
          if(key(left_lim+1).gt.key(right_lim)) then
             call swap_elem(key,left_lim+1,right_lim)
          endif
          if(key(left_lim).gt.key(left_lim+1)) then
             call swap_elem(key,left_lim,left_lim+1)
          endif

          !Index of left limit of partition
          i=left_lim+1
          !Index of right limit of partition
          j=right_lim

          !Get partition values
          ak=key(left_lim+1)

          !Partition elements
          partitioned=.false.
          do while(.not.partitioned)
             !Search from left for first element bigger than ak
             i=i+1
             do while(key(i).lt.ak)
                i=i+1
             end do

             !Search from right for first element smaller than ak
             j=j-1
             do while(key(j).gt.ak)
                j=j-1
             end do

             !If first smaller element before first larger element then
             !partioning has finished --> Can now sort
             !Otherwise swap elements and try again
             if(j.lt.i)then
                !Indicate partioning complete
                partitioned=.true.

                !Insert partioning element in correct place
                key(left_lim+1)=key(j) ; key(j)=ak
             else
                call swap_elem(key,i,j)
             endif
          end do

          !Update stack counter
          stack_count=stack_count+2

          !At this point i is the start of the left section, which goes up to right_lim
          !and the right section goes from left_lim up to j
          !Update stack_inds values to record section done
          !If left section is larger push to stack and do right
          !otherwise push right to stack and do left
          if(right_lim-i+1.ge.j-left_lim) then
             !Store detail of left section
             stack_inds(stack_count)=right_lim
             stack_inds(stack_count-1)=i

             !Set right lim to one before start of left section
             right_lim=j-1
          else
             !Store details of right section
             stack_inds(stack_count)=j-1
             stack_inds(stack_count-1)=left_lim

             !Set left lim to one after end of right stack
             left_lim=i
          endif
       end if
    enddo
  end subroutine r_quicksort