module quicksort

  ! from Numerical Recipies in Fortran 77, pp. 326-327
  ! http://lib-www.lanl.gov/numerical/

contains
  
  SUBROUTINE sort2(n,arr,brr)
    INTEGER n,M,NSTACK 
    REAL arr(n),brr(n)
    PARAMETER (M=7,NSTACK=50)
    ! Sorts an array arr(1:n) into ascending order using Quicksort, 
    ! while making the corresponding rearrangement of the array brr(1:n).
    INTEGER i,ir,j,jstack,k,l,istack(NSTACK) 
    REAL a,b,temp 
    jstack=0 
    l=1 
    ir=n 
1   if(ir-l.lt.M)then ! Insertion sort when subarray small enough.
       do j=l+1,ir 
          a=arr(j) 
          b=brr(j)
          do i=j-1,l,-1 
             if(arr(i).le.a)goto 2 
             arr(i+1)=arr(i) 
             brr(i+1)=brr(i)
          enddo
          i=l-1
2         arr(i+1)=a 
          brr(i+1)=b 
       enddo
       if(jstack.eq.0)return 
       ir=istack(jstack) ! Pop stack and begin a new round of partitioning. 
       l=istack(jstack-1) 
       jstack=jstack-2 
    else 
       ! Choose median of left, center and right elements as partitioning element a. 
       ! Also rearrange so that a(l) <= a(l+1) <= a(ir).
       k=(l+ir)/2 
       temp=arr(k) 
       arr(k)=arr(l+1) 
       arr(l+1)=temp 
       temp=brr(k) 
       brr(k)=brr(l+1) 
       brr(l+1)=temp 
       if(arr(l).gt.arr(ir))then 
          temp=arr(l) 
          arr(l)=arr(ir) 
          arr(ir)=temp 
          temp=brr(l) 
          brr(l)=brr(ir) 
          brr(ir)=temp 
       endif
       if(arr(l+1).gt.arr(ir))then 
          temp=arr(l+1) 
          arr(l+1)=arr(ir) 
          arr(ir)=temp 
          temp=brr(l+1) 
          brr(l+1)=brr(ir) 
          brr(ir)=temp 
       endif
       if(arr(l).gt.arr(l+1))then 
          temp=arr(l) 
          arr(l)=arr(l+1) 
          arr(l+1)=temp 
          temp=brr(l) 
          brr(l)=brr(l+1) 
          brr(l+1)=temp 
       endif
       i=l+1 ! Initialize pointers for partitioning.
       j=ir 
       a=arr(l+1) ! Partitioning element.
       b=brr(l+1)
3      continue ! Beginning of innermost loop. 
         i=i+1 ! Scan up to find element > a.
       if(arr(i).lt.a)goto 3 
4      continue 
         j=j-1 ! Scan down to find element < a. 
       if(arr(j).gt.a)goto 4 
       if(j.lt.i)goto 5 ! Pointers crossed. Exit with partitioning complete.
       temp=arr(i) ! Exchange elements of both arrays.
       arr(i)=arr(j) 
       arr(j)=temp 
       temp=brr(i) 
       brr(i)=brr(j) 
       brr(j)=temp 
       goto 3 ! End of innermost loop. 
5      arr(l+1)=arr(j) ! Insert partitioning element in both arrays. 
       arr(j)=a 
       brr(l+1)=brr(j) 
       brr(j)=b 
       jstack=jstack+2
       ! Push pointers to larger subarray on stack, process smaller subarray immediately.
       if(jstack.gt.NSTACK)pause 'NSTACK too small in sort2'
       if(ir-i+1.ge.j-l)then
          istack(jstack)=ir 
          istack(jstack-1)=i 
          ir=j-1 
       else 
          istack(jstack)=j-1 
          istack(jstack-1)=l 
          l=i 
       endif
    endif
    goto 1 
  END SUBROUTINE sort2

end module quicksort