FUNCTION GAMLN(X)                                                 GAM   10
C
C         A CDC 6600 SUBROUTINE
C
C     AUTHORS
C         D.E. AMOS AND S.L. DANIEL
C         ALBUQUERQUE, NEW MEXICO, 87115
C         JANUARY, 1975
C
C     REFERENCES
C         ABRAMOWITZ, M. AND STEGUN, I.A. HANDBOOK OF MATHEMATICAL
C         FUNCTIONS. NBS APPLIED MATHEMATICS SERIES 55, U.S. GOVERNMENT
C         PRINTING OFFICE, WASHINGTON, D.C., CHAPTER 6.
C
C         AMOS, D.E., DANIEL, S.L. AND WESTON, M.K. CDC 6600
C         SUBROUTINES IBESS AND JBESS FOR BESSEL FUNCTIONS
C         I/SUB(NU)/(X) AND J/SUB(NU)/(X), X.GE.0, NU.GE.0.
C         ACM TRANS. MATH. SOFTWARE, MARCH, 1977.
C
C         HART, J.F., ET. AL. COMPUTER APPROXIMATIONS, WILEY, NEW YORK.
C         PP. 130-136, 1968.
C
C     ABSTRACT
C         GAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR
C         X.GT.0. A RATIONAL CHEBYSHEV APPROXIMATION IS USED ON
C         8.LT.X.LT.1000., THE ASYMPTOTIC EXPANSION FOR X.GE.1000. AND
C         BACKWARD RECURSION FOR 0.LT.X.LT.8 FOR NON-INTEGRAL X. FOR
C         X=1.,...,8., GAMLN IS SET TO NATURAL LOGS OF FACTORIALS.
C
C     DESCRIPTION OF ARGUMENTS
C
C         INPUT
C           X      - X.GT.0
C
C         OUTPUT
C           GAMLN  - NATURAL LOG OF THE GAMMA FUNCTION AT X
C
C     ERROR CONDITIONS
C         IMPROPER INPUT ARGUMENT - A FATAL ERROR
      DOUBLE PRECISION DXZ, RXZ, XK
      DIMENSION G(8), P(5), Q(2)
      DATA XLIM1,XLIM2,RTWPIL/    8.  ,  1000.   , 9.18938533204673E-01/
      DATA  G              / 8.52516136106541E+00, 6.57925121201010E+00,
     1 4.78749174278205E+00, 3.17805383034795E+00, 1.79175946922806E+00,
     2 6.93147180559945E-01,         2(0.)       /
      DATA  P              / 7.66345188000000E-04,-5.94095610520000E-04,
     1 7.93643110484500E-04,-2.77777775657725E-03, 8.33333333333169E-02/
      DATA  Q              /-2.77777777777778E-03, 8.33333333333333E-02/
      IF (X) 140, 140, 10
   10 DX = X - XLIM1
      IF (DX) 20, 110, 80
   20 IF (X-1.) 30, 130, 40
   30 XZ = X + 8.
      TX = X
      FK = -.5
      NDX = 7
      GO TO 50
   40 DX = ABS(DX)
      NDX = DX
      DNDX = NDX
      NDX = NDX + 1
      IF ((DNDX-DX).EQ.0.) GO TO 120
      XZ = X + DNDX + 1.
      TX = 1.
      FK = .5
   50 DXZ = XZ
      RXZ = 1.D+0/DXZ
      RX = RXZ
      RXX = RX*RX
      XK = 1.D+0
      DO 60 I=1,NDX
        XK = XK - RXZ
        SXK = XK
        TX = TX*SXK
   60 CONTINUE
      SUM = (X-FK)*ALOG(XZ) - ALOG(TX) - XZ
      PX = P(1)
      DO 70 I=2,5
        PX = PX*RXX + P(I)
   70 CONTINUE
      GAMLN = PX*RX + SUM + RTWPIL
      RETURN
   80 RX = 1./X
      RXX = RX*RX
      IF ((X-XLIM2).LT.0.) GO TO 90
      PX = Q(1)*RXX + Q(2)
      GAMLN = PX*RX + (X-.5)*ALOG(X) - X + RTWPIL
      RETURN
   90 PX = P(1)
      SUM = (X-.5)*ALOG(X) - X
      DO 100 I=2,5
        PX = PX*RXX + P(I)
  100 CONTINUE
      GAMLN = PX*RX + SUM + RTWPIL
      RETURN
  110 GAMLN = G(1)
      RETURN
  120 GAMLN = G(NDX)
      RETURN
  130 GAMLN = G(8)
      RETURN
  140 PRINT 99999, X
      STOP
99999 FORMAT (50H ARGUMENT FOR GAMLN IS LESS THAN OR EQUAL TO ZERO,,
     * 3H X=,E25.14)
      END
      SUBROUTINE IBESS(KODE, ALPHA, N, X, Y)                            IBE   10
C
C         A CDC 6600 SUBROUTINE
C
C     AUTHORS
C         D.E. AMOS AND S.L. DANIEL
C         SANDIA LABORATORIES
C         JANUARY, 1975
C
C     REFERENCES
C         ABRAMOWITZ, M. AND STEGUN, I.A. HANDBOOK OF MATHEMATICAL
C         FUNCTIONS. NBS APPLIED MATHEMATICS SERIES 55, U.S. GOVERNMENT
C         PRINTING OFFICE, WASHINGTON, D.C., CHAPTERS 9 AND 10.
C
C         AMOS, D.E., DANIEL, S.L. AND WESTON, M.K. CDC 6600
C         SUBROUTINES IBESS AND JBESS FOR BESSEL FUNCTIONS
C         I/SUB(NU)/(X) AND J/SUB(NU)/(X), X.GE.0, NU.GE.0.
C         ACM TRANS. MATH. SOFTWARE, MARCH, 1977.
C
C         OLVER, F.W.J. TABLES FOR BESSEL FUNCTIONS OF MODERATE OR
C         LARGE ORDERS. NPL MATHEMATICAL TABLES, VOL 6. HER MAJESTY-S
C         STATIONERY OFFICE, LONDON, 1962.
C
C         OLVER, F.W.J. THE ASYMPTOTIC EXPANSION OF BESSEL FUNCTIONS OF
C         LARGE ORDER. PHIL. TRANS. A, 247, PP. 328-368, 1954.
C
C     ABSTRACT
C         IBESS COMPUTES AN N MEMBER SEQUENCE OF I BESSEL FUNCTIONS
C         I/SUB(ALPHA+K-1)/(X), K=1,...,N OR SCALED BESSEL FUNCTIONS
C         EXP(-X)*I/SUB(ALPHA*K-1)/(X), K=1,...,N FOR NON-NEGATIVE ALPHA
C         AND X. A COMBINATION OF THE POWER SERIES, THE ASYMPTOTIC
C         EXPANSION FOR X TO INFINITY, AND THE UNIFORM ASYMPTOTIC
C         EXPANSION FOR NU TO INFINITY ARE APPLIED OVER SUBDIVISIONS OF
C         THE (NU,X) PLANE. FOR VALUES NOT COVERED BY ONE OF THESE
C         FORMULAE, THE ORDER IS INCREMENTED BY AN INTEGER SO THAT ONE
C         OF THESE FORMULAE APPLY. BACKWARD RECURSION IS USED TO REDUCE
C         ORDERS BY INTEGER VALUES. THE ASYMPTOTIC EXPANSION FOR X TO
C         INFINITY IS USED ONLY WHEN THE ENTIRE SEQUENCE (SPECIFICALLY
C         THE LAST MEMBER) LIES WITHIN THE REGION COVERED BY THE
C         EXPANSION. LEADING TERMS OF THESE EXPANSIONS ARE USED TO TEST
C         FOR OVER OR UNDERFLOW WHERE APPROPRIATE. IF A SEQUENCE IS
C         REQUESTED AND THE LAST MEMBER WOULD UNDERFLOW, THE RESULT IS
C         SET TO ZERO AND THE NEXT LOWER ORDER TRIED, ETC., UNTIL A
C         MEMBER COMES ON SCALE OR ALL ARE SET TO ZERO. AN OVERFLOW
C         CANNOT OCCUR WITH SCALING. IBESS CALLS FUNCTION GAMLN.
C
C     DESCRIPTION OF ARGUMENTS
C
C         INPUT
C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
C                    KODE=1 RETURNS
C                           Y(K)=        I/SUB(ALPHA+K-1)/(X),
C                                K=1,...,N
C                    KODE=2 RETURNS
C                           Y(K)=EXP(-X)*I/SUB(ALPHA+K-1)/(X),
C                                K=1,...,N
C           ALPHA  - ORDER OF FIRST MEMBER OF THE SEQUENCE, ALPHA.GE.0
C           N      - NUMBER OF MEMBERS IN THE SEQUENCE, N.GE.1
C           X      - X.GE.0
C
C         OUTPUT
C           Y      - A VECTOR WHOSE FIRST N COMPONENTS CONTAIN
C                    VALUES FOR I/SUB(ALPHA+K-1)/(X) OR SCALED
C                    VALUES FOR EXP(-X)*I/SUB(ALPHA+K-1)/(X),
C                    K=1,...,N DEPENDING ON KODE
C
C     ERROR CONDITIONS
C         IMPROPER INPUT ARGUMENTS - A FATAL ERROR
C         OVERFLOW WITH KODE=1 - A FATAL ERROR
C         UNDERFLOW - A NON-FATAL ERROR
C
      DOUBLE PRECISION DX, TRX, DTM, DFN
      DIMENSION Y(1), TEMP(3)
      DIMENSION C(11,10)
      DIMENSION C1(11,8), C2(11,2)
      EQUIVALENCE (C(1,1),C1(1,1))
      EQUIVALENCE (C(1,9),C2(1,1))
      DATA ELIM,TOL        /        667.       ,       1.E-15          /
      DATA RTPI,RTTPI      / 1.59154943091895E-01, 3.98942280401433E-01/
      DATA CE              / 3.45387763900000E+01/
      DATA INLIM           /          80         /
      DATA C1              /-2.08333333333333E-01, 1.25000000000000E-01,
     1         9(0.)       , 3.34201388888889E-01,-4.01041666666667E-01,
     2 7.03125000000000E-02,         8(0.)       ,-1.02581259645062E+00,
     3 1.84646267361111E+00,-8.91210937500000E-01, 7.32421875000000E-02,
     4         7(0.)       , 4.66958442342625E+00,-1.12070026162230E+01,
     5 8.78912353515625E+00,-2.36408691406250E+00, 1.12152099609375E-01,
     6         6(0.)       ,-2.82120725582002E+01, 8.46362176746007E+01,
     7-9.18182415432400E+01, 4.25349987453885E+01,-7.36879435947963E+00,
     8 2.27108001708984E-01,         5(0.)       , 2.12570130039217E+02,
     9-7.65252468141182E+02, 1.05999045252800E+03,-6.99579627376133E+02,
     A 2.18190511744212E+02,-2.64914304869516E+01, 5.72501420974731E-01,
     B         4(0.)       ,-1.91945766231841E+03, 8.06172218173731E+03,
     C-1.35865500064341E+04, 1.16553933368645E+04,-5.30564697861340E+03,
     D 1.20090291321635E+03,-1.08090919788395E+02, 1.72772750258446E+00,
     E         3(0.)       , 2.02042913309661E+04,-9.69805983886375E+04,
     F 1.92547001232532E+05,-2.03400177280416E+05, 1.22200464983017E+05,
     G-4.11926549688976E+04, 7.10951430248936E+03,-4.93915304773088E+02,
     H 6.07404200127348E+00,         2(0.)       /
      DATA C2              /-2.42919187900551E+05, 1.31176361466298E+06,
     1-2.99801591853811E+06, 3.76327129765640E+06,-2.81356322658653E+06,
     2 1.26836527332162E+06,-3.31645172484564E+05, 4.52187689813627E+04,
     3-2.49983048181121E+03, 2.43805296995561E+01,         1(0.)       ,
     4 3.28446985307204E+06,-1.97068191184322E+07, 5.09526024926646E+07,
     5-7.41051482115327E+07, 6.63445122747290E+07,-3.75671766607634E+07,
     6 1.32887671664218E+07,-2.78561812808645E+06, 3.08186404612662E+05,
     7-1.38860897537170E+04, 1.10017140269247E+02/
