If you use the p^(n+1)-p^n formulation of either Chorin or IPCS, you can apply the pressure condition (although it is not natural in the sense of boundary conditions).
In your code you are applying the pressure boundary condition to the \phi. This means that
p^n+1 = p^n + bcp on any boundary with a dirichletbc, which is not correct.
What you can do is to solve for \phi, with zero Dirichlet conditions on this boundary.
Then after calling
you could assign bcp to p_.vector using dolfinx.fem.petsc.set_bc(p_.vector, bcp).