P_integral Subroutine

  SUBROUTINE P_integral(z,P)
    REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: P
    REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: z
    REAL(rp) :: a
    INTEGER :: ll,ss

    ss = SIZE(z)

    P = 0.0_rp

    do ll=1_idef,ss
       IF (z(ll) .LT. 0.5_rp) THEN
          a = (2.16_rp/2.0_rp**(2.0_rp/3.0_rp))*z(ll)**(1.0_rp/3.0_rp)
          P(ll) = IntBesselK(z(ll),a) + IntK(5.0_rp/3.0_rp,a)
       ELSE IF ((z(ll) .GE. 0.5_rp).AND.(z(ll) .LT. 2.5_rp)) THEN
          a = 0.72_rp*(z(ll) + 1.0_rp)
          P(ll) = IntBesselK(z(ll),a) + IntK(5.0_rp/3.0_rp,a)
       ELSE
          P(ll) = IntK(5.0_rp/3.0_rp,z(ll))
       END IF
    end do
  END SUBROUTINE P_integral