      SUBROUTINE EXPLI(INT,ARG,RESULT,IERR)
C
      REAL A(6),B(5),C(7),D(7),E(7),F(7),P1(8),Q1(8),P2(8),Q2(7)
     1        ,P3(8),Q3(7),P4(8),Q4(7),P0(6),Q0(6),PX(9),QX(9)
     2        ,FRAC,SUMP,SUMQ,T,W,X,X0,XX0,XMX0,Y,DEXP40,XMAX
     3        ,XMIN,EI,ARG,RESULT
      INTEGER I,INT,IERR
      DOUBLE PRECISION DX0
C
C     THIS FUNPACK SUBROUTINE COMPUTES SINGLE PRECISION VALUES OF THE
C     EXPONENTIAL INTEGRALS
C
C            EI(X),  E-SUB-1(X) = -EI(-X),  AND  EXP(-X)*EI(X)
C
C     WHERE
C
C             INTEGRAL (FROM T=-INFINITY TO T=X) (EXP(T)/T),  X .GT. 0,
C     EI(X) =
C             -INTEGRAL (FROM T=-X TO T=INFINITY) (EXP(-T)/T),
C                                                         X .LT. 0,
C
C     AND WHERE THE FIRST INTEGRAL IS A PRINCIPAL VALUE INTEGRAL. THE
C     ARGUMENTS INT, ARG, AND RESULT HAVE THE FOLLOWING USAGE ...
C
C            INT              ARG             RESULT
C             1            X .NE. 0          EI(X)
C             2            X .GT. 0          E-SUB-1(X)
C             3            X .NE. 0          EXP(-X)*EI(X)
C
C     THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV
C     APPROXIMATIONS PUBLISHED IN MATH. COMP. 22, 641-649(1968), AND
C     MATH. COMP. 23, 289-303(1969) BY CODY AND THACHER.
C
C     OTHER SUBPROGRAMS REQUIRED
C
C               ABS, ALOG, DBLE, EXP
C
C        ------------
C
C     ERROR MONITORING
C
C        THE PARAMETER IERR IS A VARIABLE THAT IS SET BY THE ROUTINE.
C     IF NO ERRORS ARE DETECTED THEN IERR IS SET TO 0. THE FOLLOWING
C     TABLE INDICATES THE TYPES OF ERRORS THAT MAY BE ENCOUNTERED IN
C     THE ROUTINE AND THE FUNCTION VALUES SUPPLIED IN EACH CASE.
C
C     IERR    ERROR     ARGUMENT          FUNCTION VALUES FOR
C                        RANGE        EI(X)  EXP(-X)*EI(X)  E-SUB-1(X)
C      1    UNDERFLOW   X .LT. XMIN     0          -          0
C      2    OVERFLOW    X .GT. XMAX     T          -          -
C      3    ILLEGAL X     X = 0         T          T          T
C      4    ILLEGAL X    X .LT. 0       -          -          T
C
C     T INDICATES THAT THE ROUTINE TERMINATES WITHOUT ASSIGNING A VALUE
C     TO THE FUNCTION.
C
C        ----------
C
C     THIS SUBROUTINE WAS WRITTEN AT ARGONNE NATIONAL LABORATORY. THE
C     FORMAT OF THE SUBROUTINE WAS MODIFIED BY A.H. MORRIS AT NAVAL
C     SURFACE WEAPONS CENTER, DAHLGREN, VIRGINIA.
C
C     ----------------------------------------------------------------
C
C     XMIN AND XMAX ARE MACHINE DEPENDENT CONSTANTS FOR DETECTING
C     UNDERFLOW AND OVERFLOW (SEE THE DEFINITION OF IERR ABOVE).
C     IF B IS THE BASE OF THE FLOATING POINT ARITHMETIC BEING USED
C     AND EMIN AND EMAX THE SMALLEST AND LARGEST EXPONENTS
C     (SEE THE FUNCTION I1MACH), THEN THE FOLLOWING SETTINGS ARE
C     APPROPRIATE ...
C
C           B      EMIN      EMAX        XMIN         XMAX
C
C           2      -127       127        -84.2         88.0
C           2      -128       127        -84.9         88.0
C           8       -50        76       -101.4        163.1
C          16       -64        63       -175.0        179.8
C           2      -256       255       -172.9        181.9
C           2      -974      1070       -669.3        748.2
C           2     -8192      8191      -5670.3       5686.2
C
C     IN PARTICULAR, THE FOLLOWING SETTINGS ARE APPROPRIATE FOR THE
C     CDC 6000-7000 SERIES COMPUTERS ...
C
      DATA XMIN/-669.31/, XMAX/748.28/
