      SUBROUTINE QCKSRT(N,ARR)
      PARAMETER (M=7,NSTACK=50,FM=7875.,FA=211.,FC=1663.
     *    ,FMI=1.2698413E-4)
      DIMENSION ARR(N),ISTACK(NSTACK)
      JSTACK=0
      L=1
      IR=N
      FX=0.
10    IF(IR-L.LT.M)THEN
        DO 13 J=L+1,IR
          A=ARR(J)
          DO 11 I=J-1,1,-1
            IF(ARR(I).LE.A)GO TO 12
            ARR(I+1)=ARR(I)
11        CONTINUE
          I=0
12        ARR(I+1)=A
13      CONTINUE
        IF(JSTACK.EQ.0)RETURN
        IR=ISTACK(JSTACK)
        L=ISTACK(JSTACK-1)
        JSTACK=JSTACK-2
      ELSE
        I=L
        J=IR
        FX=MOD(FX*FA+FC,FM)
        IQ=L+(IR-L+1)*(FX*FMI)
        A=ARR(IQ)
        ARR(IQ)=ARR(L)
20      CONTINUE
21        IF(J.GT.0)THEN
            IF(A.LT.ARR(J))THEN
              J=J-1
              GO TO 21
            ENDIF
          ENDIF
          IF(J.LE.I)THEN
            ARR(I)=A
            GO TO 30
          ENDIF
          ARR(I)=ARR(J)
          I=I+1
22        IF(I.LE.N)THEN
            IF(A.GT.ARR(I))THEN
              I=I+1
              GO TO 22
            ENDIF
          ENDIF
          IF(J.LE.I)THEN
            ARR(J)=A
            I=J
            GO TO 30
          ENDIF
          ARR(J)=ARR(I)
          J=J-1
        GO TO 20
30      JSTACK=JSTACK+2
        IF(JSTACK.GT.NSTACK)PAUSE 'NSTACK must be made larger.'
        IF(IR-I.GE.I-L)THEN
          ISTACK(JSTACK)=IR
          ISTACK(JSTACK-1)=I+1
          IR=I-1
        ELSE
          ISTACK(JSTACK)=I-1
          ISTACK(JSTACK-1)=L
          L=I+1
        ENDIF
      ENDIF
      GO TO 10
      END
      SUBROUTINE QGAUS(FUNC,A,B,SS)
      DIMENSION X(5),W(5)
      DATA X/.1488743389,.4333953941,.6794095682,.8650633666,.9739065285
     */
      DATA W/.2955242247,.2692667193,.2190863625,.1494513491,.0666713443
     */
      XM=0.5*(B+A)
      XR=0.5*(B-A)
      SS=0
      DO 11 J=1,5
        DX=XR*X(J)
        SS=SS+W(J)*(FUNC(XM+DX)+FUNC(XM-DX))
11    CONTINUE
      SS=XR*SS
      RETURN
      END
      SUBROUTINE QROMB(FUNC,A,B,SS)
      PARAMETER(EPS=1.E-6,JMAX=20,JMAXP=JMAX+1,K=5,KM=4)
      DIMENSION S(JMAXP),H(JMAXP)
      H(1)=1.
      DO 11 J=1,JMAX
        CALL TRAPZD(FUNC,A,B,S(J),J)
        IF (J.GE.K) THEN
          L=J-KM
          CALL POLINT(H(L),S(L),K,0.,SS,DSS)
          IF (ABS(DSS).LT.EPS*ABS(SS)) RETURN
        ENDIF
        S(J+1)=S(J)
        H(J+1)=0.25*H(J)
11    CONTINUE
      PAUSE 'Too many steps.'
      END
      SUBROUTINE QROMO(FUNC,A,B,SS,CHOOSE)
      PARAMETER (EPS=1.E-6,JMAX=14,JMAXP=JMAX+1,KM=4,K=KM+1)
      DIMENSION S(JMAXP),H(JMAXP)
      H(1)=1.
      DO 11 J=1,JMAX
        CALL CHOOSE(FUNC,A,B,S(J),J)
        IF (J.GE.K) THEN
          CALL POLINT(H(J-KM),S(J-KM),K,0.0,SS,DSS)
          IF (ABS(DSS).LT.EPS*ABS(SS)) RETURN
        ENDIF
        S(J+1)=S(J)
        H(J+1)=H(J)/9.
11    CONTINUE
      PAUSE 'Too many steps.'
      END
      SUBROUTINE QROOT(P,N,B,C,EPS)
      PARAMETER (NMAX=20,ITMAX=20,TINY=1.0E-6)
      DIMENSION P(N),Q(NMAX),D(3),REM(NMAX),QQ(NMAX)
      D(3)=1.
      DO 12 ITER=1,ITMAX
        D(2)=B
        D(1)=C
        CALL POLDIV(P,N,D,3,Q,REM)
        S=REM(1)
        R=REM(2)
        CALL POLDIV(Q,N-1,D,3,QQ,REM)
        SC=-REM(1)
        RC=-REM(2)
        DO 11 I=N-1,1,-1
          Q(I+1)=Q(I)
11      CONTINUE
        Q(1)=0.
        CALL POLDIV(Q,N,D,3,QQ,REM)
        SB=-REM(1)
        RB=-REM(2)
        DIV=1./(SB*RC-SC*RB)
        DELB=(R*SC-S*RC)*DIV
        DELC=(-R*SB+S*RB)*DIV
        B=B+DELB
        C=C+DELC
        IF((ABS(DELB).LE.EPS*ABS(B).OR.ABS(B).LT.TINY)
     *      .AND.(ABS(DELC).LE.EPS*ABS(C)
     *      .OR.ABS(C).LT.TINY)) RETURN
12    CONTINUE
      PAUSE 'too many iterations in QROOT'
      END
      SUBROUTINE QSIMP(FUNC,A,B,S)
      PARAMETER (EPS=1.E-6, JMAX=20)
      OST=-1.E30
      OS= -1.E30
      DO 11 J=1,JMAX
        CALL TRAPZD(FUNC,A,B,ST,J)
        S=(4.*ST-OST)/3.
        IF (ABS(S-OS).LT.EPS*ABS(OS)) RETURN
        OS=S
        OST=ST
11    CONTINUE
      PAUSE 'Too many steps.'
      END
      SUBROUTINE QTRAP(FUNC,A,B,S)
      PARAMETER (EPS=1.E-6, JMAX=20)
      OLDS=-1.E30
      DO 11 J=1,JMAX
        CALL TRAPZD(FUNC,A,B,S,J)
        IF (ABS(S-OLDS).LT.EPS*ABS(OLDS)) RETURN
        OLDS=S
11    CONTINUE
      PAUSE 'Too many steps.'
      END
      SUBROUTINE QUAD3D(X1,X2,SS)
      EXTERNAL H
      CALL QGAUSX(H,X1,X2,SS)
      RETURN
      END
C     
      FUNCTION F(ZZ)
      EXTERNAL FUNC
      COMMON /XYZ/ X,Y,Z
      Z=ZZ
      F=FUNC(X,Y,Z)
      RETURN
      END
C     
      FUNCTION G(YY)
      EXTERNAL Z1,Z2,F
      COMMON /XYZ/ X,Y,Z
      Y=YY
      CALL QGAUSZ(F,Z1(X,Y),Z2(X,Y),SS)
      G=SS
      RETURN
      END
C     
      FUNCTION H(XX)
      EXTERNAL Y1,Y2,G
      COMMON /XYZ/ X,Y,Z
      X=XX
      CALL QGAUSY(G,Y1(X),Y2(X),SS)
      H=SS
      RETURN
      END