C TEST INPUT ARGUMENTS
      KT = 1
C TEST INPUT ARGUMENTS
      IF (N-1) 580, 10, 20
   10 KT = 2
   20 NN = N
      IF (KODE.LT.1 .OR. KODE.GT.2) GO TO 560
      IF (X) 590, 30, 80
   30 IF (ALPHA) 570, 40, 50
   40 Y(1) = 1.
      IF (N.EQ.1) RETURN
      I1 = 2
      GO TO 60
   50 I1 = 1
   60 DO 70 I=I1,N
        Y(I) = 0.
   70 CONTINUE
      RETURN
   80 CONTINUE
      IF (ALPHA.LT.0.) GO TO 570
      DFN = DBLE(FLOAT(N)) + DBLE(ALPHA) - 1.D+0
      FNU = DFN
      IN = 0
      XO2 = X*.5
      SXO2 = XO2*XO2
      ETX = KODE - 1
      SX = ETX*X
C DECISION TREE FOR REGION WHERE SERIES, ASYMPTOTIC EXPANSION FOR X
C TO INFINITY AND ASYMPTOTIC EXPANSION FOR NU TO INFINITY ARE
C APPLIED.
      IF (SXO2.LE.(FNU+1.)) GO TO 90
      IF (X.LE.12.) GO TO 110
      FN = 0.55*FNU*FNU
      FN = AMAX1(17.,FN)
      IF (X.GE.FN) GO TO 430
      NS = AMAX1(36.-FNU,0.)
      DFN = DFN + DBLE(FLOAT(NS))
      FN = DFN
      IS = KT
      KM = N - 1 + NS
      IF (KM.GT.0) IS = 3
      GO TO 120
   90 FN = FNU
      FNP1 = FN + 1.
      XO2L = ALOG(XO2)
      IS = KT
      IF (X.LE.0.5) GO TO 230
      NS = 0
  100 DFN = DFN + DBLE(FLOAT(NS))
      FN = DFN
      FNP1 = FN + 1.
      IS = KT
      IF (N-1+NS.GT.0) IS = 3
      GO TO 230
  110 XO2L = ALOG(XO2)
      NS = SXO2 - FNU
      GO TO 100
  120 CONTINUE
C OVERFLOW TEST ON UNIFORM ASYMPTOTIC EXPANSION
      IF (KODE.EQ.2) GO TO 130
      IF (ALPHA.LT.1.) GO TO 150
      Z = X/ALPHA
      RA = SQRT(1.+Z*Z)
      GLN = ALOG((1.+RA)/Z)
      T = RA*(1.-ETX) + ETX/(Z+RA)
      ARG = ALPHA*(T-GLN)
      IF (ARG.GT.ELIM) GO TO 600
      IF (KM.EQ.0) GO TO 140
  130 CONTINUE
C UNDERFLOW TEST ON UNIFORM ASYMPTOTIC EXPANSION
      Z = X/FN
      RA = SQRT(1.+Z*Z)
      GLN = ALOG((1.+RA)/Z)
      T = RA*(1.-ETX) + ETX/(Z+RA)
      ARG = FN*(T-GLN)
  140 IF (ARG.LT.-ELIM) GO TO 280
      GO TO 190
  150 IF (X.GT.ELIM) GO TO 600
      GO TO 130
C UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY
  160 IF (KM.NE.0) GO TO 170
      Y(1) = TEMP(3)
      RETURN
  170 TEMP(1) = TEMP(3)
      IN = NS
      KT = 1
  180 CONTINUE
      IS = 2
      DFN = DFN - 1.D+0
      FN = DFN
      Z = X/FN
      RA = SQRT(1.+Z*Z)
      GLN = ALOG((1.+RA)/Z)
      T = RA*(1.-ETX) + ETX/(Z+RA)
      ARG = FN*(T-GLN)
  190 COEF = EXP(ARG)
      T = 1./RA
      T2 = T*T
      T = T/FN
      S2 = 1.
      AP = 1.
      DO 210 K=1,10
        KP1 = K + 1
        S1 = C(1,K)
        DO 200 J=2,KP1
          S1 = S1*T2 + C(J,K)
  200   CONTINUE
        AP = AP*T
        AK = AP*S1
        S2 = S2 + AK
        IF (ABS(AK).LT.TOL) GO TO 220
  210 CONTINUE
  220 CONTINUE
      TEMP(IS) = SQRT(T*RTPI)*COEF*S2
      GO TO (180, 350, 500), IS
C SERIES FOR (X/2)**2.LE.NU+1
  230 CONTINUE
      GLN = GAMLN(FNP1)
      ARG = FN*XO2L - GLN - SX
      IF (ARG.LT.-ELIM) GO TO 300
      EARG = EXP(ARG)
  240 CONTINUE
      S = 1.
      AK = 3.
      T2 = 1.
      T = 1.
      S1 = FN
      DO 250 K=1,17
        S2 = T2 + S1
        T = T*SXO2/S2
        S = S + T
        IF (ABS(T).LT.TOL) GO TO 260
        T2 = T2 + AK
        AK = AK + 2.
        S1 = S1 + FN
  250 CONTINUE
  260 CONTINUE
      TEMP(IS) = S*EARG
      GO TO (270, 350, 490), IS
  270 EARG = EARG*FN/XO2
      DFN = DFN - 1.D+0
      FN = DFN
      IS = 2
      GO TO 240
C SET UNDERFLOW VALUE AND UPDATE PARAMETERS
  280 Y(NN) = 0.
      NN = NN - 1
      DFN = DFN - 1.D+0
      FN = DFN
      IF (NN-1) 340, 290, 130
  290 KT = 2
      IS = 2
      GO TO 130
  300 Y(NN) = 0.
      NN = NN - 1
      FNP1 = FN
      DFN = DFN - 1.D+0
      FN = DFN
      IF (NN-1) 340, 310, 320
  310 KT = 2
      IS = 2
  320 IF (SXO2.LE.FNP1) GO TO 330
      GO TO 130
  330 ARG = ARG - XO2L + ALOG(FNP1)
      IF (ARG.LT.-ELIM) GO TO 300
      GO TO 230
  340 NZ = N - NN
      PRINT 99994, NZ, KODE, ALPHA, N, X
      RETURN
C BACKWARD RECURSION SECTION
  350 CONTINUE
      NZ = N - NN
      IF (NZ.NE.0) PRINT 99994, NZ, KODE, ALPHA, N, X
  360 GO TO (370, 420), KT
  370 CONTINUE
      S1 = TEMP(1)
      S2 = TEMP(2)
      DX = X
      TRX = 2.D+0/DX
      DTM = DFN*TRX
      TM = DTM
      IF (IN.EQ.0) GO TO 390
