C     this routine will return the integer array ORDER containing the
C     indexes of the elements of X arranged in order of increasing X
C     values.  It uses a bubble sort, which is quite efficient if X is
C     already nearly arranged in increasing order.

      SUBROUTINE BUBBLE(X,ORDER,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 X(*)
      LOGICAL STUCK
      INTEGER ORDER(*), BUB
      DO 10 I=1,N
10       ORDER(I)=I
      BUB=N
1     CALL SINK(BUB,X,ORDER,N,STUCK)
      IF (STUCK) THEN
         BUB=BUB-1
         IF(BUB.EQ.1) THEN
            RETURN
         ELSE
            GO TO 1
         ENDIF
      ELSE
         IF (BUB.LT.N) BUB=BUB+1
         GO TO 1
      ENDIF
      END

      SUBROUTINE SINK(BUB,X,ORDER,N,STUCK)
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER ORDER,BUB,BUBLOC
      LOGICAL STUCK
      DIMENSION X(*),ORDER(*)
      BUBLOC=BUB
10    CONTINUE
      IF (X(ORDER(BUBLOC)).LT.X(ORDER(BUBLOC-1))) THEN
         IORD=ORDER(BUBLOC)
         ORDER(BUBLOC)=ORDER(BUBLOC-1)
         ORDER(BUBLOC-1)=IORD
         BUBLOC=BUBLOC-1
         IF (BUBLOC.EQ.1) THEN
            STUCK=BUBLOC.EQ.BUB
            RETURN
         ELSE
            GO TO 10
         ENDIF
      ELSE
         STUCK=BUBLOC.EQ.BUB
         RETURN
      ENDIF
      END
