The dynamics of Oceanic flows change drastically with scale. While mesoscales (10-100km) are well-understood and containing 90% of the ocean's kinetic energy and characterized by geostrophically balanced flows, the submesoscale (<10km) is very different. This regime is characterized by unbalanced flows that create fine-scale structures requiring high resolution simulations, hence computationally challenging. This has in general motivated the development of new parametrizations and data-driven surrogate modeling. In this work, we develop a new data-driven stochastic parametrization for submesoscale dynamics using an unsupervised conditional diffusion model. We focus on a two-vertical-mode ocean model, in which the evolution of the barotropic (vertically averaged) vorticity is coupled to baroclinic instabilities governed by multiple nonlinear equations. Instead of explicitly evolving the baroclinic modes, we replace their influence with an abstract field in a new surrogate barotropic vorticity equation. The abstract field modifies the surrogate flow such that the turbulent trajectory resembles the turbulence coming from the original two-mode-model equations. This term is sampled conditionally from the resolved barotropic state, without access to or supervision from baroclinic variables. The resulting learning problem is inherently unsupervised, as the target stochastic forcing is not observable and hence cannot be storable as target data. To ensure numerical stability and physical realism, we incorporate control mechanisms that regulate the injected variance and enforce consistency with key flow statistics. We evaluate the surrogate model by comparing it with the full two-mode system using energy spectra, fluxes, and probability distributions. The surrogate model generalizes across forcings, initial conditions, and flow regimes. The surrogate model significantly accelerates simulations while retaining essential submesoscale characteristics.