C BACKWARD RECUR TO INDEX ALPHA+NN-1
      DO 380 I=1,IN
        S = S2
        S2 = TM*S2 + S1
        S1 = S
        DTM = DTM - TRX
        TM = DTM
  380 CONTINUE
      Y(NN) = S1
      IF (NN.EQ.1) RETURN
      Y(NN-1) = S2
      IF (NN.EQ.2) RETURN
      GO TO 400
  390 CONTINUE
C BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
      Y(NN) = S1
      Y(NN-1) = S2
      IF (NN.EQ.2) RETURN
  400 K = NN + 1
      DO 410 I=3,NN
        K = K - 1
        Y(K-2) = TM*Y(K-1) + Y(K)
        DTM = DTM - TRX
        TM = DTM
  410 CONTINUE
      RETURN
  420 Y(1) = TEMP(2)
      RETURN
C ASYMPTOTIC EXPANSION FOR X TO INFINITY
  430 CONTINUE
      EARG = RTTPI/SQRT(X)
      IF (KODE.EQ.2) GO TO 440
      IF (X.GT.ELIM) GO TO 600
      EARG = EARG*EXP(X)
  440 ETX = 8.*X
      IS = KT
      IN = 0
  450 DX = DFN + DFN
      DTM = DX*DX
      S1 = ETX
      TRX = S1
      DX = -(DTM-1.D+0)/TRX
      T = DX
      TRX = 1.D+0 + DX
      S = TRX
      S2 = 1.
      AK = 8.
      DO 460 K=1,25
        S1 = S1 + ETX
        S2 = S2 + AK
        DX = S2
        TRX = DTM - DX
        AP = TRX
        T = -T*AP/S1
        S = S + T
        IF (ABS(T).LE.TOL) GO TO 470
        AK = AK + 8.
  460 CONTINUE
  470 TEMP(IS) = S*EARG
      GO TO (480, 360), IS
  480 IS = 2
      DFN = DFN - 1.D+0
      GO TO 450
C BACKWARD RECURSION WITH NORMALIZATION BY
C ASYMPTOTIC EXPANSION FOR NU TO INFINITY OR POWER SERIES.
  490 CONTINUE
C COMPUTATION OF LAST ORDER FOR SERIES NORMALIZATION
      KM = AMAX1(3.-FN,0.)
      TFN = FN + FLOAT(KM)
      TA = (GLN+TFN-0.9189385332-0.0833333333/TFN)/(TFN+0.5)
      TA = XO2L - TA
      TB = -(1.-1./TFN)/TFN
      IN = CE/(-TA+SQRT(TA*TA-CE*TB)) + 1.5
      IN = IN + KM
      GO TO 510
  500 CONTINUE
C COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION
      IN = CE/(GLN+SQRT(GLN*GLN+T*CE)) + 1.5
      IF (IN.GT.INLIM) GO TO 160
  510 DX = FLOAT(IN)
      DTM = DFN + DX
      DX = X
      TRX = 2.D+0/DX
      DTM = DTM*TRX
      TM = DTM
      TA = 0.
      TB = TOL
      KK = 1
  520 CONTINUE
C BACKWARD RECUR UNINDEXED
      DO 530 I=1,IN
        S = TB
        TB = TM*TB + TA
        TA = S
        DTM = DTM - TRX
        TM = DTM
  530 CONTINUE
C NORMALIZATION
      IF (KK.NE.1) GO TO 540
      TA = (TA/TB)*TEMP(3)
      TB = TEMP(3)
      KK = 2
      IN = NS
      IF (NS.NE.0) GO TO 520
  540 Y(NN) = TB
      NZ = N - NN
      IF (NZ.NE.0) PRINT 99994, NZ, KODE, ALPHA, N, X
      IF (NN.EQ.1) RETURN
      TB = TM*TB + TA
      K = NN - 1
      Y(K) = TB
      IF (NN.EQ.2) RETURN
      DTM = DTM - TRX
      TM = DTM
      KM = K - 1
C BACKWARD RECUR INDEXED
      DO 550 I=1,KM
        Y(K-1) = TM*Y(K) + Y(K+1)
        DTM = DTM - TRX
        TM = DTM
        K = K - 1
  550 CONTINUE
      RETURN
  560 PRINT 99999, KODE, ALPHA, N, X
      STOP
  570 PRINT 99998, KODE, ALPHA, N, X
      STOP
  580 PRINT 99997, KODE, ALPHA, N, X
      STOP
  590 PRINT 99996, KODE, ALPHA, N, X
      STOP
  600 PRINT 99995, KODE, ALPHA, N, X
      STOP
99999 FORMAT (51H0IBESS CALLED WITH SCALING OPTION, KODE, NOT 1 OR 2
     * /6H KODE=,I2,7H ALPHA=,E25.14,3H N=,I6,3H X=,E25.14)
99998 FORMAT (51H0IBESS CALLED WITH THE ORDER, ALPHA, LESS THAN ZERO
     * /6H KODE=,I2,7H ALPHA=,E25.14,3H N=,I6,3H X=,E25.14)
99997 FORMAT (34H0IBESS CALLED WITH N LESS THAN ONE/6H KODE=,
     * I2,7H ALPHA=,E25.14,3H N=,I6,3H X=,E25.14)
99996 FORMAT (35H0IBESS CALLED WITH X LESS THAN ZERO/6H KODE=,
     * I2,7H ALPHA=,E25.14,3H N=,I6,3H X=,E25.14)
99995 FORMAT (42H0OVERFLOW IN IBESS, X TOO LARGE FOR KODE=1/6H KODE=,
     * I2,7H ALPHA=,E25.14,3H N=,I6,3H X=,E25.14)
99994 FORMAT (25H0UNDERFLOW IN IBESS, LAST,I6,20H VALUE(S) OF Y ARRAY,
     * 17H WERE SET TO ZERO/6H KODE=,I2,7H ALPHA=,E25.14,3H N=,I6,
     * 3H X=,E25.14)
      END
      SUBROUTINE JAIRY(X, RX, C, AI, DAI)                               JAI   10
