      SUBROUTINE DQPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX)
C     ..................................................................
C
C 1.        DQPSRT
C           ORDERING ROUTINE
C              STANDARD FORTRAN SUBROUTINE
C              DOUBLE PRECISION VERSION
C
C 2.        PURPOSE
C              THIS ROUTINE MAINTAINS THE DESCENDING ORDERING
C              IN THE LIST OF THE LOCAL ERROR ESTIMATES RESULTING FROM
C              THE INTERVAL SUBDIVISION PROCESS. AT EACH CALL TWO ERROR
C              ESTIMATES ARE INSERTED USING THE SEQUENTIAL SEARCH
C              METHOD, TOP-DOWN FOR THE LARGEST ERROR ESTIMATE
C              AND BOTTOM-UP FOR THE SMALLEST ERROR ESTIMATE.
C
C 3.        CALLING SEQUENCE
C              CALL DQPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX)
C
C           PARAMETERS (MEANING AT OUTPUT)
C              LIMIT  - INTEGER
C                       MAXIMUM NUMBER OF ERROR ESTIMATES THE LIST
C                       CAN CONTAIN
C
C              LAST   - INTEGER
C                       NUMBER OF ERROR ESTIMATES CURRENTLY
C                       IN THE LIST
C
C              MAXERR - INTEGER
C                       MAXERR POINTS TO THE NRMAX-TH LARGEST ERROR
C                       ESTIMATE CURRENTLY IN THE LIST
C
C              ERMAX  - DOUBLE PRECISION
C                       NRMAX-TH LARGEST ERROR ESTIMATE
C                       ERMAX = ELIST(MAXERR)
C
C              ELIST  - DOUBLE PRECISION
C                       VECTOR OF DIMENSION LAST CONTAINING
C                       THE ERROR ESTIMATES
C
C              IORD   - INTEGER
C                       VECTOR OF DIMENSION LAST, THE FIRST K
C                       ELEMENTS OF WHICH CONTAIN POINTERS
C                       TO THE ERROR ESTIMATES, SUCH THAT
C                       ELIST(IORD(1)),... , ELIST(IORD(K))
C                       FORM A DECREASING SEQUENCE, WITH
C                       K = LAST IF LAST.LE.(LIMIT/2+2), AND
C                       K = LIMIT+1-LAST OTHERWISE
C
C              NRMAX  - INTEGER
C                       MAXERR = IORD(NRMAX)
C
C 4.        NO SUBROUTINES OR FUNCTIONS NEEDED
C
C     ..................................................................
C
      DOUBLE PRECISION ELIST,ERMAX,ERRMAX,ERRMIN
      INTEGER I,IBEG,IDO,IORD,ISUCC,J,JBND,JUPBN,K,LAST,LIMIT,MAXERR,
     *  NRMAX
      DIMENSION ELIST(LAST),IORD(LAST)
C
C           CHECK WHETHER THE LIST CONTAINS MORE THAN
C           TWO ERROR ESTIMATES.
C
C***FIRST EXECUTABLE STATEMENT DQPSRT
      IF(LAST.GT.2) GO TO 10
      IORD(1) = 1
      IORD(2) = 2
      GO TO 90
C
C           THIS PART OF THE ROUTINE IS ONLY EXECUTED
C           IF, DUE TO A DIFFICULT INTEGRAND, SUBDIVISION
C           INCREASED THE ERROR ESTIMATE. IN THE NORMAL CASE
C           THE INSERT PROCEDURE SHOULD START AFTER THE
C           NRMAX-TH LARGEST ERROR ESTIMATE.
C
   10 ERRMAX = ELIST(MAXERR)
      IF(NRMAX.EQ.1) GO TO 30
      IDO = NRMAX-1
      DO 20 I = 1,IDO
        ISUCC = IORD(NRMAX-1)
C ***JUMP OUT OF DO-LOOP
        IF(ERRMAX.LE.ELIST(ISUCC)) GO TO 30
        IORD(NRMAX) = ISUCC
        NRMAX = NRMAX-1
   20    CONTINUE
C
C           COMPUTE THE NUMBER OF ELEMENTS IN THE LIST TO
C           BE MAINTAINED IN DESCENDING ORDER. THIS NUMBER
C           DEPENDS ON THE NUMBER OF SUBDIVISIONS STILL
C           ALLOWED.
C
   30 JUPBN = LAST
      IF(LAST.GT.(LIMIT/2+2)) JUPBN = LIMIT+3-LAST
      ERRMIN = ELIST(LAST)
C
C           INSERT ERRMAX BY TRAVERSING THE LIST TOP-DOWN,
C           STARTING COMPARISON FROM THE ELEMENT ELIST(IORD(NRMAX+1)).
C
      JBND = JUPBN-1
      IBEG = NRMAX+1
      IF(IBEG.GT.JBND) GO TO 50
      DO 40 I=IBEG,JBND
        ISUCC = IORD(I)
C ***JUMP OUT OF DO-LOOP
        IF(ERRMAX.GE.ELIST(ISUCC)) GO TO 60
        IORD(I-1) = ISUCC
   40 CONTINUE
   50 IORD(JBND) = MAXERR
      IORD(JUPBN) = LAST
      GO TO 90
C
C           INSERT ERRMIN BY TRAVERSING THE LIST BOTTOM-UP.
C
   60 IORD(I-1) = MAXERR
      K = JBND
      DO 70 J=I,JBND
        ISUCC = IORD(K)
C ***JUMP OUT OF DO-LOOP
        IF(ERRMIN.LT.ELIST(ISUCC)) GO TO 80
        IORD(K+1) = ISUCC
        K = K-1
   70 CONTINUE
      IORD(I) = LAST
      GO TO 90
   80 IORD(K+1) = LAST
C
C           SET MAXERR AND ERMAX.
C
   90 MAXERR = IORD(NRMAX)
      ERMAX = ELIST(MAXERR)
      RETURN
      END
