SUBROUTINE SORT (ARRAY, IINDEX, N, IRANK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C     **************************************************
C     Sort an array with the HEAPSORT algorithm 
C     (cf. W. H. Press et al., "Numerical Recipes", 
C     Cambridge(1986), p. 229 - 235)
C     **************************************************
      DIMENSION ARRAY(N), IINDEX(N), IRANK(N)
C     Initialize the index array
      DO J = 1,N
      IINDEX(J) = J
      ENDDO
C     Index an array ARRAY of length N, i. e. outputs the 
C     array IINDEX such that ARRAY(IINDEX(J)) is in ascending 
C     order for J=1,2,...,N.  ARRAY itself is not changed.
      IF (N.EQ.1) GOTO 30
      L = N/2+1
      IR = N
 10   IF (L.GT.1) THEN
      L = L-1
      INDEXT = IINDEX(L)
      Q = ARRAY(INDEXT)
      ELSE
      INDEXT = IINDEX(IR)
      Q = ARRAY(INDEXT)
      IINDEX(IR) = IINDEX(1)
      IR = IR-1
      IF (IR.EQ.1) THEN
      IINDEX(1) = INDEXT
      GOTO 30
      ENDIF
      ENDIF
      I = L
      J = L+L
 20   IF (J.LE.IR) THEN
      IF (J.LT.IR) THEN
      IF (ARRAY(IINDEX(J)).LT.ARRAY(IINDEX(J+1))) J=J+1
      ENDIF
      IF (Q.LT.ARRAY(IINDEX(J))) THEN
      IINDEX(I) = IINDEX(J)
      I = J
      J = J+J
      ELSE
      J = IR+1
      ENDIF
      GOTO 20
      ENDIF
      IINDEX(I) = INDEXT
      GOTO 10
 30   DO J=1,N
      IRANK(IINDEX(J)) = J
      ENDDO
      RETURN
      END