C
C     CDC 6600 ROUTINE
C     1-2-74
C
C                  JAIRY COMPUTES THE AIRY FUNCTION AI(X)
C                   AND ITS DERIVATIVE DAI(X) FOR JBESS
C                                   INPUT
C
C         X - ARGUMENT, COMPUTED BY JBESS, X.LE.(ELIM2*1.5)**(2./3.)
C        RX - RX=SQRT(ABS(X)), COMPUTED BY JBESS
C         C - C=2.*(ABS(X)**1.5)/3., COMPUTED BY JBESS
C
C                                  OUTPUT
C
C        AI - VALUE OF FUNCTION AI(X)
C       DAI - VALUE OF THE DERIVATIVE DAI(X)
C
C                                WRITTEN BY
C
C                                D. E. AMOS
C                               S. L. DANIEL
C                               M. K. WESTON
C
      DIMENSION AJP(19), AJN(19), A(15), B(15)
      DIMENSION AK1(14), AK2(23), AK3(14)
      DIMENSION DAJP(19), DAJN(19), DA(15), DB(15)
      DIMENSION DAK1(14), DAK2(24), DAK3(14)
      DATA N1, N2, N3, N4 /14,23,19,15/
      DATA M1, M2, M3, M4 /12,21,17,13/
      DATA FPI12,CON2,CON3,CON4,CON5             / 1.30899693899575E+00,
     1 5.03154716196777E+00, 3.80004589867293E-01, 8.33333333333333E-01,
     2 8.66025403784439E-01/
      DATA  AK1            / 2.20423090987793E-01,-1.25290242787700E-01,
     1 1.03881163359194E-02, 8.22844152006343E-04,-2.34614345891226E-04,
     2 1.63824280172116E-05, 3.06902589573189E-07,-1.29621999359332E-07,
     3 8.22908158823668E-09, 1.53963968623298E-11,-3.39165465615682E-11,
     4 2.03253257423626E-12,-1.10679546097884E-14,-5.16169497785080E-15/
      DATA  AK2            / 2.74366150869598E-01, 5.39790969736903E-03,
     1-1.57339220621190E-03, 4.27427528248750E-04,-1.12124917399925E-04,
     2 2.88763171318904E-05,-7.36804225370554E-06, 1.87290209741024E-06,
     3-4.75892793962291E-07, 1.21130416955909E-07,-3.09245374270614E-08,
     4 7.92454705282654E-09,-2.03902447167914E-09, 5.26863056595742E-10,
     5-1.36704767639569E-10, 3.56141039013708E-11,-9.31388296548430E-12,
     6 2.44464450473635E-12,-6.43840261990955E-13, 1.70106030559349E-13,
     7-4.50760104503281E-14, 1.19774799164811E-14,-3.19077040865066E-15/
      DATA  AK3            / 2.80271447340791E-01,-1.78127042844379E-03,
     1 4.03422579628999E-05,-1.63249965269003E-06, 9.21181482476768E-08,
     2-6.52294330229155E-09, 5.47138404576546E-10,-5.24408251800260E-11,
     3 5.60477904117209E-12,-6.56375244639313E-13, 8.31285761966247E-14,
     4-1.12705134691063E-14, 1.62267976598129E-15,-2.46480324312426E-16/
      DATA  AJP            / 7.78952966437581E-02,-1.84356363456801E-01,
     1 3.01412605216174E-02, 3.05342724277608E-02,-4.95424702513079E-03,
     2-1.72749552563952E-03, 2.43137637839190E-04, 5.04564777517082E-05,
     3-6.16316582695208E-06,-9.03986745510768E-07, 9.70243778355884E-08,
     4 1.09639453305205E-08,-1.04716330588766E-09,-9.60359441344646E-11,
     5 8.25358789454134E-12, 6.36123439018768E-13,-4.96629614116015E-14,
     6-3.29810288929615E-15, 2.35798252031104E-16/
      DATA  AJN            / 3.80497887617242E-02,-2.45319541845546E-01,
     1 1.65820623702696E-01, 7.49330045818789E-02,-2.63476288106641E-02,
     2-5.92535597304981E-03, 1.44744409589804E-03, 2.18311831322215E-04,
     3-4.10662077680304E-05,-4.66874994171766E-06, 7.15218807277160E-07,
     4 6.52964770854633E-08,-8.44284027565946E-09,-6.44186158976978E-10,
     5 7.20802286505285E-11, 4.72465431717846E-12,-4.66022632547045E-13,
     6-2.67762710389189E-14, 2.36161316570019E-15/
      DATA   A             / 4.90275424742791E-01, 1.57647277946204E-03,
     1-9.66195963140306E-05, 1.35916080268815E-07, 2.98157342654859E-07,
     2-1.86824767559979E-08,-1.03685737667141E-09, 3.28660818434328E-10,
     3-2.57091410632780E-11,-2.32357655300677E-12, 9.57523279048255E-13,
     4-1.20340828049719E-13,-2.90907716770715E-15, 4.55656454580149E-15,
     5-9.99003874810259E-16/
      DATA   B             / 2.78593552803079E-01,-3.52915691882584E-03,
     1-2.31149677384994E-05, 4.71317842263560E-06,-1.12415907931333E-07,
     2-2.00100301184339E-08, 2.60948075302193E-09,-3.55098136101216E-11,
     3-3.50849978423875E-11, 5.83007187954202E-12,-2.04644828753326E-13,
     4-1.10529179476742E-13, 2.87724778038775E-14,-2.88205111009939E-15,
     5-3.32656311696166E-16/
      DATA N1D, N2D, N3D, N4D /14,24,19,15/
      DATA M1D, M2D, M3D, M4D /12,22,17,13/
      DATA  DAK1           / 2.04567842307887E-01,-6.61322739905664E-02,
     1-8.49845800989287E-03, 3.12183491556289E-03,-2.70016489829432E-04,
     2-6.35636298679387E-06, 3.02397712409509E-06,-2.18311195330088E-07,
     3-5.36194289332826E-10, 1.13098035622310E-09,-7.43023834629073E-11,
     4 4.28804170826891E-13, 2.23810925754539E-13,-1.39140135641182E-14/
      DATA  DAK2           / 2.93332343883230E-01,-8.06196784743112E-03,
     1 2.42540172333140E-03,-6.82297548850235E-04, 1.85786427751181E-04,
     2-4.97457447684059E-05, 1.32090681239497E-05,-3.49528240444943E-06,
     3 9.24362451078835E-07,-2.44732671521867E-07, 6.49307837648910E-08,
     4-1.72717621501538E-08, 4.60725763604656E-09,-1.23249055291550E-09,
     5 3.30620409488102E-10,-8.89252099772401E-11, 2.39773319878298E-11,
     6-6.48013921153450E-12, 1.75510132023731E-12,-4.76303829833637E-13,
     7 1.29498241100810E-13,-3.52679622210430E-14, 9.62005151585923E-15,
     8-2.62786914342292E-15/
      DATA  DAK3           / 2.84675828811349E-01, 2.53073072619080E-03,
     1-4.83481130337976E-05, 1.84907283946343E-06,-1.01418491178576E-07,
     2 7.05925634457153E-09,-5.85325291400382E-10, 5.56357688831339E-11,
     3-5.90889094779500E-12, 6.88574353784436E-13,-8.68588256452194E-14,
     4 1.17374762617213E-14,-1.68523146510923E-15, 2.55374773097056E-16/
      DATA  DAJP           / 6.53219131311457E-02,-1.20262933688823E-01,
     1 9.78010236263823E-03, 1.67948429230505E-02,-1.97146140182132E-03,
     2-8.45560295098867E-04, 9.42889620701976E-05, 2.25827860945475E-05,
     3-2.29067870915987E-06,-3.76343991136919E-07, 3.45663933559565E-08,
     4 4.29611332003007E-09,-3.58673691214989E-10,-3.57245881361895E-11,
     5 2.72696091066336E-12, 2.26120653095771E-13,-1.58763205238303E-14,
     6-1.12604374485125E-15, 7.31327529515367E-17/
      DATA  DAJN           / 1.08594539632967E-02, 8.53313194857091E-02,
     1-3.15277068113058E-01,-8.78420725294257E-02, 5.53251906976048E-02,
     2 9.41674060503241E-03,-3.32187026018996E-03,-4.11157343156826E-04,
     3 1.01297326891346E-04, 9.87633682208396E-06,-1.87312969812393E-06,
     4-1.50798500131468E-07, 2.32687669525394E-08, 1.59599917419225E-09,
     5-2.07665922668385E-10,-1.24103350500302E-11, 1.39631765331043E-12,
     6 7.39400971155740E-14,-7.32887475627500E-15/
      DATA   DA            / 4.91627321104601E-01, 3.11164930427489E-03,
     1 8.23140762854081E-05,-4.61769776172142E-06,-6.13158880534626E-08,
     2 2.87295804656520E-08,-1.81959715372117E-09,-1.44752826642035E-10,
     3 4.53724043420422E-11,-3.99655065847223E-12,-3.24089119830323E-13,
     4 1.62098952568741E-13,-2.40765247974057E-14, 1.69384811284491E-16,
     5 8.17900786477396E-16/
      DATA   DB            /-2.77571356944231E-01, 4.44212833419920E-03,
     1-8.42328522190089E-05,-2.58040318418710E-06, 3.42389720217621E-07,
     2-6.24286894709776E-09,-2.36377836844577E-09, 3.16991042656673E-10,
     3-4.40995691658191E-12,-5.18674221093575E-12, 9.64874015137022E-13,
     4-4.90190576608710E-14,-1.77253430678112E-14, 5.55950610442662E-15,
     5-7.11793337579530E-16/
      IF (X.LT.0.) GO TO 90
      IF (C.GT.5.) GO TO 60
      IF (X.GT.1.2) GO TO 30
      T = (X+X-1.2)*CON4
      TT = T + T
      J = N1
      F1 = AK1(J)
      F2 = 0.
      DO 10 I=1,M1
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + AK1(J)
        F2 = TEMP1
   10 CONTINUE
      AI = T*F1 - F2 + AK1(1)
      J = N1D
      F1 = DAK1(J)
      F2 = 0.
      DO 20 I=1,M1D
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + DAK1(J)
        F2 = TEMP1
   20 CONTINUE
      DAI = -(T*F1-F2+DAK1(1))
      RETURN
   30 CONTINUE
      T = (X+X-CON2)*CON3
      TT = T + T
      J = N2
      F1 = AK2(J)
      F2 = 0.
      DO 40 I=1,M2
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + AK2(J)
        F2 = TEMP1
   40 CONTINUE
      RTRX = SQRT(RX)
      EC = EXP(-C)
      AI = EC*(T*F1-F2+AK2(1))/RTRX
      J = N2D
      F1 = DAK2(J)
      F2 = 0.
      DO 50 I=1,M2D
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + DAK2(J)
        F2 = TEMP1
   50 CONTINUE
      DAI = -EC*(T*F1-F2+DAK2(1))*RTRX
      RETURN
   60 CONTINUE
      T = 10./C - 1.
      TT = T + T
      J = N1
      F1 = AK3(J)
      F2 = 0.
      DO 70 I=1,M1
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + AK3(J)
        F2 = TEMP1
   70 CONTINUE
      RTRX = SQRT(RX)
      EC = EXP(-C)
      AI = EC*(T*F1-F2+AK3(1))/RTRX
      J = N1D
      F1 = DAK3(J)
      F2 = 0.
      DO 80 I=1,M1D
        J = J - 1
        TEMP1 = F1
        F1 = TT*F1 - F2 + DAK3(J)
        F2 = TEMP1
   80 CONTINUE
      DAI = -RTRX*EC*(T*F1-F2+DAK3(1))
      RETURN
   90 CONTINUE
      IF (C.GT.5.) GO TO 120
      T = .4*C - 1.
      TT = T + T
      J = N3
      F1 = AJP(J)
      E1 = AJN(J)
      F2 = 0.
      E2 = 0.
      DO 100 I=1,M3
        J = J - 1
        TEMP1 = F1
        TEMP2 = E1
        F1 = TT*F1 - F2 + AJP(J)
        E1 = TT*E1 - E2 + AJN(J)
        F2 = TEMP1
        E2 = TEMP2
  100 CONTINUE
      AI = (T*E1-E2+AJN(1)) - X*(T*F1-F2+AJP(1))
      J = N3D
      F1 = DAJP(J)
      E1 = DAJN(J)
      F2 = 0.
      E2 = 0.
      DO 110 I=1,M3D
        J = J - 1
        TEMP1 = F1
        TEMP2 = E1
        F1 = TT*F1 - F2 + DAJP(J)
        E1 = TT*E1 - E2 + DAJN(J)
        F2 = TEMP1
        E2 = TEMP2
  110 CONTINUE
      DAI = X*X*(T*F1-F2+DAJP(1)) + (T*E1-E2+DAJN(1))
      RETURN
  120 CONTINUE
      T = 10./C - 1.
      TT = T + T
      J = N4
      F1 = A(J)
      E1 = B(J)
      F2 = 0.
      E2 = 0.
      DO 130 I=1,M4
        J = J - 1
        TEMP1 = F1
        TEMP2 = E1
        F1 = TT*F1 - F2 + A(J)
        E1 = TT*E1 - E2 + B(J)
        F2 = TEMP1
        E2 = TEMP2
  130 CONTINUE
      TEMP1 = T*F1 - F2 + A(1)
      TEMP2 = T*E1 - E2 + B(1)
      RTRX = SQRT(RX)
      CV = C - FPI12
      CCV = COS(CV)
      SCV = SIN(CV)
      AI = (TEMP1*CCV-TEMP2*SCV)/RTRX
      J = N4D
      F1 = DA(J)
      E1 = DB(J)
      F2 = 0.
      E2 = 0.
      DO 140 I=1,M4D
        J = J - 1
        TEMP1 = F1
        TEMP2 = E1
        F1 = TT*F1 - F2 + DA(J)
        E1 = TT*E1 - E2 + DB(J)
        F2 = TEMP1
        E2 = TEMP2
  140 CONTINUE
      TEMP1 = T*F1 - F2 + DA(1)
      TEMP2 = T*E1 - E2 + DB(1)
      E1 = CCV*CON5 + .5*SCV
      E2 = SCV*CON5 - .5*CCV
      DAI = (TEMP1*E1-TEMP2*E2)*RTRX
      RETURN
      END
      SUBROUTINE JBESS(ALPHA, N, X, Y)                                  JBE   10
