itsl0 Subroutine

   subroutine itsl0(x,Tl0)

      real(wp),intent(in) :: x  !! Upper limit  ( x ≥ 0 )
      real(wp),intent(out) :: Tl0 !! Integration of L0(t) from 0 to x

      real(wp) :: a(18) , a0 , a1 , af , r , rd , s , s0 , ti
      integer :: k

      real(wp),parameter :: two_pi = 2.0_wp * pi

      r = 1.0_wp
      if ( x<=20.0_wp ) then
         s = 0.5_wp
         do k = 1 , 100
            rd = 1.0_wp
            if ( k==1 ) rd = 0.5_wp
            r = r*rd*k/(k+1.0_wp)*(x/(2.0_wp*k+1.0_wp))**2
            s = s + r
            if ( abs(r/s)<1.0e-12_wp ) exit
         enddo
         Tl0 = 2.0_wp/pi*x*x*s
      else
         s = 1.0_wp
         do k = 1 , 10
            r = r*k/(k+1.0_wp)*((2.0_wp*k+1.0_wp)/x)**2
            s = s + r
            if ( abs(r/s)<1.0e-12_wp ) exit
         enddo
         s0 = -s/(pi*x*x) + 2.0_wp/pi*(log(2.0_wp*x)+gamma)
         a0 = 1.0_wp
         a1 = 5.0_wp/8.0_wp
         a(1) = a1
         do k = 1 , 10
            af = ((1.5_wp*(k+0.50_wp)*(k+5.0_wp/6.0_wp)*a1-0.5_wp*(k+0.5_wp) &
                 **2*(k-0.5_wp)*a0))/(k+1.0_wp)
            a(k+1) = af
            a0 = a1
            a1 = af
         enddo
         ti = 1.0_wp
         r = 1.0_wp
         do k = 1 , 11
            r = r/x
            ti = ti + a(k)*r
         enddo
         Tl0 = ti/sqrt(two_pi*x)*exp(x) + s0
      endif

   end subroutine itsl0