C
C     ----------------------------------------------------------------
C
C     VALUE OF EXP(40.0)
C
C     ----------
      DATA DEXP40/.235385266837020E+18/
C     ----------
C
C     ZERO OF EI(X)
C
C     ----------
      DATA X0/.372507410781366/,DX0/.37250741078136663446199186658D0/
C     ----------
C
C     COEFFICIENTS FOR R(5,4) APPROXIMATION,
C     USED FOR -1.0 .LE. X .LT. 0.0
C
C     ----------
      DATA A(1)/-5.77215664901531321E-01/,
     1     A(2)/ 7.58833087029942888E-01/,
     2     A(3)/ 1.25660818982053153E-01/,
     3     A(4)/ 2.04158408934305020E-02/,
     4     A(5)/ 8.25035122466538239E-04/,
     5     A(6)/ 9.62949813453924099E-06/
      DATA B(1)/ 1.0E0/,
     1     B(2)/ 4.17810755380398068E-01/,
     2     B(3)/ 7.30228560396799462E-02/,
     3     B(4)/ 6.42720224671077755E-03/,
     4     B(5)/ 2.45134203588368906E-04/
C     ----------
C
C     COEFFICIENTS FOR R(6,6) APPROXIMATION,
C     USED FOR -4.0 .LE. X .LT. -1.0
C
C     ----------
      DATA C(1)/ 4.6562710797509566E-07/,
     1     C(2)/ 9.9997957705159497E-01/,
     2     C(3)/ 9.0416155694632866E+00/,
     3     C(4)/ 2.4378408879131673E+01/,
     4     C(5)/ 2.3019255939133346E+01/,
     5     C(6)/ 6.9052252278444357E+00/,
     6     C(7)/ 4.3096783946938878E-01/
      DATA D(1)/ 1.0E0/,
     1     D(2)/ 1.0041164382905448E+01/,
     2     D(3)/ 3.2426421069513805E+01/,
     3     D(4)/ 4.1280784189142434E+01/,
     4     D(5)/ 2.0449478501379417E+01/,
     5     D(6)/ 3.3190921359330160E+00/,
     6     D(7)/ 1.0340013040487398E-01/
C     ----------
C
C     COEFFICIENTS FOR R(6,6) APPROXIMATION,
C     USED FOR X .LT. -4.0
C
C     ----------
      DATA E(1)/-9.99999999998446910E-01/,
     1     E(2)/-2.66271060431811452E+01/,
     2     E(3)/-2.41055827097014847E+02/,
     3     E(4)/-8.95927957772936814E+02/,
     4     E(5)/-1.29885688746484101E+03/,
     5     E(6)/-5.45374158883132873E+02/,
     6     E(7)/-5.66575206533868735E+00/
      DATA F(1)/ 1.0E0/,
     1     F(2)/ 2.86271060422191896E+01/,
     2     F(3)/ 2.92310039388533247E+02/,
     3     F(4)/ 1.33278537748257231E+03/,
     4     F(5)/ 2.77761949509162962E+03/,
     5     F(6)/ 2.40401713225908947E+03/,
     6     F(7)/ 6.31657483280800229E+02/
