This article presents a massively parallel and robust strategy to perform the simulation of turbulent incompressible two-phase flows on unstructured grids in complex geometries. This strategy relies on a combination of a narrow-band accurate conservative level set (ACLS)/ghost-fluid framework with isotropic adaptive mesh refinement. This combination enables to accurately capture interface dynamics and topology, and the small physical scales at the liquid-gas interface are resolved at an affordable cost. The ACLS method, even if not strictly mass-conserving, ensures the exact conservation of a smoothed phase indicator, which minimizes the liquid mass conservation errors. In the accurate conservative level set framework, presented first in Desjardins et al. (2008) [25], the interface is defined as the iso-contour of a hyperbolic tangent function, which is advected by the fluid, and then reshaped using a reinitialization equation. Several forms of this reinitialization exist: the original ACLS form proposed by Desjardins et al. involves numerical estimation of the hyperbolic tangent gradient, which is difficult to compute accurately on unstructured meshes. It is thus susceptible to induce artificial deformation of the interface. A new form has been recently proposed in Chiodi and Desjardins (2017) [29], which takes advantage of a mapping onto a classical distance level set while much better preserving the interface shape. Nevertheless, the implementation of this new form on unstructured grids requires special care. In this work, a robust implementation of this new form on unstructured meshes is proposed and implemented in the YALES2 low-Mach flow solver. In order to compute interface normals and curvature, the signed-distance function is reconstructed in parallel at nodes in the narrow band around the interface using a Geometric-Projection Marker Method (GPMM). This method relies on the triangulation of the level set iso-contour and exact geometric projection to the closest surface elements. Spatial convergence, robustness and efficiency of the overall procedure are firstly demonstrated through classical interface transport test cases and two-phase flow examples. Eventually, to emphasize the significant computational gain using adaptive mesh refinement and the ability to compute complex turbulent flows with large density ratios, two Large-Eddy Simulations (LES) of atomizing liquid jets in air are presented, each one at various resolutions. The first one is a low-pressure water jet in quiescent air from a compound nozzle with full computation of the internal injector flow, while the latter is a high-pressure kerosene jet in crossflow. Both simulations are validated against experiments, demonstrating the potential of the method to access a deep numerical insight into jet instabilities and internal flow dynamics with 3D unstructured meshes.