This paper describes the implementation of a 3D parallel and Cartesian level set (LS) method coupled with a volume of fluid (VOF) method into the commercial CFD code FLUENT for modeling the gas-liquid interface in bubbly flow. Both level set and volume of fluid methods belong to the so called "one" fluid methods, where a single set of conservation equations is solved and the interface is captured via a scalar function. Since both LS and VOF have advantages and disadvantages, our aim is to couple these two methods to obtain a method, which is superior to both standalone LS and VOF and verify it versus a selection of test cases. VOF is already available in FLUENT, so we implemented an LS method into FLUENT via user defined functions. The level set function is used to compute the surface tension contribution to the momentum equations, via curvature and its normal to the interface, using the Brackbill method while the volume of fluid function is used to capture the interface itself. A re-initialization equation is implemented and solved at every time step using a fifth-order weighted essentially nonoscillatory scheme for the spatial derivative, and a first-order Euler method for time integration. The coupling effect is introduced by solving at the end of each time step an equation, which connects the volume fractions with the level set function. The verification of parasitic currents and interfacial deformation due to numerical error is assessed in comparison to original VOF scheme. Validation is presented for free rising bubbles of different diameters for Morton numbers ranging from 10(2) to 10(-11). [DOI: 10.1115/1.4002166]