Summary: A stochastic subgrid model for large-eddy simulation of atomizing spray is developed. Following Kolmogorov’s concept of viewing solid particle-breakup as a discrete random process, atomization of liquid blobs at high relative liquid-to-gas velocity is considered in the framework of uncorrelated breakup events, independent of the initial droplet size. Kolmogorov’s discrete model of breakup is rewritten in the form of differential Fokker-Planck equation for the PDF of droplet radii. Along with the Lagrangian tracking of spray dynamics, the size and number density of the newly produced droplets is governed by the evolution of this PDF in the space of droplet-radius. The parameters of the model are obtained dynamically by relating them to the local Weber number with two-way coupling between the gas and liquid phases. Computations of spray are performed for the representative conditions encountered in idealized diesel and gas-turbine engine configurations. A broad spectrum of droplet sizes is obtained at each location with co-existence of large and small droplets. A novel numerical algorithm capable of simultaneously simulating individual droplets as well as a group of droplets with similar properties commonly known as parcels is proposed and compared with standard parcels-approach usually employed in the computations of multiphase flows. The present approach is shown to be computationally efficient and captures the complex fragmentary process of liquid atomization.