In this work, we study underlay radar-massive MIMO cellular coexistence in LoS/near-LoS channels, where both systems have 3D beamforming capabilities. Using mathematical tools from stochastic geometry, we derive an upper bound on the average interference power at the radar due to the 3D massive MIMO cellular downlink under the worst-case `cell-edge beamforming' conditions. To overcome the technical challenges imposed by asymmetric and arbitrarily large cells, we devise a novel construction in which each Poisson Voronoi (PV) cell is bounded by its circumcircle to bound the effect of the random cell shapes on average interference. Since this model is intractable for further analysis due to the correlation between adjacent PV cells' shapes and sizes, we propose a tractable nominal interference model, where we model each PV cell as a circular disk with an area equal to the average area of the typical cell. We quantify the gap in the average interference power between these two models and show that the upper bound is tight for realistic deployment parameters. We also compare them with a more practical but intractable MU-MIMO scheduling model to show that our worst-case interference models show the same trends and do not deviate significantly from realistic scheduler models. Under the nominal interference model, we characterize the interference distribution using the dominant interferer approximation by deriving the equi-interference contour expression when the typical receiver uses 3D beamforming. Finally, we use tractable expressions for the interference distribution to characterize radar's spatial probability of false alarm/detection in a quasi-static target tracking scenario. Our results reveal useful trends in the average interference as a function of the deployment parameters (BS density, exclusion zone radius, antenna height, transmit power of each BS, etc.).