C
C         A CDC 6600 SUBROUTINE
C
C     AUTHORS
C         D.E. AMOS, S.L. DANIEL, AND M.K. WESTON
C         SANDIA LABORATORIES
C         JANUARY, 1975
C
C     REFERENCES
C         ABRAMOWITZ, M. AND STEGUN, I.A. HANDBOOK OF MATHEMATICAL
C         FUNCTIONS. NBS APPLIED MATHEMATICS SERIES 55, U.S. GOVERNMENT
C         PRINTING OFFICE, WASHINGTON, D.C., CHAPTERS 9 AND 10.
C
C         AMOS, D.E., DANIEL, S.L. AND WESTON, M.K. CDC 6600
C         SUBROUTINES IBESS AND JBESS FOR BESSEL FUNCTIONS
C         I/SUB(NU)/(X) AND J/SUB(NU)/(X), X.GE.0, NU.GE.0.
C         ACM TRANS. MATH. SOFTWARE, MARCH, 1977.
C
C         OLVER, F.W.J. TABLES FOR BESSEL FUNCTIONS OF MODERATE OR
C         LARGE ORDERS. NPL MATHEMATICAL TABLES, VOL 6. HER MAJESTY-S
C         STATIONERY OFFICE, LONDON, 1962.
C
C         OLVER, F.W.J. THE ASYMPTOTIC EXPANSION OF BESSEL FUNCTIONS OF
C         LARGE ORDER. PHIL. TRANS. A, 247, PP. 328-368, 1954.
C
C     ABSTRACT
C         JBESS COMPUTES AN N MEMBER SEQUENCE OF J BESSEL FUNCTIONS
C         J/SUB(ALPHA+K-1)/(X), K=1,...,N FOR NON-NEGATIVE ALPHA AND X.
C         A COMBINATION OF THE POWER SERIES, THE ASYMPTOTIC EXPANSION
C         FOR X TO INFINITY AND THE UNIFORM ASYMPTOTIC EXPANSION FOR
C         NU TO INFINITY ARE APPLIED OVER SUBDIVISIONS OF THE (NU,X)
C         PLANE. FOR VALUES OF (NU,X) NOT COVERED BY ONE OF THESE
C         FORMULAE, THE ORDER IS INCREMENTED OR DECREMENTED BY INTEGER
C         VALUES INTO A REGION WHERE ONE OF THE FORMULAE APPLY. BACKWARD
C         RECURSION IS APPLIED TO REDUCE ORDERS BY INTEGER VALUES EXCEPT
C         WHERE THE ENTIRE SEQUENCE LIES IN THE OSCILLATORY REGION. IN
C         THIS CASE FORWARD RECURSION IS STABLE AND VALUES FROM THE
C         ASYMPTOTIC EXPANSION FOR X TO INFINITY START THE RECURSION
C         WHEN IT IS EFFICIENT TO DO SO. LEADING TERMS OF THE SERIES AND
C         UNIFORM EXPANSION ARE TESTED FOR UNDERFLOW. IF A SEQUENCE IS
C         REQUESTED AND THE LAST MEMBER WOULD UNDERFLOW, THE RESULT IS
C         SET TO ZERO AND THE NEXT LOWER ORDER TRIED, ETC., UNTIL A
C         MEMBER COMES ON SCALE OR ALL MEMBERS ARE SET TO ZERO. OVERFLOW
C         CANNOT OCCUR. JBESS CALLS SUBROUTINE JAIRY AND FUNCTION GAMLN.
C
C     DESCRIPTION OF ARGUMENTS
C
C         INPUT
C           ALPHA  - ORDER OF FIRST MEMBER OF THE SEQUENCE, ALPHA.GE.0
C           N      - NUMBER OF MEMBERS IN THE SEQUENCE, N.GE.1
C           X      - X.GE.0
C
C         OUTPUT
C           Y      - A VECTOR WHOSE FIRST N COMPONENTS CONTAIN
C                    VALUES FOR J/SUB(ALPHA+K-1)/(X), K=1,...,N
C
C     ERROR CONDITIONS
C         IMPROPER INPUT ARGUMENTS - A FATAL ERROR
C         UNDERFLOW  - A NON-FATAL ERROR
C
      DOUBLE PRECISION DX, TRX, DTM, DFN
      DIMENSION Y(1)
      DIMENSION C(11,10), ALFA(26,4), BETA(26,5)
      DIMENSION C1(11,8), C2(11,2), ALFA1(26,2), ALFA2(26,2)
      DIMENSION BETA1(26,2), BETA2(26,2), BETA3(26,1)
      DIMENSION GAMA(26), TEMP(3), KMAX(5), AR(8), BR(10), UPOL(10)
      DIMENSION FNULIM(2), PP(4)
      DIMENSION CR(10), DR(10)
      EQUIVALENCE (C(1,1),C1(1,1))
      EQUIVALENCE (C(1,9),C2(1,1))
      EQUIVALENCE (ALFA(1,1),ALFA1(1,1))
      EQUIVALENCE (ALFA(1,3),ALFA2(1,1))
      EQUIVALENCE (BETA(1,1),BETA1(1,1))
      EQUIVALENCE (BETA(1,3),BETA2(1,1))
      EQUIVALENCE (BETA(1,5),BETA3(1,1))
      DATA ELIM1,ELIM2,TOL /   667.    ,    644.    ,    1.E-15        /
C
      DATA PP              / 8.72909153935547E+00, 2.65693932265030E-01,
     1 1.24578576865586E-01, 7.70133747430388E-04/
C TOLS=LN(1.E-3)
      DATA TOLS            /-6.90775527898214E+00/
      DATA CON1,CON2,CON548/ 6.66666666666667E-01, 3.33333333333333E-01,
     1 1.04166666666667E-01/
      DATA RTWO,PDF,RTTP,PIDT                    / 1.34839972492648E+00,
     1 7.85398163397448E-01, 7.97884560802865E-01, 1.57079632679490E+00/
      DATA FNULIM          /         100.        ,         60.         /