C     ----------
C
C     COEFFICIENTS FOR R(5,5) APPROXIMATION FOR
C     LN(X/X0),  ABS(1-X/X0) .LT. .1
C
C     ----------
      DATA P0(1)/ 7.62694308249173762E+01/,
     1     P0(2)/ 1.70134167575098874E+02/,
     2     P0(3)/ 1.29869336902938160E+02/,
     3     P0(4)/ 3.89354857282565053E+01/,
     4     P0(5)/ 3.73025533400687462E+00/,
     5     P0(6)/ 2.78735890716569633E-02/
      DATA Q0(1)/ 7.62694308249173762E+01/,
     1     Q0(2)/ 2.08268882987557562E+02/,
     2     Q0(3)/ 2.08580634788411149E+02/,
     3     Q0(4)/ 9.28701998328388450E+01/,
     4     Q0(5)/ 1.74518115695293285E+01/,
     5     Q0(6)/ 1.0E+00/
C     ----------
C
C     COEFFICIENTS FOR R(7,7) APPROXIMATION,
C     IN CHEBYSHEV POLYNOMIAL FORM, USED FOR
C     0.0 .LT. X .LT. 6.0
C
C     ----------
      DATA P1(1)/-8.66937339951069561E+00/,
     1     P1(2)/-5.49142265521085145E+02/,
     2     P1(3)/-4.21001615357069934E+03/,
     3     P1(4)/-2.49301393458647588E+05/,
     4     P1(5)/-1.19623669349246868E+05/,
     5     P1(6)/-2.21744627758845378E+07/,
     6     P1(7)/ 3.89280421311201405E+06/,
     7     P1(8)/-1.95773036904547774E+08/
      DATA Q1(1)/ 3.41718750000000000E+01/,
     1     Q1(2)/-1.60708926587220852E+03/,
     2     Q1(3)/ 3.57300298058508106E+04/,
     3     Q1(4)/-4.83547436162163509E+05/,
     4     Q1(5)/ 4.28559624611749039E+06/,
     5     Q1(6)/-2.49033375740540327E+07/,
     6     Q1(7)/ 8.91925767575612109E+07/,
     7     Q1(8)/-8.26271498626054555E+07/
C     ----------
C
C     COEFFICIENTS FOR R(7,7) APPROXIMATION,
C     IN J-FRACTION FORM, USED FOR
C     6.0 .LE. X .LT. 12.0
C
C     ----------
      DATA P2(1)/-2.18086381520723706E+00/,
     1     P2(2)/-2.19010233854880688E+01/,
     2     P2(3)/ 9.30816385662165146E+00/,
     3     P2(4)/ 2.50762811293559830E+01/,
     4     P2(5)/-3.31842531997221118E+01/,
     5     P2(6)/ 6.01217990830080478E+01/,
     6     P2(7)/-4.32531132878134577E+01/,
     7     P2(8)/ 1.00443109228077910E+00/
      DATA Q2(1)/ 3.93707701852715001E+00/,
     1     Q2(2)/ 3.00892648372915199E+02/,
     2     Q2(3)/-6.25041161671875543E+00/,
     3     Q2(4)/ 1.00367439516725773E+03/,
     4     Q2(5)/ 1.43256738121937600E+01/,
     5     Q2(6)/ 2.73624119889328058E+03/,
     6     Q2(7)/ 5.27468851962907854E-01/
C     ----------
C
C     COEFFICIENTS FOR R(7,7) APPROXIMATION,
C     IN J-FRACTION FORM, USED FOR
C     12.0 .LE. X .LT. 24.0
C
C     ----------
      DATA P3(1)/-3.48334653602852613E+00/,
     1     P3(2)/-1.86545454883398835E+01/,
     2     P3(3)/-8.28561994140641315E+00/,
     3     P3(4)/-3.23467330305403459E+01/,
     4     P3(5)/ 1.79601688769251642E+01/,
     5     P3(6)/ 1.75656315469614422E+00/,
     6     P3(7)/-1.95022321289659817E+00/,
     7     P3(8)/ 9.99994296074708285E-01/
      DATA Q3(1)/ 6.95000655887433978E+01/,
     1     Q3(2)/ 5.72837193837323721E+01/,
     2     Q3(3)/ 2.57776384238439867E+01/,
     3     Q3(4)/ 7.60761148007734583E+02/,
     4     Q3(5)/ 2.89516727925135050E+01/,
     5     Q3(6)/-3.43942266899869969E+00/,
     6     Q3(7)/ 1.00083867402639122E+00/
