fermi Function

  elemental function fermi(x,beta,limit)
    real(8),intent(in)          :: x, beta
    real(8),optional,intent(in) :: limit
    real(8)                     :: fermi,arg,limit_
    limit_ = 200d0 ; if(present(limit))limit_=abs(limit)
    arg = x*beta
    if(arg < -limit_)then
       fermi = 1d0
    elseif(arg > limit_)then
       fermi = 0d0
    else
       fermi = 1d0/(1d0+exp(arg))
    endif
  end function fermi