C CE=-ALOG(TOL) , TCE=-0.75*ALOG(TOL)
      DATA CE , TCE        / 3.45387763949107E+01, 2.59040822961830E+01/
      DATA INLIM           /         150         /
      DATA   AR            / 8.35503472222222E-02, 1.28226574556327E-01,
     1 2.91849026464140E-01, 8.81627267443758E-01, 3.32140828186277E+00,
     2 1.49957629868626E+01, 7.89230130115865E+01, 4.74451538868264E+02/
      DATA   BR            /-1.45833333333333E-01,-9.87413194444444E-02,
     1-1.43312053915895E-01,-3.17227202678414E-01,-9.42429147957120E-01,
     2-3.51120304082635E+00,-1.57272636203680E+01,-8.22814390971859E+01,
     3-4.92355370523671E+02,-3.31621856854797E+03/
      DATA C1              /-2.08333333333333E-01, 1.25000000000000E-01,
     1         9(0.)       , 3.34201388888889E-01,-4.01041666666667E-01,
     2 7.03125000000000E-02,         8(0.)       ,-1.02581259645062E+00,
     3 1.84646267361111E+00,-8.91210937500000E-01, 7.32421875000000E-02,
     4         7(0.)       , 4.66958442342625E+00,-1.12070026162230E+01,
     5 8.78912353515625E+00,-2.36408691406250E+00, 1.12152099609375E-01,
     6         6(0.)       ,-2.82120725582002E+01, 8.46362176746007E+01,
     7-9.18182415432400E+01, 4.25349987453885E+01,-7.36879435947963E+00,
     8 2.27108001708984E-01,         5(0.)       , 2.12570130039217E+02,
     9-7.65252468141182E+02, 1.05999045252800E+03,-6.99579627376133E+02,
     A 2.18190511744212E+02,-2.64914304869516E+01, 5.72501420974731E-01,
     B         4(0.)       ,-1.91945766231841E+03, 8.06172218173731E+03,
     C-1.35865500064341E+04, 1.16553933368645E+04,-5.30564697861340E+03,
     D 1.20090291321635E+03,-1.08090919788395E+02, 1.72772750258446E+00,
     E         3(0.)       , 2.02042913309661E+04,-9.69805983886375E+04,
     F 1.92547001232532E+05,-2.03400177280416E+05, 1.22200464983017E+05,
     G-4.11926549688976E+04, 7.10951430248936E+03,-4.93915304773088E+02,
     H 6.07404200127348E+00,         2(0.)       /
      DATA C2              /-2.42919187900551E+05, 1.31176361466298E+06,
     1-2.99801591853811E+06, 3.76327129765640E+06,-2.81356322658653E+06,
     2 1.26836527332162E+06,-3.31645172484564E+05, 4.52187689813627E+04,
     3-2.49983048181121E+03, 2.43805296995561E+01,         1(0.)       ,
     4 3.28446985307204E+06,-1.97068191184322E+07, 5.09526024926646E+07,
     5-7.41051482115327E+07, 6.63445122747290E+07,-3.75671766607634E+07,
     6 1.32887671664218E+07,-2.78561812808645E+06, 3.08186404612662E+05,
     7-1.38860897537170E+04, 1.10017140269247E+02/
      DATA ALFA1           /-4.44444444444444E-03,-9.22077922077922E-04,
     1-8.84892884892885E-05, 1.65927687832450E-04, 2.46691372741793E-04,
     2 2.65995589346255E-04, 2.61824297061501E-04, 2.48730437344656E-04,
     3 2.32721040083232E-04, 2.16362485712365E-04, 2.00738858762752E-04,
     4 1.86267636637545E-04, 1.73060775917876E-04, 1.61091705929016E-04,
     5 1.50274774160908E-04, 1.40503497391270E-04, 1.31668816545923E-04,
     6 1.23667445598253E-04, 1.16405271474738E-04, 1.09798298372713E-04,
     7 1.03772410422993E-04, 9.82626078369363E-05, 9.32120517249503E-05,
     8 8.85710852478712E-05, 8.42963105715700E-05, 8.03497548407791E-05,
     9 6.93735541354589E-04, 2.32241745182922E-04,-1.41986273556691E-05,
     A-1.16444931672049E-04,-1.50803558053049E-04,-1.55121924918096E-04,
     B-1.46809756646466E-04,-1.33815503867491E-04,-1.19744975684254E-04,
     C-1.06184319207974E-04,-9.37699549891194E-05,-8.26923045588193E-05,
     D-7.29374348155221E-05,-6.44042357721016E-05,-5.69611566009369E-05,
     E-5.04731044303562E-05,-4.48134868008883E-05,-3.98688727717599E-05,
     F-3.55400532972042E-05,-3.17414256609022E-05,-2.83996793904175E-05,
     G-2.54522720634871E-05,-2.28459297164725E-05,-2.05352753106481E-05,
     H-1.84816217627666E-05,-1.66519330021394E-05/
      DATA ALFA2           /-3.54211971457744E-04,-1.56161263945159E-04,
     1 3.04465503594936E-05, 1.30198655773243E-04, 1.67471106699712E-04,
     2 1.70222587683593E-04, 1.56501427608595E-04, 1.36339170977445E-04,
     3 1.14886692029825E-04, 9.45869093034688E-05, 7.64498419250898E-05,
     4 6.07570334965197E-05, 4.74394299290509E-05, 3.62757512005344E-05,
     5 2.69939714979225E-05, 1.93210938247939E-05, 1.30056674793963E-05,
     6 7.82620866744497E-06, 3.59257485819352E-06, 1.44040049814252E-07,
     7-2.65396769697939E-06,-4.91346867098486E-06,-6.72739296091248E-06,
     8-8.17269379678658E-06,-9.31304715093561E-06,-1.02011418798016E-05,
     9 3.78194199201773E-04, 2.02471952761816E-04,-6.37938506318862E-05,
     A-2.38598230603006E-04,-3.10916256027362E-04,-3.13680115247576E-04,
     B-2.78950273791323E-04,-2.28564082619141E-04,-1.75245280340847E-04,
     C-1.25544063060690E-04,-8.22982872820208E-05,-4.62860730588116E-05,
     D-1.72334302366962E-05, 5.60690482304602E-06, 2.31395443148287E-05,
     E 3.62642745856794E-05, 4.58006124490189E-05, 5.24595294959114E-05,
     F 5.68396208545815E-05, 5.94349820393104E-05, 6.06478527578422E-05,
     G 6.08023907788436E-05, 6.01577894539460E-05, 5.89199657344698E-05,
     H 5.72515823777593E-05, 5.52804375585853E-05/
      DATA BETA1           / 1.79988721413553E-02, 5.59964911064388E-03,
     1 2.88501402231133E-03, 1.80096606761054E-03, 1.24753110589199E-03,
     2 9.22878876572938E-04, 7.14430421727287E-04, 5.71787281789705E-04,
     3 4.69431007606482E-04, 3.93232835462917E-04, 3.34818889318298E-04,
     4 2.88952148495752E-04, 2.52211615549573E-04, 2.22280580798883E-04,
     5 1.97541838033063E-04, 1.76836855019718E-04, 1.59316899661821E-04,
     6 1.44347930197334E-04, 1.31448068119965E-04, 1.20245444949303E-04,
     7 1.10449144504599E-04, 1.01828770740567E-04, 9.41998224204238E-05,
     8 8.74130545753834E-05, 8.13466262162801E-05, 7.59002269646219E-05,
     9-1.49282953213429E-03,-8.78204709546389E-04,-5.02916549572035E-04,
     A-2.94822138512746E-04,-1.75463996970783E-04,-1.04008550460816E-04,
     B-5.96141953046458E-05,-3.12038929076098E-05,-1.26089735980230E-05,
     C-2.42892608575730E-07, 8.05996165414274E-06, 1.36507009262147E-05,
     D 1.73964125472926E-05, 1.98672978842134E-05, 2.14463263790823E-05,
     E 2.23954659232457E-05, 2.28967783814713E-05, 2.30785389811178E-05,
     F 2.30321976080909E-05, 2.28236073720349E-05, 2.25005881105292E-05,
     G 2.20981015361991E-05, 2.16418427448104E-05, 2.11507649256221E-05,
     H 2.06388749782171E-05, 2.01165241997082E-05/
      DATA BETA2           / 5.52213076721293E-04, 4.47932581552385E-04,
     1 2.79520653992021E-04, 1.52468156198447E-04, 6.93271105657044E-05,
     2 1.76258683069991E-05,-1.35744996343269E-05,-3.17972413350427E-05,
     3-4.18861861696693E-05,-4.69004889379141E-05,-4.87665447413787E-05,
     4-4.87010031186735E-05,-4.74755620890087E-05,-4.55813058138628E-05,
     5-4.33309644511266E-05,-4.09230193157750E-05,-3.84822638603221E-05,
     6-3.60857167535411E-05,-3.37793306123367E-05,-3.15888560772110E-05,
     7-2.95269561750807E-05,-2.75978914828336E-05,-2.58006174666884E-05,
     8-2.41308356761280E-05,-2.25823509518346E-05,-2.11479656768913E-05,
     9-4.74617796559960E-04,-4.77864567147321E-04,-3.20390228067038E-04,
     A-1.61105016119962E-04,-4.25778101285435E-05, 3.44571294294968E-05,
     B 7.97092684075675E-05, 1.03138236708272E-04, 1.12466775262204E-04,
     C 1.13103642108481E-04, 1.08651634848774E-04, 1.01437951597662E-04,
     D 9.29298396593364E-05, 8.40293133016090E-05, 7.52727991349134E-05,
     E 6.69632521975731E-05, 5.92564547323195E-05, 5.22169308826976E-05,
     F 4.58539485165361E-05, 4.01445513891487E-05, 3.50481730031328E-05,
     G 3.05157995034347E-05, 2.64956119950516E-05, 2.29363633690998E-05,
     H 1.97893056664022E-05, 1.70091984636413E-05/
      DATA BETA3           / 7.36465810572578E-04, 8.72790805146194E-04,
     1 6.22614862573135E-04, 2.85998154194304E-04, 3.84737672879366E-06,
     2-1.87906003636972E-04,-2.97603646594555E-04,-3.45998126832656E-04,
     3-3.53382470916038E-04,-3.35715635775049E-04,-3.04321124789040E-04,
     4-2.66722723047613E-04,-2.27654214122820E-04,-1.89922611854562E-04,
     5-1.55058918599094E-04,-1.23778240761874E-04,-9.62926147717644E-05,
     6-7.25178327714425E-05,-5.22070028895634E-05,-3.50347750511901E-05,
     7-2.06489761035552E-05,-8.70106096849767E-06, 1.13698686675100E-06,
     8 9.16426474122779E-06, 1.56477785428873E-05, 2.08223629482467E-05/
      DATA GAMA            / 6.29960524947437E-01, 2.51984209978975E-01,
     1 1.54790300415656E-01, 1.10713062416159E-01, 8.57309395527395E-02,
     2 6.97161316958684E-02, 5.86085671893714E-02, 5.04698873536311E-02,
     3 4.42600580689155E-02, 3.93720661543510E-02, 3.54283195924455E-02,
     4 3.21818857502098E-02, 2.94646240791158E-02, 2.71581677112934E-02,
     5 2.51768272973862E-02, 2.34570755306079E-02, 2.19508390134907E-02,
     6 2.06210828235646E-02, 1.94388240897881E-02, 1.83810633800683E-02,
     7 1.74293213231963E-02, 1.65685837786612E-02, 1.57865285987918E-02,
     8 1.50729501494096E-02, 1.44193250839955E-02, 1.38184805735342E-02/