C     ----------
C
C     COEFFICIENTS FOR R(7.7) APPROXIMATION,
C     IN J-FRACTION FORM, USED FOR X .GE. 24.0
C
C     ----------
      DATA P4(1)/-5.31686623494481619E+01/,
     1     P4(2)/ 8.91263822573707752E+00/,
     2     P4(3)/-1.39381360364405073E+00/,
     3     P4(4)/-3.08336269051762702E+01/,
     4     P4(5)/-7.49289167792884428E+00/,
     5     P4(6)/-5.00140345515924347E+00/,
     6     P4(7)/-3.00000016782085177E+00/,
     7     P4(8)/ 1.00000000000058393E+00/
      DATA Q4(1)/ 1.04745362652468301E+03/,
     1     Q4(2)/-6.74704580465832431E+00/,
     2     Q4(3)/ 2.95999399486831326E+02/,
     3     Q4(4)/-4.31325836146627961E+00/,
     4     Q4(5)/-7.90404992298925513E+00/,
     5     Q4(6)/-2.99996432944446454E+00/,
     6     Q4(7)/ 1.99999999924130892E+00/
C
C     ---------------------------------------------------------------
C
      X = ARG
      IERR = 0
      GO TO (100,450,100), INT
  100 IF (X) 280, 640, 110
  110 IF (X .GE. 12.E0) GO TO 200
      IF (X .GE. 6.E0) GO TO 160
C     ---------- 0.0 .LT. X .LT. 6.0.
C                RATIONAL APPROXIMATION USED IS EXPRESSED
C                IN TERMS OF CHEBYSHEV POLYNOMIALS TO
C                IMPROVE CONDITIONING  ----------
      T = X + X
      T = T / 3.0E0 - 2.0E0
      PX(1) = 0.0E0
      QX(1) = 0.0E0
      PX(2) = P1(1)
      QX(2) = Q1(1)
C
      DO 120 I = 2, 7
         PX(I+1) = T * PX(I) - PX(I-1) + P1(I)
         QX(I+1) = T * QX(I) - QX(I-1) + Q1(I)
  120 CONTINUE
C
      SUMP = .5E0 * T * PX(8) - PX(7) + P1(8)
      SUMQ = .5E0 * T * QX(8) - QX(7) + Q1(8)
      FRAC = SUMP / SUMQ
      XMX0 = DBLE(X) - DX0
      IF (ABS(XMX0) .LT. 0.037E0) GO TO 140
      XX0 = X / X0
      EI = ALOG(XX0) + XMX0 * FRAC
      IF (INT .EQ. 3) EI = EXP(-X) * EI
      GO TO 410
C     ---------- EVALUATE APPROXIMATION FOR LN(X/X0)
C                FOR X CLOSE TO X0 ----------
  140 Y = XMX0 / X0
      SUMP = ((((P0(6) * Y + P0(5)) * Y + P0(4)) * Y + P0(3)) * Y
     X          + P0(2)) * Y + P0(1)
      SUMQ = ((((Q0(6) * Y + Q0(5)) * Y + Q0(4)) * Y + Q0(3)) * Y
     X          + Q0(2)) * Y + Q0(1)
      EI = (SUMP / (SUMQ * X0) + FRAC) * XMX0
      IF (INT .EQ. 3) EI = EXP(-X) * EI
      GO TO 410
C     ---------- 6.0 .LE. X .LT. 12.0 ----------
  160 FRAC = 0.0E0
C
      DO 180 I = 1, 7
  180 FRAC = Q2(I) / (P2(I) + X + FRAC)
C
      EI = (P2(8) + FRAC) / X
      IF (INT .EQ. 3) GO TO 410
      EI = EI * EXP(X)
      GO TO 410
