function  ze1(z)
  !
  !       ====================================================
  !       Purpose: Compute complex exponential integral E1(z)
  !       Input :  z   --- Argument of E1(z)
  !       Output:  CE1 --- E1(z)
  !       ====================================================
  !
  ! Déclaration des variables en passage de paramètre
  !
  implicit none
  complex(8), intent(in) :: z
  complex(8) :: ze1
  !
  !  Déclaration des variables locales
  !
  integer(4) :: k
  real(8) :: pi,el,x,a0
  complex(8) :: cr,ct0,ct
  parameter(pi=3.141592653589793D0,el=0.5772156649015328D0)
  !
  ! traitement en fonction des différents cas
  ! - Z nul entraîne E1 infini
  ! - module de Z inférieur à 10 ou 20 : formule log+gam+somme
  ! - module de Z supérieur à 10 ou 20 : formule asymptotique
  !
   x=dreal(z)
   a0=cdabs(z)
   if (a0==0.0D0) then
     ze1 = dcmplx(1.0D+300,0.0D0)
   else if ((a0 <= 10.D0).or.(x <= 0.D0.and.a0 <= 20.D0)) then
     ze1 = dcmplx(1.0D0,0.0D0)
     cr = dcmplx(1.0D0,0.0D0)
     k=0
     do while (cdabs(cr)<=cdabs(ze1)*1.0D-15)
       k = k+1
       cr = -cr*k*z/(k+1.0D0)**2
       ze1 = ze1+cr
     end do
     ze1 = -el-cdlog(z)+z*ze1
   else
     ct0 = dcmplx(0.0D0,0.0D0)
     do k=120,1,-1
        ct0 = k/(1.0D0+k/(z+ct0))
     end do
     ct = 1.0D0/(z+ct0)
     ze1 = cdexp(-z)*ct
     if (x <= 0.D0.AND.dimag(z) == 0.0) ze1 = ze1-pi*dcmplx(0.D0,1.D0)
   end if
  !
   return
  end function