      SUBROUTINE ISUBX(PP,QQ,XX,ISUBXX,ERROR,EPSLON)
C  ISUBX FROM NWL REPORT 1949-DIDONATO AND JARNAGIN
      REAL ISUBXX      ,LNDIF,I,N,KHALF,JHALF,NMINUS,LOGAM,IXPHAF,LGP,LL
      REAL K,J,LNEPS,LAMBDA,LNEPP,ISXHH
      INTEGER ERROR
      DIMENSION W(20),Y(10)
      EQUIVALENCE (Y(1),W(11))
      DATA W(1) /6.6671344308688E-2/,  W(2) /1.4945134915058E-1/,
     1     W(3) /2.1908636251598E-1/,  W(4) /2.6926671931000E-1/,
     2     W(5) /2.9552422471475E-1/,  W(6) /2.9552422471475E-1/,
     3     W(7) /2.6926671931000E-1/,  W(8) /2.1908636251598E-1/,
     4     W(9) /1.4945134915058E-1/,  W(10)/6.6671344308688E-2/,
     5     W(11)/1.3046735791414E-2/,  W(12)/6.7468316655507E-2/,
     6     W(13)/1.6029521585049E-1/,  W(14)/2.8330230293538E-1/,
     7     W(15)/4.2556283050918E-1/,  W(16)/5.7443716949081E-1/,
     8     W(17)/7.1669769706462E-1/,  W(18)/8.3970478414951E-1/,
     9     W(19)/9.3253168334449E-1/,  W(20)/9.8695326420859E-1/
      DATA GHALF/.57236494292470/,PDTWO/1.5707963267949/
C     -------------------
      P = PP
      Q = QQ
      X = XX
      IF(X.EQ.0.0.OR.X.EQ.1.0)GO TO 44
      ISUBXX=0.0
      IF(P.LT.0.0.OR.Q.LT.0.0.OR.X.LT.0.0)GO TO 43
      IF(Q.GT.70.0.OR.Q.LT.0.5.OR.P.GT.1.E8.OR.P.LT.0.5.OR.X.GT.1.0)
     *     GO TO 42
      K=INT(P)
      J=INT(Q)
      JHALF=Q-J
      KHALF=P-K
      IF((JHALF.NE.0.0.AND.JHALF.NE.0.5).OR.(KHALF.NE.0.0.AND.
     *     KHALF.NE.0.5)) GO TO 41
C P,Q,AND X ARE VALID INPUT ARGUMENTS.
      ONEMX=1.0-X
      RONEMX=SQRT(ONEMX)
      EPP=EPSLON/2.0
      EPRIME=EPP/2.0
      LNEPP=ALOG(EPP)
      LNEPS=ALOG(EPSLON)
      PMINUS=P-1.
      IGO=0
      IF(JHALF.NE.0.0)GO TO 12
      IF(KHALF.NE.0.0)GO TO 2
      IF(P.GT.Q)GO TO 2
C INTERCHANGE P AND Q
    1 TEMP=Q
      Q=P
      P=TEMP
      TEMP=ONEMX
      ONEMX=X
      X=TEMP
      TEMP=JHALF
      JHALF=KHALF
      KHALF=TEMP
      TEMP=J
      J=K
      K=TEMP
      PMINUS=P-1.
      IGO=1
C BEGIN EQUATION 48 COMPUTATIONS.
    2 IF(P.EQ.0.5)GO TO 3
      N=INT( ONEMX*PMINUS/X+1.)
      IF(N.GT.J)N=J
      I=N-1.
      GO TO 4
    3 I=0.
      N=1.
    4 NMINUS=I
      ASUBN=P*ALOG(X)+LOGAM(N+PMINUS)+NMINUS*ALOG(ONEMX)-LOGAM(N)
     *     -LOGAM(P)
      IF(LNEPS-ALOG(J).LT.ASUBN)GO TO 5
      GO TO 11
    5 CONST1=1.0-EPSLON
      EXPASN=EXP(ASUBN)
      IF(EXPASN.LT.CONST1)GO TO 6
      ISUBXX=1.0
      GO TO 11
    6 LAMBDA=EXPASN
      EPSDJ=EPSLON/J
      SIG1=LAMBDA
      SIG2=0.0
    7 I=I+1.
      IF(I.GE.J)GO TO 8
      LAMBDA=((PMINUS+I)/I)*ONEMX*LAMBDA
      SIG2=SIG2+LAMBDA
      IF(EPSDJ.LT.LAMBDA)GO TO 7
    8 I=N
      LAMBDA=SIG1
    9 I=I-1.
      IF(I.EQ.0.0)GO TO 10
      LAMBDA=I*LAMBDA/((I+PMINUS)*ONEMX)
      SIG2=SIG2+LAMBDA
      IF(EPSDJ.LT.LAMBDA)GO TO 9
   10 ISUBXX=SIG1+SIG2
   11 IF(IGO.EQ.0)GO TO 45
      ISUBXX=1.0-ISUBXX
      IF (ISUBXX .LT. 0.0) ISUBXX = 0.0
      GO TO 45
   12 IF(P.GT.70.)GO TO 13
      IF(KHALF.EQ.0.0)GO TO 1
      GO TO 35
   13 TEMP3=1.0-RONEMX
      TEMP=ALOG(X)*PMINUS
      TEMP1=TEMP3*1.34
      TEMP2=SQRT(K)
      LGP=LOGAM(P)
      IF(ALOG(TEMP1*TEMP2)+TEMP.GT.LNEPP)GO TO 15
   14 IXPHAF=0.0
      GO TO 21
   15 IF(P.GT.100.)GO TO 16
      XTEMP=LOGAM(P+.5)
      GAMRAT=EXP(XTEMP-LGP-GHALF)
      GO TO 17
   16 XTEMP=LNDIF(P,.5)
      GAMRAT=EXP(XTEMP-GHALF)
   17 TEMP2 = SQRT(PMINUS)*EPP*0.56418958354776/GAMRAT
      TEMP3 = 1.0/PMINUS
      TEMP1=TEMP2**TEMP3
      TEMP  =1.0-TEMP1
      IF(TEMP  .LE.ONEMX)GO TO 14
      LAMBDA=SQRT(TEMP)
      IF(4.0*ONEMX.GT.TEMP)GO TO 19