C     ---------- 12.0 .LE. X .LT. 24.0 ----------
  200 IF (X .GE. 24.E0) GO TO 240
      FRAC = 0.0E0
C
      DO 220 I = 1, 7
  220 FRAC = Q3(I) / (P3(I) + X + FRAC)
C
      EI = (P3(8) + FRAC) / X
      IF (INT .EQ. 3) GO TO 410
      EI = EI * EXP(X)
      GO TO 410
C     ---------- 24.0 .LE. X ----------
  240 IF ((X .GT. XMAX) .AND. (INT .LT. 3)) GO TO 620
      Y = 1.0E0 / X
      FRAC = 0.0E0
C
      DO 260 I = 1, 7
  260 FRAC = Q4(I) / (P4(I) + X + FRAC)
C
      FRAC = P4(8) + FRAC
      EI = Y + Y * Y * FRAC
      IF (INT .EQ. 3) GO TO 410
      IF (X .GT. 150.0E0) GO TO 270
      EI = EI * EXP(X)
      GO TO 410
C     ---------- CALCULATION REFORMULATED TO AVOID
C                PREMATURE OVERFLOW ----------
  270 EI = (EI * EXP(X-40.0E0)) * DEXP40
      GO TO 410
C     ---------- ORIGINAL X WAS NEGATIVE.  CALCULATION OF
C                E-SUB-1 JOINS AT LABEL 300 ----------
  280 Y = -X
  300 W = 1.0E0 / Y
      IF (Y .GT. 4.0E0) GO TO 340
      IF (Y .GT. 1.0E0) GO TO 320
C     ---------- 0.0 .LT. -X .LE. 1.0 ----------
      EI = ALOG(Y) - (((((A(6) * Y + A(5)) * Y + A(4))
     1            * Y + A(3)) * Y + A(2)) * Y + A(1)) /
     2       ((((B(5) * Y + B(4)) * Y + B(3))
     2            * Y + B(2)) * Y + B(1))
      IF (INT .EQ. 3) EI = EI * EXP(Y)
      GO TO 400
C     ---------- 1.0 .LT. -X .LE. 4.0 ----------
  320 EI = -((((((C(7) * W + C(6)) * W + C(5)) * W + C(4))
     1          * W + C(3)) * W + C(2)) * W + C(1))  /
     2       ((((((D(7) * W + D(6)) * W + D(5)) * W + D(4))
     2          * W + D(3)) * W + D(2)) * W + D(1))
      IF  (INT .EQ. 3) GO TO 410
      EI = EI * EXP(-Y)
      GO TO 400
C     ---------- 4.0 .LT. -X ----------
  340 IF ((-ABS(X) .LT. XMIN) .AND. (INT .LT. 3)) GO TO 600
      EI = -W * (1.0E0 + W * ((((((E(7)
     1        * W + E(6)) * W + E(5)) * W + E(4))
     2          * W + E(3)) * W + E(2)) * W + E(1)) /
     3       ((((((F(7) * W + F(6)) * W + F(5))
     4          * W + F(4)) * W + F(3)) * W + F(2)) * W + F(1)))
      IF (INT .EQ. 3) GO TO 410
      EI = EI * EXP(-Y)
  400 IF (INT .EQ. 2) EI = -EI
  410 RESULT = EI
      RETURN
  450 Y = X
      IF (Y) 660, 640, 300
C     ---------- ERROR RETURN FOR X .LT. XMIN,
C                CAUSING UNDERFLOW ----------
  600 EI = 0.0E0
      IERR = 1
      GO TO 410
C     ---------- ERROR RETURN FOR X .GT. XMAX,
C                CAUSING OVERFLOW ----------
  620 IERR = 2
      RETURN
C     ---------- ERROR RETURN FOR ILLEGAL
C                ARGUMENT, X = 0 ----------
  640 IERR = 3
      RETURN
C     ---------- ERROR RETURN FOR NEGATIVE
C                ARGUMENT IN E-SUB-1 ----------
  660 IERR = 4
      RETURN
      END