C TEST INPUT ARGUMENTS
      KT = 1
      IF (N-1) 710, 10, 20
   10 KT = 2
   20 NN = N
      IF (X) 720, 30, 80
   30 IF (ALPHA) 700, 40, 50
   40 Y(1) = 1.
      IF (N.EQ.1) RETURN
      I1 = 2
      GO TO 60
   50 I1 = 1
   60 DO 70 I=I1,N
        Y(I) = 0.
   70 CONTINUE
      RETURN
   80 CONTINUE
      IF (ALPHA.LT.0.) GO TO 700
      DFN = DBLE(FLOAT(N)) + DBLE(ALPHA) - 1.D+0
      FNU = DFN
      XO2 = X*.5
      SXO2 = XO2*XO2
C DECISION TREE FOR REGION WHERE SERIES, ASYMPTOTIC EXPANSION FOR X
C TO INFINITY AND ASYMPTOTIC EXPANSION FOR NU TO INFINITY ARE
C APPLIED.
      IF (SXO2.LE.(FNU+1.)) GO TO 90
      TA = AMAX1(20.,FNU)
      IF (X.GT.TA) GO TO 120
      IF (X.GT.12.) GO TO 110
      XO2L = ALOG(XO2)
      NS = SXO2 - FNU
      GO TO 100
   90 FN = FNU
      FNP1 = FN + 1.
      XO2L = ALOG(XO2)
      IS = KT
      IF (X.LE.0.5) GO TO 330
      NS = 0
  100 DFN = DFN + DBLE(FLOAT(NS))
      FN = DFN
      FNP1 = FN + 1.
      IS = KT
      IF (N-1+NS.GT.0) IS = 3
      GO TO 330
  110 NS = AMAX1(36.-FNU,0.)
      DFN = DFN + DBLE(FLOAT(NS))
      FN = DFN
      IS = KT
      IF (N-1+NS.GT.0) IS = 3
      GO TO 130
  120 CONTINUE
      RTX = SQRT(X)
      TAU = RTWO*RTX
      TA = TAU + FNULIM(KT)
      IF (FNU.LE.TA) GO TO 480
      FN = FNU
      IS = KT
C UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY
  130 CONTINUE
      XX = X/FN
      W2 = 1. - XX*XX
      ABW2 = ABS(W2)
      RA = SQRT(ABW2)
      IF (ABW2.GT.0.2775) GO TO 220
C CASES NEAR X=FN, ABS(1.-(X/FN)**2).LE.0.2775
C COEFFICIENTS OF ASYMPTOTIC EXPANSION BY SERIES
C ZETA AND TRUNCATION FOR A(ZETA) AND B(ZETA) SERIES
C KMAX IS TRUNCATION INDEX FOR A(ZETA) AND B(ZETA) SERIES=MAX(2,SA)
      SA = 0.
      IF (ABW2.EQ.0.) GO TO 140
      SA = TOLS/ALOG(ABW2)
  140 SB = SA
      DO 150 I=1,5
        KMAX(I) = AMAX1(SA,2.)
        SA = SA + SB
  150 CONTINUE
      KB = KMAX(5)
      KLAST = KB - 1
      SA = GAMA(KB)
      DO 160 K=1,KLAST
        KB = KB - 1
        SA = SA*W2 + GAMA(KB)
  160 CONTINUE
      Z = W2*SA
      AZ = ABS(Z)
      RTZ = SQRT(AZ)
      FN13 = FN**CON2
      RTARY = RTZ*FN13
      ARY = -RTARY*RTARY
      AZ32 = AZ*RTZ*CON1
      ACZ = FN*AZ32
      IF (Z.LE.0.) GO TO 170
C TEST FOR UNDERFLOW, 1.E-280=EXP(-644.), ONE WORD LENGTH
C UP FROM UNDERFLOW LIMIT OF CDC 6600
      IF (ACZ.GT.ELIM2) GO TO 380
      ARY = -ARY
  170 PHI = SQRT(SQRT(SA+SA+SA+SA))
C B(ZETA) FOR S=0
      KB = KMAX(5)
      KLAST = KB - 1
      SB = BETA(KB,1)
      DO 180 K=1,KLAST
        KB = KB - 1
        SB = SB*W2 + BETA(KB,1)
  180 CONTINUE
      KSP1 = 1
      FN2 = FN*FN
      RFN2 = 1./FN2
      RDEN = 1.
      ASUM = 1.
      RELB = TOL*ABS(SB)
      BSUM = SB
      DO 200 KS=1,4
        KSP1 = KSP1 + 1
        RDEN = RDEN*RFN2
C A(ZETA) AND B(ZETA) FOR S=1,2,3,4
        KB = KMAX(5-KS)
        KLAST = KB - 1
        SA = ALFA(KB,KS)
        SB = BETA(KB,KSP1)
        DO 190 K=1,KLAST
          KB = KB - 1
          SA = SA*W2 + ALFA(KB,KS)
          SB = SB*W2 + BETA(KB,KSP1)
  190   CONTINUE
        TA = SA*RDEN
        TB = SB*RDEN
        ASUM = ASUM + TA
        BSUM = BSUM + TB
        IF (ABS(TA).LE.TOL .AND. ABS(TB).LE.RELB) GO TO 210
  200 CONTINUE
  210 CONTINUE
      BSUM = BSUM/(FN*FN13)
      GO TO 300
  220 CONTINUE
      TAU = 1./RA
      T2 = 1./W2
      IF (W2.GE.0.) GO TO 230
C CASES FOR (X/FN).GT.SQRT(1.2775)
      AZ32 = ABS(RA-ATAN(RA))
      ACZ = AZ32*FN
      CZ = -ACZ
      Z32 = 1.5*AZ32
      RTZ = Z32**CON2
      FN13 = FN**CON2
      RTARY = RTZ*FN13
      ARY = -RTARY*RTARY
      GO TO 240
  230 CONTINUE
C CASES FOR (X/FN).LT.SQRT(0.7225)
      AZ32 = ABS(ALOG((1.+RA)/XX)-RA)
C TEST FOR UNDERFLOW, 1.E-280 = EXP(-644.), ONE WORD LENGTH
C UP FROM UNDERFLOW LIMIT OF CDC 6600
      ACZ = AZ32*FN
      CZ = ACZ
      IF (ACZ.GT.ELIM2) GO TO 380
      Z32 = 1.5*AZ32
      RTZ = Z32**CON2
      FN13 = FN**CON2
      RTARY = RTZ*FN13
      ARY = RTARY*RTARY
  240 CONTINUE
      PHI = SQRT((RTZ+RTZ)*TAU)
      TB = 1.
      ASUM = 1.
      TFN = TAU/FN
      UPOL(1) = 1.
      UPOL(2) = (C(1,1)*T2+C(2,1))*TFN
      RCZ = CON1/CZ
      CRZ32 = CON548*RCZ
      BSUM = UPOL(2) + CRZ32
      RELB = TOL*ABS(BSUM)
      AP = TFN
      KS = 0
      KP1 = 2
      RZDEN = RCZ
      DO 280 LR=2,8,2
C COMPUTE TWO U POLYNOMIALS FOR NEXT A(ZETA) AND B(ZETA)
        LRP1 = LR + 1
        DO 260 K=LR,LRP1
          KS = KS + 1
          KP1 = KP1 + 1
          S1 = C(1,K)
          DO 250 J=2,KP1
            S1 = S1*T2 + C(J,K)
  250     CONTINUE
          AP = AP*TFN
          UPOL(KP1) = AP*S1
          CR(KS) = BR(KS)*RZDEN
          RZDEN = RZDEN*RCZ
          DR(KS) = AR(KS)*RZDEN
  260   CONTINUE
        SUMA = UPOL(LRP1)
        SUMB = UPOL(LR+2) + UPOL(LRP1)*CRZ32
        JU = LRP1
        DO 270 JR=1,LR
          JU = JU - 1
          SUMA = SUMA + CR(JR)*UPOL(JU)
          SUMB = SUMB + DR(JR)*UPOL(JU)
  270   CONTINUE
        TB = -TB
        IF (W2.GT.0.) TB = ABS(TB)
        ASUM = ASUM + SUMA*TB
        BSUM = BSUM + SUMB*TB
        IF (ABS(SUMA).LE.TOL .AND. ABS(SUMB).LE.RELB) GO TO 290
  280 CONTINUE
  290 TB = RTARY
      IF (W2.GT.0.) TB = -TB
      BSUM = BSUM/TB
  300 CONTINUE
      CALL JAIRY(ARY, RTARY, ACZ, AI, DAI)
      TEMP(IS) = PHI*(AI*ASUM+DAI*BSUM)/FN13
      GO TO (320, 450, 610), IS
  310 TEMP(1) = TEMP(3)
      KT = 1
  320 IS = 2
      DFN = DFN - 1.D+0
      FN = DFN
      GO TO 130
