To satisfy the robust and real-time requirements of power flow calculation for large-scale power systems, a globally convergent method is proposed with trust-region techniques, which shows satisfying robustness and convergence. Then, this method is combined with Newton's method to achieve a two-stage algorithm, benefiting from their different advantages. In the first stage, the proposed globally convergent method is used for searching power flow solution. When the values of state variables in a certain iteration are close enough to the real operational point, the algorithm enters the second stage to use Newton's method to achieve the solution. This two-stage algorithm can achieve an accurate solution for solvable cases, and can also achieve a least-square solution, which is an approximate solution for unsolvable cases. Numerical experiments show that the proposed globally convergent method and two-stage algorithm have better robustness and efficiency than the existing methods in previous research. They have universality for well- and ill-conditioned systems as well as the cases under ill-conditioned operational modes, heavy loads and inappropriate initial values. They can also handle different limit violations, such as reactive power limit in PV buses and voltage limit with tap changers action, which could significantly benefit real practice.