diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
index 4c8fed7b8..ca4e1ce5b 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
@@ -1900,79 +1900,60 @@ namespace {
 
         const std::vector<PhasePresence> cond = phaseCondition();
 
-        double CNVW = 0.;
-        double CNVO = 0.;
-        double CNVG = 0.;
+        double CNV[MaxNumPhases] = {0., 0., 0.};
+        double R_sum[MaxNumPhases] = {0., 0., 0.};
+        double B_avg[MaxNumPhases] = {0., 0., 0.};
+        double maxCoeff[MaxNumPhases] = {0., 0., 0.};
+        double mass_balance_residual[MaxNumPhases] = {0., 0., 0.};
+        Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> B;
+        Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> R;
+        Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> tempV;
 
-        double RW_sum = 0.;
-        double RO_sum = 0.;
-        double RG_sum = 0.;
-
-        double BW_avg = 0.;
-        double BO_avg = 0.;
-        double BG_avg = 0.;
-
-        if (active_[Water]) {
-            const int pos = pu.phase_pos[Water];
-            const ADB& tempBW = rq_[pos].b;
-            V BW = 1./tempBW.value();
-            V RW = residual_.material_balance_eq[Water].value();
-            BW_avg = BW.sum()/nc;
-            const V tempV = RW.abs()/pv;
-
-            CNVW = BW_avg * dt * tempV.maxCoeff();
-            RW_sum = RW.sum();
+        for(int idx=0; idx<MaxNumPhases; ++idx)
+        {
+            if (active_[idx]) {
+                const int pos = pu.phase_pos[idx];
+                const ADB& tempB = rq_[pos].b;
+                B.col(pu.phase_pos[idx])  = 1./tempB.value();
+                R.col(pu.phase_pos[idx]) = residual_.material_balance_eq[pu.phase_pos[idx]].value();
+                tempV.col(pu.phase_pos[idx]) = R.col(pu.phase_pos[idx]).abs()/pv;
+            }
+        }
+        // Do the global reductions
+        for(int idx=0; idx<MaxNumPhases; ++idx)
+        {
+            if (active_[idx]) {
+                B_avg[pu.phase_pos[idx]] = B.sum()/nc;
+                maxCoeff[pu.phase_pos[idx]]=tempV.col(pu.phase_pos[idx]).maxCoeff();
+                R_sum[pu.phase_pos[idx]] = R.col(pu.phase_pos[idx]).sum();
+            }else{
+                R_sum[pu.phase_pos[idx]] = B_avg[pu.phase_pos[idx]] = 0.;
+            }
         }
 
-        if (active_[Oil]) {
-            // Omit the disgas here. We should add it.
-            const int pos = pu.phase_pos[Oil];
-            const ADB& tempBO = rq_[pos].b;
-            V BO = 1./tempBO.value();
-            V RO = residual_.material_balance_eq[Oil].value();
-            BO_avg = BO.sum()/nc;
-            const V tempV = RO.abs()/pv;
-
-            CNVO = BO_avg * dt * tempV.maxCoeff();
-            RO_sum = RO.sum();
+        bool converged_MB = true;
+        bool converged_CNV = true;
+        // Finish computation
+        for(int idx=0; idx<MaxNumPhases; ++idx)
+        {
+            CNV[idx] = B_avg[pu.phase_pos[idx]] * dt * maxCoeff[pu.phase_pos[idx]];
+            mass_balance_residual[idx] = std::abs(B_avg[pu.phase_pos[idx]]*R_sum[pu.phase_pos[idx]]) * dt / pvSum;
+            converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
+            converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
         }
 
-        if (active_[Gas]) {
-            // Omit the vapoil here. We should add it.
-            const int pos = pu.phase_pos[Gas];
-            const ADB& tempBG = rq_[pos].b;
-            V BG = 1./tempBG.value();
-            V RG = residual_.material_balance_eq[Gas].value();
-            BG_avg = BG.sum()/nc;
-            const V tempV = RG.abs()/pv;
-
-            CNVG = BG_avg * dt * tempV.maxCoeff();
-            RG_sum = RG.sum();
-        }
-
-        const double mass_balance_residual_water = std::abs(BW_avg*RW_sum) * dt / pvSum;
-        const double mass_balance_residual_oil   = std::abs(BO_avg*RO_sum) * dt / pvSum;
-        const double mass_balance_residual_gas   = std::abs(BG_avg*RG_sum) * dt / pvSum;
-
-        bool converged_MB = (mass_balance_residual_water < tol_mb)
-                         && (mass_balance_residual_oil   < tol_mb)
-                         && (mass_balance_residual_gas   < tol_mb);
-
-        bool converged_CNV = (CNVW < tol_cnv) && (CNVO < tol_cnv) && (CNVG < tol_cnv);
-
-        const double residualWellFlux = infinityNorm( residual_.well_flux_eq );
-        const double residualWell     = infinityNorm( residual_.well_eq );
-
+        double residualWellFlux = infinityNorm(residual_.well_flux_eq);
+        double residualWell = infinityNorm(residual_.well_eq);
         bool converged_Well = (residualWellFlux < 1./Opm::unit::day) && (residualWell < Opm::unit::barsa);
         bool converged = converged_MB && converged_CNV && converged_Well;
 
         // 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(CNVW) || CNVW > maxResidualAllowed() ||
-            std::isnan(CNVO) || CNVO > maxResidualAllowed() ||
-            std::isnan(CNVG) || CNVG > maxResidualAllowed() ||
+        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(residualWellFlux) || residualWellFlux > maxResidualAllowed() ||
             std::isnan(residualWell)     || residualWell     > maxResidualAllowed() )
         {
@@ -1985,12 +1966,12 @@ namespace {
         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) << CNVW
-                  << std::setw(11) << CNVO
-                  << std::setw(11) << CNVG
+                  << 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) << residualWellFlux
                   << std::setw(11) << residualWell
                   << std::endl;