C SERIES FOR (X/2)**2.LE.NU+1
  330 CONTINUE
      GLN = GAMLN(FNP1)
      ARG = FN*XO2L - GLN
      IF (ARG.LT.-ELIM1) GO TO 400
      EARG = EXP(ARG)
  340 CONTINUE
      S = 1.
      AK = 3.
      T2 = 1.
      T = 1.
      S1 = FN
      DO 350 K=1,17
        S2 = T2 + S1
        T = -T*SXO2/S2
        S = S + T
        IF (ABS(T).LT.TOL) GO TO 360
        T2 = T2 + AK
        AK = AK + 2.
        S1 = S1 + FN
  350 CONTINUE
  360 CONTINUE
      TEMP(IS) = S*EARG
      GO TO (370, 450, 600), IS
  370 EARG = EARG*FN/XO2
      DFN = DFN - 1.D+0
      FN = DFN
      IS = 2
      GO TO 340
C SET UNDERFLOW VALUE AND UPDATE PARAMETERS
  380 Y(NN) = 0.
      NN = NN - 1
      DFN = DFN - 1.D+0
      FN = DFN
      IF (NN-1) 440, 390, 130
  390 KT = 2
      IS = 2
      GO TO 130
  400 Y(NN) = 0.
      NN = NN - 1
      FNP1 = FN
      DFN = DFN - 1.D+0
      FN = DFN
      IF (NN-1) 440, 410, 420
  410 KT = 2
      IS = 2
  420 IF (SXO2.LE.FNP1) GO TO 430
      GO TO 130
  430 ARG = ARG - XO2L + ALOG(FNP1)
      IF (ARG.LT.-ELIM1) GO TO 400
      GO TO 330
  440 NZ = N - NN
      PRINT 99996, NZ, ALPHA, N, X
      RETURN
C BACKWARD RECURSION SECTION
  450 CONTINUE
      NZ = N - NN
      IF (NZ.NE.0) PRINT 99996, NZ, ALPHA, N, X
      IF (KT.EQ.2) GO TO 470
C BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
      Y(NN) = TEMP(1)
      Y(NN-1) = TEMP(2)
      IF (NN.EQ.2) RETURN
      DX = X
      TRX = 2.D+0/DX
      DTM = DFN*TRX
      TM = DTM
      K = NN + 1
      DO 460 I=3,NN
        K = K - 1
        Y(K-2) = TM*Y(K-1) - Y(K)
        DTM = DTM - TRX
        TM = DTM
  460 CONTINUE
      RETURN
  470 Y(1) = TEMP(2)
      RETURN
C ASYMPTOTIC EXPANSION FOR X TO INFINITY WITH FORWARD RECURSION IN
C OSCILLATORY REGION X.GT.MAX(20, NU), PROVIDED THE LAST MEMBER
C OF THE SEQUENCE IS ALSO IN THE REGION.
  480 CONTINUE
      IN = ALPHA - TAU + 2.
      IF (IN.LE.0) GO TO 490
      INP1 = IN + 1
      DALPHA = ALPHA - FLOAT(INP1)
      KT = 1
      GO TO 500
  490 DALPHA = ALPHA
      IN = 0
  500 IS = KT
      ARG = X - PIDT*DALPHA - PDF
      SA = SIN(ARG)
      SB = COS(ARG)
      RA = RTTP/RTX
      ETX = 8.*X
  510 DX = DALPHA
      DX = DX + DX
      DTM = DX*DX
      T2 = DTM - 1.D+0
      T2 = T2/ETX
      S2 = T2
      RELB = TOL*ABS(T2)
      T1 = ETX
      S1 = 1.
      FN = 1.
      AK = 8.
      DO 520 K=1,13
        T1 = T1 + ETX
        FN = FN + AK
        DX = FN
        TRX = DTM - DX
        AP = TRX
        T2 = -T2*AP/T1
        S1 = S1 + T2
        T1 = T1 + ETX
        AK = AK + 8.
        FN = FN + AK
        DX = FN
        TRX = DTM - DX
        AP = TRX
        T2 = T2*AP/T1
        S2 = S2 + T2
        IF (ABS(T2).LE.RELB) GO TO 530
        AK = AK + 8.
  520 CONTINUE
  530 TEMP(IS) = RA*(S1*SB-S2*SA)
      GO TO (540, 550), IS
  540 DALPHA = DALPHA + 1.
      IS = 2
      TB = SA
      SA = -SB
      SB = TB
      GO TO 510
C FORWARD RECURSION SECTION
  550 IF (KT.EQ.2) GO TO 470
      S1 = TEMP(1)
      S2 = TEMP(2)
      TX = 2./X
      TM = DALPHA*TX
      IF (IN.EQ.0) GO TO 570
C FORWARD RECUR TO INDEX ALPHA
      DO 560 I=1,IN
        S = S2
        S2 = TM*S2 - S1
        TM = TM + TX
        S1 = S
  560 CONTINUE
      IF (NN.EQ.1) GO TO 590
      S = S2
      S2 = TM*S2 - S1
      TM = TM + TX
      S1 = S
  570 CONTINUE
C FORWARD RECUR FROM INDEX ALPHA TO ALPHA+N-1
      Y(1) = S1
      Y(2) = S2
      IF (NN.EQ.2) RETURN
      DO 580 I=3,NN
        Y(I) = TM*Y(I-1) - Y(I-2)
        TM = TM + TX
  580 CONTINUE
      RETURN
  590 Y(1) = S2
      RETURN
C BACKWARD RECURSION WITH NORMALIZATION BY
C ASYMPTOTIC EXPANSION FOR NU TO INFINITY OR POWER SERIES.
  600 CONTINUE
C COMPUTATION OF LAST ORDER FOR SERIES NORMALIZATION
      KM = AMAX1(3.-FN,0.)
      TFN = FN + FLOAT(KM)
      TA = (GLN+TFN-0.9189385332-0.0833333333/TFN)/(TFN+0.5)
      TA = XO2L - TA
      TB = -(1.-1.5/TFN)/TFN
      IN = CE/(-TA+SQRT(TA*TA-CE*TB)) + 1.5
      IN = IN + KM
      GO TO 650
  610 CONTINUE
C COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION
      GLN = AZ32 + RA
      IF (ARY.GT.30.) GO TO 630
      RDEN = (PP(4)*ARY+PP(3))*ARY + 1.
      RZDEN = PP(1) + PP(2)*ARY
      TA = RZDEN/RDEN
      IF (W2.LT.0.10) GO TO 620
      TB = GLN/RTARY
      GO TO 640
  620 TB = (1.259921049+0.1679894730*W2)/FN13
      GO TO 640
  630 CONTINUE
      TA = CON1*TCE/ACZ
      TA = ((0.0493827160*TA-0.1111111111)*TA+0.6666666667)*TA*ARY
      IF (W2.LT.0.10) GO TO 620
      TB = GLN/RTARY
  640 IN = TA/TB + 1.5
      IF (IN.GT.INLIM) GO TO 310
  650 DX = FLOAT(IN)
      DTM = DFN + DX
      DX = X
      TRX = 2.D+0/DX
      DTM = DTM*TRX
      TM = DTM
      TA = 0.
      TB = TOL
      KK = 1
  660 CONTINUE
C BACKWARD RECUR UNINDEXED
      DO 670 I=1,IN
        S = TB
        TB = TM*TB - TA
        TA = S
        DTM = DTM - TRX
        TM = DTM
  670 CONTINUE
C NORMALIZATION
      IF (KK.NE.1) GO TO 680
      TA = (TA/TB)*TEMP(3)
      TB = TEMP(3)
      KK = 2
      IN = NS
      IF (NS.NE.0) GO TO 660
  680 Y(NN) = TB
      NZ = N - NN
      IF (NZ.NE.0) PRINT 99996, NZ, ALPHA, N, X
      IF (NN.EQ.1) RETURN
      TB = TM*TB - TA
      DTM = DTM - TRX
      TM = DTM
      K = NN - 1
      Y(K) = TB
      IF (NN.EQ.2) RETURN
      KM = K - 1
C BACKWARD RECUR INDEXED
      DO 690 I=1,KM
        Y(K-1) = TM*Y(K) - Y(K+1)
        DTM = DTM - TRX
        TM = DTM
        K = K - 1
  690 CONTINUE
      RETURN
  700 PRINT 99999, ALPHA, N, X
      STOP
  710 PRINT 99998, ALPHA, N, X
      STOP
  720 PRINT 99997, ALPHA, N, X
      STOP
99999 FORMAT (51H0JBESS CALLED WITH THE ORDER, ALPHA, LESS THAN ZERO
     * /7H ALPHA=,E25.14,3H N=,I6,3H X=,E25.14)
99998 FORMAT (34H0JBESS CALLED WITH N LESS THAN ONE/7H ALPHA=,
     * E25.14,3H N=,I6,3H X=,E25.14)
99997 FORMAT (35H0JBESS CALLED WITH X LESS THAN ZERO/7H ALPHA=,
     * E25.14,3H N=,I6,3H X=,E25.14)
99996 FORMAT (25H0UNDERFLOW IN JBESS, LAST,I6,20H VALUE(S) OF Y ARRAY,
     * 17H WERE SET TO ZERO/7H ALPHA=,E25.14,3H N=,I6,3H X=,E25.14)
      END

.