C
C     ..................................................................
C
C        SUBROUTINE I0
C
C        PURPOSE
C            COMPUTE THE MODIFIED BESSEL FUNCTION I OF ORDER ZERO
C
C        USAGE
C            CALL I0(X,RI0)
C
C        DESCRIPTION OF PARAMETERS
C            X    -GIVEN ARGUMENT OF THE BESSEL FUNCTION I OF ORDER 0
C            RI0  -RESULTANT VALUE OF THE BESSEL FUNCTION I OF ORDER 0
C
C        REMARKS
C            LARGE VALUES OF THE ARGUMENT MAY CAUSE OVERFLOW IN THE
C            BUILTIN EXP-FUNCTION
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           POLYNOMIAL APPROXIMATIONS GIVEN BY E.E. ALLEN ARE USED FOR
C           CALCULATION.
C           FOR REFERENCE SEE
C           M. ABRAMOWITZ AND I.A. STEGUN,'HANDBOOK OF MATHEMATICAL
C           FUNCTIONS', U.S. DEPARTMENT OF COMMERCE, NATIONAL BUREAU OF
C           STANDARDS APPLIED MATHEMATICS SERIES, 1966, P.378.
C
C     ..................................................................
C
      SUBROUTINE I0(X,RI0)
      RI0=ABS(X)
      IF(RI0-3.75)1,1,2
    1 Z=X*X*7.111111E-2
      RI0=((((( 4.5813E-3*Z+3.60768E-2)*Z+2.659732E-1)*Z+1.206749E0)*Z
     1+3.089942E0)*Z+3.515623E0)*Z+1.
      RETURN
    2 Z=3.75/RI0
      RI0= EXP(RI0)/SQRT(RI0)*((((((((3.92377E-3*Z-1.647633E-2)*Z
     1+2.635537E-2)*Z-2.057706E-2)*Z+9.16281E-3)*Z-1.57565E-3)*Z
     2+2.25319E-3)*Z+1.328592E-2)*Z+3.989423E-1)
      RETURN
      END