module tsort
! This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/TestSortWSB.f90

contains

Function Sort(z,n)
! Sort is returned as the vector which orders z(1), z(2),...,z(n)
! in ascending order.
  implicit none
  integer, INTENT(IN) :: n
  real, dimension(n),INTENT(INOUT) :: z
  integer, dimension(n) :: Sort
  integer, allocatable, dimension(:) :: iz
  integer :: i, ierr
  allocate (iz(n),stat=ierr)
  if (ierr.ne.0) then
     print *, "??ierr=",ierr," unable to allocate iz."
     stop
  endif
  Do i=1,n
   iz(i)=i
  end do
! Quick_Sort is Walt Grainerd's code modified by Alan Miller.
! http://members.ozemail.com.au/~milleraj/misc/qsort.f90
  call Quick_Sort(z,iz)
  Sort=z
  deallocate(iz)
End Function Sort

RECURSIVE SUBROUTINE quick_sort(list, order)

! Quick sort routine from:
! Brainerd, W.S., Goldberg, C.H. & Adams, J.C. (1990) "Programmer's Guide to
! Fortran 90", McGraw-Hill  ISBN 0-07-000248-7, pages 149-150.
! Modified by Alan Miller to include an associated integer array which gives
! the positions of the elements in the original order.

IMPLICIT NONE
REAL, DIMENSION (:), INTENT(IN OUT)  :: list
INTEGER, DIMENSION (:), INTENT(OUT)  :: order

! Local variable
INTEGER :: i

DO i = 1, SIZE(list)
  order(i) = i
END DO

CALL quick_sort_1(1, SIZE(list))

CONTAINS

RECURSIVE SUBROUTINE quick_sort_1(left_end, right_end)

INTEGER, INTENT(IN) :: left_end, right_end

!     Local variables
INTEGER             :: i, j, itemp
REAL                :: reference, temp
INTEGER, PARAMETER  :: max_simple_sort_size = 6

IF (right_end < left_end + max_simple_sort_size) THEN
  ! Use interchange sort for small lists
  CALL interchange_sort(left_end, right_end)

ELSE
  ! Use partition ("quick") sort
  reference = list((left_end + right_end)/2)
  i = left_end - 1; j = right_end + 1

  DO
    ! Scan list from left end until element >= reference is found
    DO
      i = i + 1
      IF (list(i) >= reference) EXIT
    END DO
    ! Scan list from right end until element <= reference is found
    DO
      j = j - 1
      IF (list(j) <= reference) EXIT
    END DO


    IF (i < j) THEN
      ! Swap two out-of-order elements
      temp = list(i); list(i) = list(j); list(j) = temp
      itemp = order(i); order(i) = order(j); order(j) = itemp
    ELSE IF (i == j) THEN
      i = i + 1
      EXIT
    ELSE
      EXIT
    END IF
  END DO

  IF (left_end < j) CALL quick_sort_1(left_end, j)
  IF (i < right_end) CALL quick_sort_1(i, right_end)
END IF

END SUBROUTINE quick_sort_1


SUBROUTINE interchange_sort(left_end, right_end)

INTEGER, INTENT(IN) :: left_end, right_end

!     Local variables
INTEGER             :: i, j, itemp
REAL                :: temp

DO i = left_end, right_end - 1
  DO j = i+1, right_end
    IF (list(i) > list(j)) THEN
      temp = list(i); list(i) = list(j); list(j) = temp
      itemp = order(i); order(i) = order(j); order(j) = itemp
    END IF
  END DO
END DO

END SUBROUTINE interchange_sort

END SUBROUTINE quick_sort
end module tsort

program testsort
use tsort
real ,dimension (10) :: x
real ,dimension (10) :: y
data x/12.0,4.0,3.0,6.0,9.0,7.0,5.0,8.0,11.0,10.0/
print *, "Original Order: ", x
y=Sort(x,10)
print *, "Sorted Order: ",y
end program testsort