psi_spec Subroutine

   subroutine psi_spec(x,Ps)

      real(wp),intent(in) :: x  !! Argument of `psi(x)`
      real(wp),intent(out) :: Ps !! `psi(x)`

      real(wp) :: a1 , a2 , a3 , a4 , a5 , a6 , a7 , a8 , &
                  s , x2 , xa
      integer k , n

      xa = abs(x)
      s = 0.0_wp
      if ( x==int(x) .and. x<=0.0_wp ) then
         Ps = 1.0e+300_wp
         return
      elseif ( xa==int(xa) ) then
         n = int(xa)
         do k = 1 , n - 1
            s = s + 1.0_wp/k
         enddo
         Ps = -gamma + s
      elseif ( xa+0.5_wp==int(xa+0.5_wp) ) then
         n = int(xa - 0.5_wp)
         do k = 1 , n
            s = s + 1.0/(2.0_wp*k-1.0_wp)
         enddo
         Ps = -gamma + 2.0_wp*s - 1.386294361119891_wp
      else
         if ( xa<10.0_wp ) then
            n = 10 - int(xa)
            do k = 0 , n - 1
               s = s + 1.0_wp/(xa+k)
            enddo
            xa = xa + n
         endif
         x2 = 1.0_wp/(xa*xa)
         a1 = -.8333333333333e-01_wp
         a2 = .83333333333333333e-02_wp
         a3 = -.39682539682539683e-02_wp
         a4 = .41666666666666667e-02_wp
         a5 = -.75757575757575758e-02_wp
         a6 = .21092796092796093e-01_wp
         a7 = -.83333333333333333e-01_wp
         a8 = .4432598039215686e0_wp
         Ps = log(xa) - 0.5_wp/xa +                                  &
              x2*(((((((a8*x2+a7)*x2+a6)*x2+a5)*x2+a4)*x2+a3)*x2+a2) &
              *x2+a1)
         Ps = Ps - s
      endif
      if ( x<0.0_wp ) Ps = Ps - pi*cos(pi*x)/sin(pi*x) - 1.0_wp/x

   end subroutine psi_spec