Numerical modeling of primary breakup of a liquid jet under the influence of strong aerodynamic and turbulence effects is very challenging, especially for high liquid/gas,density ratio 0 (1000). A robust algorithm for large eddy simulation (LES) of two-phase flows is presented here. A coupled level set and volume of fluid (CLSVOF) technique is applied as the interface-tracking method in order to combine the advantages of level set and volume of fluid methods. The governing equations are discretized by introducing an extrapolated liquid velocity to minimize the interface momentum error. Since experimental studies on breakup of a single liquid drop in uniform gas flow are well documented, this test case is first used to validate the developed two-phase flow LES method. It is shown that the predicted drop breakup agrees quantitatively well with experiments for different Weber numbers. The solver is then applied to simulate primary breakup of liquid jets, which are more relevant to industrial applications. By simulating single round water jet atomization in high-speed coaxial airflow, it is found that the predicted liquid core breakup lengths at different air/liquid velocities agree closely with measured data, but only when appropriate turbulent inflow conditions are specified.