An accurate and robust numerical method has been developed to simulate turbulent two-phase flows. The phase interface is tracked by the level-set method to capture frequent topological changes due to breaking or merging. Because of the broad-band characteristics of length scales in two-phase flow, a Lagrangian drop breakup model has been developed, which is coupled to the level-set method. In this approach, small subgrid droplets produced from resolved ligaments are then transferred from the level-set representation to the Lagrangian particles. The further secondary atomization is handled by a stochastic breakup model. When pinching-off of ligaments is not resolved on the level-set grid, a capillary breakup model is used to predict the drop size distribution from the pinching off and inserted as Lagrangian drops. This method improves the mass conservation as well as reducing the computational cost.