The present work focuses on applying the Dynamically Orthogonal Primitive Equations (DO-PE) for realistic high-resolution stochastic ocean forecasting in regions with complex ocean dynamics. In the first part, we identify and test a streamlined process to create multi-region initial conditions for the DO-PE framework, starting from temporally and spatially sparse historical data. The process presented allows us to start from a relatively small but relevant set of measured temperature and salinity historical vertical profiles (on the order of hundreds) and to generate a massive set of initial conditions (on the order of millions) in a stochastic subspace, while still ensuring that the initial statistics respect the physical processes, modeled complex dynamics, and uncertain initial conditions of the examined domain. To illustrate the methodology, two practical examples-one in the Gulf of Mexico and another in the Alboran Sea--are provided, along with a review of the ocean dynamics for each region. In the second part, we present a case study of three massive stochastic DO-PE forecasts, corresponding to ensembles of one million members, in the Gulf of Mexico region. We examine the effect of adding more dynamic DO modes (i.e., stochastic dimensions) and show that it tends to statistical convergence along with an enhancement of the uncertainty captured by the DO forecast realizations, both by increasing the variance of already existing features as well as by adding new uncertain features. We also use this case study to validate the DO-PE methodology for realistic high-resolution probabilistic ocean forecasting. We show good accuracy against equivalent deterministic simulations, starting from the same initial conditions and simulated with the same assumptions, setup, and original ocean model equations. Importantly, by comparing the reduced-order realizations against their deterministic counterparts, we show that the errors due to the DO subspace truncation are much smaller and growing slower than the fields themselves are evolving in time, both in the Root Mean Square Error (RMSE) sense as well as in the 3D multivariate ocean field sense. Based on these observations, we conclude that the DO-PE realizations closely match their full-order equivalents, thus enabling massive forecast ensembles with practically low numerical errors at a tractable computational cost.