      FUNCTION ERF(X)
C     ******************************************************************
C     EVALUATION OF THE REAL ERROR FUNCTION
C     ******************************************************************
      DIMENSION A(4),B(4),P(8),Q(8),R(5),S(5)
      DATA A(1)/-1.65581836870402E-4/, A(2)/3.25324098357738E-2/,
     1     A(3)/1.02201136918406E-1/,  A(4)/1.12837916709552E00/
      DATA B(1)/4.64988945913179E-3/,  B(2)/7.01333417158511E-2/,
     1     B(3)/4.23906732683201E-1/,  B(4)/1.00000000000000E00/
      DATA P(1)/-1.36864857382717E-7/, P(2)/5.64195517478974E-1/,
     1     P(3)/7.21175825088309E00/,  P(4)/4.31622272220567E01/,
     2     P(5)/1.52989285046940E02/,  P(6)/3.39320816734344E02/,
     3     P(7)/4.51918953711873E02/,  P(8)/3.00459261020162E02/
      DATA Q(1)/1.00000000000000E00/,  Q(2)/1.27827273196294E01/,
     1     Q(3)/7.70001529352295E01/,  Q(4)/2.77585444743988E02/,
     2     Q(5)/6.38980264465631E02/,  Q(6)/9.31354094850610E02/,
     3     Q(7)/7.90950925327898E02/,  Q(8)/3.00459260956983E02/
      DATA R(1)/2.10144126479064E00/,  R(2)/2.62370141675169E01/,
     1     R(3)/2.13688200555087E01/,  R(4)/4.65807828718470E00/,
     2     R(5)/2.82094791773523E-1/
      DATA S(1)/9.41537750555460E01/,  S(2)/1.87114811799590E02/,
     1     S(3)/9.90191814623914E01/,  S(4)/1.80124575948747E01/,
     2     S(5)/1.00000000000000E00/
      DATA C/5.64189583547756E-1/
C     -------------------
      AX=ABS(X)
      X2=AX*AX
      IF (AX.GE.0.5) GO TO 10
      TOP=((A(1)*X2+A(2))*X2+A(3))*X2+A(4)
      BOT=((B(1)*X2+B(2))*X2+B(3))*X2+B(4)
      ERF=X*TOP/BOT
      RETURN
C
   10 IF (AX.GT.4.0) GO TO 20
      TOP=((((((P(1)*AX+P(2))*AX+P(3))*AX+P(4))*AX+P(5))*AX
     *                 +P(6))*AX+P(7))*AX+P(8)
      BOT=((((((Q(1)*AX+Q(2))*AX+Q(3))*AX+Q(4))*AX+Q(5))*AX
     *                 +Q(6))*AX+Q(7))*AX+Q(8)
      ERF=1.0-EXP(-X2)*TOP/BOT
      IF (X.LT.0.0) ERF=-ERF
      RETURN
C
   20 ERF=1.0
      IF (AX.GE.5.54) GO TO 21
      T=1.0/X2
      TOP=(((R(1)*T+R(2))*T+R(3))*T+R(4))*T+R(5)
      BOT=(((S(1)*T+S(2))*T+S(3))*T+S(4))*T+S(5)
      ERF=C-TOP/(X2*BOT)
      ERF=1.0-EXP(-X2)*ERF/AX
   21 IF (X.LT.0.0) ERF=-ERF
      RETURN
      END
