SUBROUTINE indexx_dp(arr,index)
  USE nrtype; USE nrutil, ONLY : arth,assert_eq,nrerror,swap
  IMPLICIT NONE
  REAL(DP), DIMENSION(:), INTENT(IN) :: arr
  INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: index
  INTEGER(I4B), PARAMETER :: NN=15, NSTACK=50
  REAL(DP) :: a
  INTEGER(I4B) :: n,k,i,j,indext,jstack,l,r
  INTEGER(I4B), DIMENSION(NSTACK) :: istack
  n=assert_eq(size(index),size(arr),'indexx_dp')
  index=arth(1,1,n)
  jstack=0
  l=1
  r=n
  do
     if (r-l < NN) then
        do j=l+1,r
           indext=index(j)
           a=arr(indext)
           do i=j-1,l,-1
              if (arr(index(i)) <= a) exit
              index(i+1)=index(i)
           end do
           index(i+1)=indext
        end do
        if (jstack == 0) RETURN
        r=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
     else
        k=(l+r)/2
        call swap(index(k),index(l+1))
        call icomp_xchg(index(l),index(r))
        call icomp_xchg(index(l+1),index(r))
        call icomp_xchg(index(l),index(l+1))
        i=l+1
        j=r
        indext=index(l+1)
        a=arr(indext)
        do
           do
              i=i+1
              if (arr(index(i)) >= a) exit
           end do
           do
              j=j-1
              if (arr(index(j)) <= a) exit
           end do
           if (j < i) exit
           call swap(index(i),index(j))
        end do
        index(l+1)=index(j)
        index(j)=indext
        jstack=jstack+2
        if (jstack > NSTACK) call nrerror('indexx: NSTACK too small')
        if (r-i+1 >= j-l) then
           istack(jstack)=r
           istack(jstack-1)=i
           r=j-1
        else
           istack(jstack)=j-1
           istack(jstack-1)=l
           l=i
        end if
     end if
  end do
CONTAINS
  !BL
  SUBROUTINE icomp_xchg(i,j)
    INTEGER(I4B), INTENT(INOUT) :: i,j
    INTEGER(I4B) :: swp
    if (arr(j) < arr(i)) then
       swp=i
       i=j
       j=swp
    end if
  END SUBROUTINE icomp_xchg
END SUBROUTINE indexx_dp
 
SUBROUTINE indexx_sp(arr,index)
  USE nrtype; USE nrutil, ONLY : arth,assert_eq,nrerror,swap
  IMPLICIT NONE
  REAL(SP), DIMENSION(:), INTENT(IN) :: arr
  INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: index
  INTEGER(I4B), PARAMETER :: NN=15, NSTACK=50
  REAL(SP) :: a
  INTEGER(I4B) :: n,k,i,j,indext,jstack,l,r
  INTEGER(I4B), DIMENSION(NSTACK) :: istack
  n=assert_eq(size(index),size(arr),'indexx_sp')
  index=arth(1,1,n)
  jstack=0
  l=1
  r=n
  do
     if (r-l < NN) then
        do j=l+1,r
           indext=index(j)
           a=arr(indext)
           do i=j-1,l,-1
              if (arr(index(i)) <= a) exit
              index(i+1)=index(i)
           end do
           index(i+1)=indext
        end do
        if (jstack == 0) RETURN
        r=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
     else
        k=(l+r)/2
        call swap(index(k),index(l+1))
        call icomp_xchg(index(l),index(r))
        call icomp_xchg(index(l+1),index(r))
        call icomp_xchg(index(l),index(l+1))
        i=l+1
        j=r
        indext=index(l+1)
        a=arr(indext)
        do
           do
              i=i+1
              if (arr(index(i)) >= a) exit
           end do
           do
              j=j-1
              if (arr(index(j)) <= a) exit
           end do
           if (j < i) exit
           call swap(index(i),index(j))
        end do
        index(l+1)=index(j)
        index(j)=indext
        jstack=jstack+2
        if (jstack > NSTACK) call nrerror('indexx: NSTACK too small')
        if (r-i+1 >= j-l) then
           istack(jstack)=r
           istack(jstack-1)=i
           r=j-1
        else
           istack(jstack)=j-1
           istack(jstack-1)=l
           l=i
        end if
     end if
  end do
CONTAINS
  !BL
  SUBROUTINE icomp_xchg(i,j)
    INTEGER(I4B), INTENT(INOUT) :: i,j
    INTEGER(I4B) :: swp
    if (arr(j) < arr(i)) then
       swp=i
       i=j
       j=swp
    end if
  END SUBROUTINE icomp_xchg
END SUBROUTINE indexx_sp

SUBROUTINE indexx_i4b(iarr,index)
  USE nrtype; USE nrutil, ONLY : arth,assert_eq,nrerror,swap
  IMPLICIT NONE
  INTEGER(I4B), DIMENSION(:), INTENT(IN) :: iarr
  INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: index
  INTEGER(I4B), PARAMETER :: NN=15, NSTACK=50
  INTEGER(I4B) :: a
  INTEGER(I4B) :: n,k,i,j,indext,jstack,l,r
  INTEGER(I4B), DIMENSION(NSTACK) :: istack
  n=assert_eq(size(index),size(iarr),'indexx_sp')
  index=arth(1,1,n)
  jstack=0
  l=1
  r=n
  do
     if (r-l < NN) then
        do j=l+1,r
           indext=index(j)
           a=iarr(indext)
           do i=j-1,1,-1
              if (iarr(index(i)) <= a) exit
              index(i+1)=index(i)
           end do
           index(i+1)=indext
        end do
        if (jstack == 0) RETURN
        r=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
     else
        k=(l+r)/2
        call swap(index(k),index(l+1))
        call icomp_xchg(index(l),index(r))
        call icomp_xchg(index(l+1),index(r))
        call icomp_xchg(index(l),index(l+1))
        i=l+1
        j=r
        indext=index(l+1)
        a=iarr(indext)
        do
           do
              i=i+1
              if (iarr(index(i)) >= a) exit
           end do
           do
              j=j-1
              if (iarr(index(j)) <= a) exit
           end do
           if (j < i) exit
           call swap(index(i),index(j))
        end do
        index(l+1)=index(j)
        index(j)=indext
        jstack=jstack+2
        if (jstack > NSTACK) call nrerror('indexx: NSTACK too small')
        if (r-i+1 >= j-l) then
           istack(jstack)=r
           istack(jstack-1)=i
           r=j-1
        else
           istack(jstack)=j-1
           istack(jstack-1)=l
           l=i
        end if
     end if
  end do
CONTAINS
  !BL
  SUBROUTINE icomp_xchg(i,j)
    INTEGER(I4B), INTENT(INOUT) :: i,j
    INTEGER(I4B) :: swp
    if (iarr(j) < iarr(i)) then
       swp=i
       i=j
       j=swp
    end if
  END SUBROUTINE icomp_xchg
END SUBROUTINE indexx_i4b
