diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp
index 103811dca..fb501ac38 100644
--- a/opm/autodiff/BlackoilModelBase_impl.hpp
+++ b/opm/autodiff/BlackoilModelBase_impl.hpp
@@ -2384,7 +2384,7 @@ namespace detail {
 
         const int nc = Opm::AutoDiffGrid::numCells(grid_);
         const int nw = localWellsActive() ? wells().number_of_wells : 0;
-        const int np = fluid_.numPhases();
+        const int np = numPhases();
         assert(int(rq_.size()) == np);
 
         const V pv = geo_.poreVolume();
@@ -2438,40 +2438,52 @@ namespace detail {
         // Residual in Pascal can have high values and still be ok.
         const double maxWellResidualAllowed = 1000.0 * maxResidualAllowed();
 
-        // if one of the residuals is NaN, throw exception, so that the solver can be restarted
-        if ( std::isnan(mass_balance_residual[Water]) || mass_balance_residual[Water] > maxResidualAllowed() ||
-            std::isnan(mass_balance_residual[Oil])   || mass_balance_residual[Oil]   > maxResidualAllowed() ||
-            std::isnan(mass_balance_residual[Gas])   || mass_balance_residual[Gas]   > maxResidualAllowed() ||
-            std::isnan(CNV[Water]) || CNV[Water] > maxResidualAllowed() ||
-            std::isnan(CNV[Oil]) || CNV[Oil] > maxResidualAllowed() ||
-            std::isnan(CNV[Gas]) || CNV[Gas] > maxResidualAllowed() ||
-            std::isnan(well_flux_residual[Water]) || well_flux_residual[Water] > maxResidualAllowed() ||
-            std::isnan(well_flux_residual[Oil]) || well_flux_residual[Oil] > maxResidualAllowed() ||
-            std::isnan(well_flux_residual[Gas]) || well_flux_residual[Gas] > maxResidualAllowed() ||
-            std::isnan(residualWell) || residualWell > maxWellResidualAllowed )
-        {
-            OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN or too large!");
+        for (int idx = 0; idx < np; ++idx) {
+            if (std::isnan(mass_balance_residual[idx])
+                || std::isnan(CNV[idx])
+                || std::isnan(well_flux_residual[idx])) {
+                OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName(idx));
+            }
+            if (mass_balance_residual[idx] > maxResidualAllowed()
+                || CNV[idx] > maxResidualAllowed()
+                || well_flux_residual[idx] > maxResidualAllowed()) {
+                OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName(idx));
+            }
+        }
+        if (std::isnan(residualWell) || residualWell > maxWellResidualAllowed) {
+            OPM_THROW(Opm::NumericalProblem, "NaN or too large residual for well control equation");
         }
 
         if ( terminal_output_ )
         {
             // Only rank 0 does print to std::cout
             if (iteration == 0) {
-                std::cout << "\nIter  MB(WATER)   MB(OIL)    MB(GAS)       CNVW       CNVO       CNVG   W-FLUX(W)  W-FLUX(O)  W-FLUX(G)\n";
+                // std::cout << "\nIter  MB(WATER)   MB(OIL)    MB(GAS)       CNVW       CNVO       CNVG   W-FLUX(W)  W-FLUX(O)  W-FLUX(G)\n";
+                std::cout << "\nIter";
+                for (int idx = 0; idx < np; ++idx) {
+                    std::cout << "   MB(" << phaseName(idx).substr(0, 3) << ") ";
+                }
+                for (int idx = 0; idx < np; ++idx) {
+                    std::cout << "    CNV(" << phaseName(idx).substr(0, 1) << ") ";
+                }
+                for (int idx = 0; idx < np; ++idx) {
+                    std::cout << "  W-FLUX(" << phaseName(idx).substr(0, 1) << ")";
+                }
+                std::cout << '\n';
             }
             const std::streamsize oprec = std::cout.precision(3);
             const std::ios::fmtflags oflags = std::cout.setf(std::ios::scientific);
-            std::cout << std::setw(4) << iteration
-                      << std::setw(11) << mass_balance_residual[Water]
-                      << std::setw(11) << mass_balance_residual[Oil]
-                      << std::setw(11) << mass_balance_residual[Gas]
-                      << std::setw(11) << CNV[Water]
-                      << std::setw(11) << CNV[Oil]
-                      << std::setw(11) << CNV[Gas]
-                      << std::setw(11) << well_flux_residual[Water]
-                      << std::setw(11) << well_flux_residual[Oil]
-                      << std::setw(11) << well_flux_residual[Gas]
-                      << std::endl;
+            std::cout << std::setw(4) << iteration;
+            for (int idx = 0; idx < np; ++idx) {
+                std::cout << std::setw(11) << mass_balance_residual[idx];
+            }
+            for (int idx = 0; idx < np; ++idx) {
+                std::cout << std::setw(11) << CNV[idx];
+            }
+            for (int idx = 0; idx < np; ++idx) {
+                std::cout << std::setw(11) << well_flux_residual[idx];
+            }
+            std::cout << std::endl;
             std::cout.precision(oprec);
             std::cout.flags(oflags);
         }
@@ -2490,7 +2502,7 @@ namespace detail {
 
         const int nc = Opm::AutoDiffGrid::numCells(grid_);
         const int nw = localWellsActive() ? wells().number_of_wells : 0;
-        const int np = fluid_.numPhases();
+        const int np = numPhases();
 
         const V pv = geo_.poreVolume();
         std::vector<double> R_sum(np);