FLICA-4: a three-dimensional two-phase flow computer code with advanced numerical methods for nuclear applications. This paper is devoted to new numerical methods developed for three-dimensional two-phase flow calculations. These methods are finite volume numerical methods. They are based on an extension of Roe’s approximate Riemann solver to define convective fluxes versus mean cell quantities [Godunov, S.K., 1959, Math. Sb. 47, 217; Roe, P.L., 1981, Approximate Riemanns solvers parameter vectors and difference scheme. J. Comp. Phys. 43, 357–372; Toumi, I., 1992, A weak formulation of Roe’s approximate Riemann solver. J. Comp. Phys. 102, 360–373]. To go forward in time, a linearized conservative implicit integrating step is used [Yee, H.C., 1987. NASA TM-89464], together with a Newton iterative method. We also present here some improvements performed to obtain a fully implicit solution method that provides fast running steady state calculations. This kind of numerical method, which is used widely for fluid dynamic calculations, has proved to be very efficient for the numerical solution to two-phase flow problems. This numerical method has been implemented for the three-dimensional thermal-hydraulic code FLICA-4 that is mainly dedicated to core thermal-hydraulic transient and steady-state analysis [Toumi, I., Caruge, D., 1998. An implicit second order method for 3D two phase flow calculations. Nucl. Sci. Eng. 130, 213–225; Raymond, P., Toumi, I., 1992. Numerical method for three-dimensional steady-state two-phase flow calculation, NURETH-5, Salt Lake City]. Hereafter, elements of physical validation against hydraulic and two-phase flow rod bundle experiments are presented. We will also find some results obtained for the EPR reactor running in a steady-state at 60% of nominal power with three pumps out of four, and a thermal-hydraulic core analysis for a 1300 MW PWR at low flow Steam-Line-Break conditions.