      SUBROUTINE CERF (MO, A, W)
C     ******************************************************************
C     FORTRAN SUBROUTINE FOR COMPLEX ERROR FUNCTION
C     ******************************************************************
C     MO = MODE OF OPERATION
C     A  = ARGUMENT (COMPLEX NUMBER)
C     W  = FUNCTION (COMPLEX NUMBER)
C     -------------------
      COMPLEX A, W
      DIMENSION AZ(2), QZ(2), SZ(2), EF(2)
      DIMENSION TS(2), SM(2), TM(2), QF(2), CD(18), CE(18)
      DATA CD(1) /0.00000000000000E00/,  CD(2) /2.08605856013476E-2/,
     1     CD(3) /8.29806940495687E-2/,  CD(4) /1.85421653326079E-1/,
     2     CD(5) /3.27963479382361E-1/,  CD(6) /5.12675279912828E-1/,
     3     CD(7) /7.45412958045105E-1/,  CD(8) /1.03695067418297E00/,
     4     CD(9) /1.40378061255437E00/,  CD(10)/1.86891662214001E00/,
     5     CD(11)/2.46314830523929E00/,  CD(12)/3.22719383737352E00/,
     6     CD(13)/4.21534348280013E00/,  CD(14)/5.50178873151549E00/,
     7     CD(15)/7.19258966683102E00/,  CD(16)/9.45170208076408E00/,
     8     CD(17)/1.25710718314784E+1/,  CD(18)/1.72483537216334E+1/
      DATA CE(1) /8.15723083324096E-2/,  CE(2) /1.59285285253437E-1/,
     1     CE(3) /1.48581625614499E-1/,  CE(4) /1.33219670836245E-1/,
     2     CE(5) /1.15690392878957E-1/,  CE(6) /9.78580959447535E-2/,
     3     CE(7) /8.05908834297624E-2/,  CE(8) /6.40204538609872E-2/,
     4     CE(9) /4.81445242767885E-2/,  CE(10)/3.33540658473295E-2/,
     5     CE(11)/2.05548099470193E-2/,  CE(12)/1.07847403887506E-2/,
     6     CE(13)/4.55634892214219E-3/,  CE(14)/1.43984458138925E-3/,
     7     CE(15)/3.07056139834171E-4/,  CE(16)/3.78156541168541E-5/,
     8     CE(17)/2.05173509616121E-6/,  CE(18)/2.63564823682747E-8/
C     -------------------
      AZ(1)=REAL(A)
      AZ(2)=AIMAG(A)
  001 ZS=AZ(1)*AZ(1)+AZ(2)*AZ(2)
      SZ(1)=AZ(1)*AZ(1)-AZ(2)*AZ(2)
      SZ(2)=2.0*AZ(1)*AZ(2)
      QZ(1)=+AZ(1)
      QZ(2)=+AZ(2)
      SN=+1.0
      IF(AZ(1))002,003,003
  002 QZ(1)=-AZ(1)
      QZ(2)=-AZ(2)
      SN=-1.0
  003 IF(ZS-1.0)014,014,004
  004 IF(ZS-38.0)005,006,006
  005 IF(SZ(1)+0.064*SZ(2)*SZ(2))014,014,012
C
C              ASYMPTOTIC EXPANSION
C
  006 SM(1)=0.0
      SM(2)=0.0
      QF(1)=+SZ(1)/(ZS*ZS)
      QF(2)=-SZ(2)/(ZS*ZS)
      QM=0.564189583547756*EXP(-SZ(1))
      TS(1)=QM*COS(-SZ(2))
      TS(2)=QM*SIN(-SZ(2))
      TM(1)=+(TS(1)*QZ(1)+TS(2)*QZ(2))/ZS
      TM(2)=-(TS(1)*QZ(2)-TS(2)*QZ(1))/ZS
      PM=-0.5
      GO TO 008
  007 PM=PM+1.0
      TS(1)=TM(1)*QF(1)-TM(2)*QF(2)
      TS(2)=TM(1)*QF(2)+TM(2)*QF(1)
      TM(1)=-PM*TS(1)
      TM(2)=-PM*TS(2)
      IF(ABS(SM(1))+ABS(TM(1)).NE.ABS(SM(1)))GO TO 008
      IF(ABS(SM(2))+ABS(TM(2)).EQ.ABS(SM(2)))GO TO 009
  008 SM(1)=SM(1)+TM(1)
      SM(2)=SM(2)+TM(2)
      IF(PM.LT.36.5)GO TO 007
  009 IF (QZ(1).LT.0.01) GO TO 019
      EF(1)=SM(1)
      EF(2)=SM(2)
  010 IF(MO.EQ.0)GO TO 011
      W=CMPLX(EF(1),EF(2))
      IF (SN.EQ.1.0) RETURN
      W=CMPLX(2.0-EF(1),-EF(2))
      RETURN
  011 EF(1)=SN*(1.0-EF(1))
      EF(2)=-SN*EF(2)
      W=CMPLX(EF(1),EF(2))
      RETURN
C
C         RATIONAL FUNCTION APPROXIMATION
C
  012 SM(1)=0.0
      SM(2)=0.0
      QM=0.564189583547756*EXP(-SZ(1))
      TS(1)=QM*COS(-SZ(2))
      TS(2)=QM*SIN(-SZ(2))
      QF(1)=TS(1)*QZ(1)-TS(2)*QZ(2)
      QF(2)=TS(1)*QZ(2)+TS(2)*QZ(1)
      DO 013 I=1,18
      TS(1)=SZ(1)+CD(I)
      TS(2)=SZ(2)
      SS=TS(1)*TS(1)+TS(2)*TS(2)
      TM(1)=+CE(I)*TS(1)/SS
      TM(2)=-CE(I)*TS(2)/SS
      SM(1)=SM(1)+TM(1)
      SM(2)=SM(2)+TM(2)
  013 CONTINUE
      EF(1)=QF(1)*SM(1)-QF(2)*SM(2)
      EF(2)=QF(1)*SM(2)+QF(2)*SM(1)
      GO TO 010
C
C            TAYLOR SERIES EXPANSION
C
  014 TM(1)=1.128379167095513*QZ(1)
      TM(2)=1.128379167095513*QZ(2)
      SM(1)=TM(1)
      SM(2)=TM(2)
      PM=0.0
  015 PM=PM+1.0
      DM=2.0*PM+1.0
      TS(1)=TM(1)*SZ(1)-TM(2)*SZ(2)
      TS(2)=TM(1)*SZ(2)+TM(2)*SZ(1)
      TM(1)=-TS(1)/PM
      TM(2)=-TS(2)/PM
      TS(1)=TM(1)/DM
      TS(2)=TM(2)/DM
      IF(ABS(SM(1))+ABS(TS(1)).NE.ABS(SM(1)))GO TO 016
      IF(ABS(SM(2))+ABS(TS(2)).EQ.ABS(SM(2)))GO TO 017
  016 SM(1)=SM(1)+TS(1)
      SM(2)=SM(2)+TS(2)
      GO TO 015
  017 IF(MO.NE.0)GO TO 018
      W=CMPLX(SN*SM(1),SN*SM(2))
      RETURN
  018 EF(1)=1.0-SN*SM(1)
      EF(2)=-SN*SM(2)
      W=CMPLX(EF(1),EF(2))
      RETURN
C
C         MODIFIED ASYMPTOTIC EXPANSION
C
  019 SN=-SN
      GO TO 017
      END
