PO and PN have been newly identified in several oxygen-rich circumstellar envelopes, using the Submillimeter Telescope of the Arizona Radio Observatory. The J=5→4 and J=6→5 transitions of PN near 235 and 282 GHz, and the lambda doublets originating in the J=5.5→4.5 and J=6.5→5.5 lines of PO at 240 and 284 GHz, have been detected toward the shells of asymptotic giant branch (AGB) stars TX Cam and R Cas. A similar set of lines has been observed toward the supergiant NML Cyg, and new transitions of these two molecules were also measured toward the AGB star IK Tau. Along with the previous data from VY Canis Majoris (VY CMa), these spectral lines were analyzed using the non-local thermodynamic equilibrium (non-LTE) circumstellar modeling code, ESCAPADE. For the AGB stars, peak abundances found for PN and PO were f∼(1-2)×10 −8 and (0.5-1)×10−7 , respectively, while those for the supergiants were f(PN)∼(0.3-0.7)×10 −8 and f(PO)∼(5-7)×10 −8 . PN was well modeled with a spherical radial distribution, suggesting formation near the stellar photosphere, perhaps enhanced by shocks. PO was best reproduced by a shell model, indicating a photochemical origin, except for VY CMa. Overall, the abundance of PO is a factor of 5-20 greater than that of PN. This study suggests that phosphorus-bearing molecules are common in O-rich envelopes, and that a significant amount of phosphorus (>20%) remains in the gas phase.