C BEGIN EQUATION  96 COMPUTATIONS.
      TEMP1=0.0
      DO 18 L=1,10
   18 TEMP1=TEMP1+(-Y(L)*Y(L)*ONEMX+1.)**PMINUS*W(L)
      IXPHAF=TEMP1*GAMRAT*(-RONEMX)+1.
      GO TO 21
   19 TEMP1=LAMBDA-RONEMX
C BEGIN EQUATION  95 COMPUTATIONS.
      TEMP=0.0
      DO 20 L=1,10
   20 TEMP=(-( Y(L)*TEMP1 + RONEMX)**2 +1.)**PMINUS*W(L)+TEMP
      IXPHAF =TEMP*GAMRAT*TEMP1+EPRIME
   21 IF(J.EQ.0.0)GO TO 32
C EQUATION 59 COMPUTATIONS.
      N=INT(PMINUS*ONEMX/X+.5)
      IF(N.GE.J)GO TO 23
      IF(N.LT.2.)GO TO 22
      I=N-1.
      GO TO 24
   22 I=0.
      N=1.
      GO TO 24
   23 N=J
      I=J-1.
   24 NMINUS=N-1.
      TEMP=N-.5
      IF(P.GT.100.)GO TO 25
      TEMP=LOGAM(P+TEMP)-LGP
      GO TO 26
   25 TEMP=LNDIF(P,TEMP)
   26 TEMP= TEMP-LOGAM(N+.5)+P*ALOG(X)+ALOG(RONEMX)+NMINUS*ALOG(ONEMX)
      TEMP2=-ALOG(J)+LNEPP
      CONST1=1.0-LNEPP
      IF(TEMP.LE.TEMP2)GO TO 32
      SIG1=EXP(TEMP)
      LAMBDA=SIG1
      IF(SIG1+IXPHAF .LT.CONST1)GO TO 27
      ISUBXX=1.0
      GO TO 45
   27 SIG2=0.0
   28 I=I+1.
      IF(I.LT.J)GO TO 34
   29 LAMBDA=SIG1
      I=N-1.
   30 IF(I.GT.0.)GO TO 33
   31 ISUBXX=SIG1+SIG2
   32 ISUBXX=ISUBXX+IXPHAF
      GO TO 45
   33 TEMP=I+P-.5
      LAMBDA=(I+.5)/(ONEMX*TEMP)*LAMBDA
      SIG2=LAMBDA+SIG2
      IF(EPP/J.GE.LAMBDA)GO TO 31
      I=I-1.
      GO TO 30
   34 TEMP=I+.5
      LAMBDA=(TEMP-1.0+P)*ONEMX*LAMBDA/TEMP
      SIG2=SIG2+LAMBDA
      IF(LAMBDA*J.GT.EPP)GO TO 28
      GO TO 29
   35 TEMP= SQRT(X)
      ISXHH=ATAN(    TEMP /RONEMX)/PDTWO
      IF(K.EQ.0.0)GO TO 39
      TEMP=-TEMP
      LL=K
      IF(J.NE.0.0)LGP=LOGAM(P)
   36 I=0.0
      TEMP2=1.0
      IXPHAF=0.0
   37 I=I+1.
      IF(I.EQ.LL)GO TO 38
      TWOI=I+I
      TEMP2=X*(TWOI/(TWOI+1.))*TEMP2
      IXPHAF=TEMP2+IXPHAF
      GO TO 37
   38 IXPHAF=(IXPHAF+1.)*TEMP*RONEMX/PDTWO+ISXHH
      IF(K.NE.0.0)GO TO 21
      ISUBXX=IXPHAF
      GO TO 45
   39 IF(J.NE.0.0)GO TO 40
      ISUBXX=ISXHH
      GO TO 45
   40 LL=J
      TEMP2=X
      X=ONEMX
      ONEMX=TEMP2
      GO TO 36
   41 ERROR=4
      GO TO 46
   42 ERROR=3
      GO TO 46
   43 ERROR=2
      GO TO 46
   44 ISUBXX=X
   45 ERROR=1
   46 RETURN
      END
