c     Produce an index indx to the array arr(n) so that when accessed via
c     the index the elements are in ascending numerical order.
c     
      SUBROUTINE  indexxx(n,arr,indx) 
c     
      IMPLICIT none
      INTEGER n,NSWITCH
      PARAMETER(NSWITCH=60)
      REAL arr(n)
      INTEGER indx(n)
c     
      if (n.gt.NSWITCH) then    !if n is large use the Numerical recipes routine
         call indexx(n,arr,indx) !heapsort from numerical recipes 
      else                      !simple shell sort which is faster for small n but scales as n^3/2 in worst case
         call indexsh(n,arr,indx) !rather than the n log(n) of heapsort
      end if
c     
      return
      end
c     
*****************************************************************************
*****************************************************************************
c     Indexing or array arr(n) using shell sort
c     
      SUBROUTINE  indexsh(n,arr,indx)
c     
      IMPLICIT none
      INTEGER n,m
      REAL arr(n),aln2i,tiny
      PARAMETER(TINY=1.0e-05,ALN2I=1.0/0.69314718)
      INTEGER indx(n),nn,lognb2,l,i,tindx,j,k
c     
      do i=1,n
         indx(i)=i
      end do
      if(n.ge.2) then
         lognb2=int(log(float(n))*ALN2I+TINY)
         m=n
         do nn=1,lognb2
            m=m/2
            k=n-m
            do j=1,k
               i=j
 3             continue
               l=i+m
               if(arr(indx(l)).lt.arr(indx(i))) then
                  tindx=indx(i)
                  indx(i)=indx(l)
                  indx(l)=tindx
                  i=i-m
                  if (i.ge.1) goto 3
               end if
            end do
         end do
      end if
c     
      return
      END
c     
*****************************************************************************
*****************************************************************************
c     Numerical Recipes heapsort
	SUBROUTINE INDEXX(N,ARRIN,INDX)
        IMPLICIT NONE
        INTEGER N
        INTEGER INDX(N),J,L,I,IR,INDXT
	REAL ARRIN(N),Q
	DO J=1,N
		INDX(J)=J
	ENDDO
        if(n.le.1) return
	L=N/2+1
	IR=N
10	CONTINUE
		IF(L.GT.1)THEN
			L=L-1
			INDXT=INDX(L)
			Q=ARRIN(INDXT)
		ELSE
			INDXT=INDX(IR)
			Q=ARRIN(INDXT)
			INDX(IR)=INDX(1)
			IR=IR-1
			IF(IR.EQ.1)THEN
				INDX(1)=INDXT
				RETURN
			ENDIF
		ENDIF
		I=L
		J=L+L
20		IF(J.LE.IR)THEN
			IF(J.LT.IR)THEN
				IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
			ENDIF
			IF(Q.LT.ARRIN(INDX(J)))THEN
				INDX(I)=INDX(J)
				I=J
				J=J+J
			ELSE
				J=IR+1
			ENDIF
		GO TO 20
		ENDIF
		INDX(I)=INDXT
	GO TO 10
	END