diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp
index f32aa46fb..2d505ffb8 100644
--- a/opm/autodiff/BlackoilModelBase.hpp
+++ b/opm/autodiff/BlackoilModelBase.hpp
@@ -540,6 +540,7 @@ namespace Opm {
                              int nc) const;
 
         double dpMaxRel() const { return param_.dp_max_rel_; }
+        double dbhpMaxRel() const {return param_.dbhp_max_rel_; }
         double dsMax() const { return param_.ds_max_; }
         double drMaxRel() const { return param_.dr_max_rel_; }
         double maxResidualAllowed() const { return param_.max_residual_allowed_; }
diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp
index 42f00cf07..f1a94c761 100644
--- a/opm/autodiff/BlackoilModelBase_impl.hpp
+++ b/opm/autodiff/BlackoilModelBase_impl.hpp
@@ -1064,7 +1064,7 @@ typedef Eigen::Array<double,
                 ADB::V total_residual_v = total_residual.value();
                 const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix());
                 assert(dx.size() == total_residual_v.size());
-                asImpl().wellModel().updateWellState(dx.array(), dpMaxRel(), well_state);
+                asImpl().wellModel().updateWellState(dx.array(), dbhpMaxRel(), well_state);
             }
             // We have to update the well controls regardless whether there are local
             // wells active or not as parallel logging will take place that needs to
@@ -1141,6 +1141,7 @@ typedef Eigen::Array<double,
         const V null;
         assert(null.size() == 0);
         const V zero = V::Zero(nc);
+        const V ones = V::Constant(nc,1.0);
 
         // Extract parts of dx corresponding to each part.
         const V dp = subset(dx, Span(nc));
@@ -1165,7 +1166,6 @@ typedef Eigen::Array<double,
         const V p = (p_old - dp_limited).max(zero);
         std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
 
-
         // Saturation updates.
         const Opm::PhaseUsage& pu = fluid_.phaseUsage();
         const DataBlock s_old = Eigen::Map<const DataBlock>(& reservoir_state.saturation()[0], nc, np);
@@ -1213,7 +1213,6 @@ typedef Eigen::Array<double,
             so = so_old - step * dso;
         }
 
-        // Appleyard chop process.
         if (active_[Gas]) {
             auto ixg = sg < 0;
             for (int c = 0; c < nc; ++c) {
@@ -1255,45 +1254,23 @@ typedef Eigen::Array<double,
             }
         }
 
-        //const V sumSat = sw + so + sg;
-        //sw = sw / sumSat;
-        //so = so / sumSat;
-        //sg = sg / sumSat;
-
-        // Update the reservoir_state
-        if (active_[Water]) {
-            for (int c = 0; c < nc; ++c) {
-                reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
-            }
-        }
-
-        if (active_[Gas]) {
-            for (int c = 0; c < nc; ++c) {
-                reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
-            }
-        }
-
-        if (active_[ Oil ]) {
-            for (int c = 0; c < nc; ++c) {
-                reservoir_state.saturation()[c*np + pu.phase_pos[ Oil ]] = so[c];
-            }
-        }
-
         // Update rs and rv
         const double drmaxrel = drMaxRel();
         V rs;
         if (has_disgas_) {
             const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
             const V drs = isRs_ * dxvar;
-            const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drmaxrel);
+            const V drs_limited = sign(drs) * drs.abs().min( (rs_old.abs()*drmaxrel).max( ones*1e-6));
             rs = rs_old - drs_limited;
+            rs = rs.max(zero);
         }
         V rv;
         if (has_vapoil_) {
             const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
             const V drv = isRv_ * dxvar;
-            const V drv_limited = sign(drv) * drv.abs().min(rv_old.abs()*drmaxrel);
+            const V drv_limited = sign(drv) * drv.abs().min( (rv_old.abs()*drmaxrel).max( ones*1e-6));
             rv = rv_old - drv_limited;
+            rv = rv.max(zero);
         }
 
         // Sg is used as primal variable for water only cells.
@@ -1318,11 +1295,16 @@ typedef Eigen::Array<double,
             for (int c = 0; c < nc; ++c) {
                 if (useSg[c]) {
                     rs[c] = rsSat[c];
+                    if (watOnly[c]) {
+                        so[c] = 0;
+                        sg[c] = 0;
+                        rs[c] = 0;
+                    }
                 } else {
                     hydroCarbonState[c] = HydroCarbonState::OilOnly;
                 }
             }
-
+            rs = rs.min(rsSat);
         }
 
         // phase transitions so <-> rv
@@ -1345,14 +1327,37 @@ typedef Eigen::Array<double,
             for (int c = 0; c < nc; ++c) {
                 if (useSg[c]) {
                     rv[c] = rvSat[c];
+                    if (watOnly[c]) {
+                        so[c] = 0;
+                        sg[c] = 0;
+                        rv[c] = 0;
+                    }
                 } else {
                     hydroCarbonState[c] = HydroCarbonState::GasOnly;
                 }
             }
-
+            rv = rv.min(rvSat);
         }
 
         // Update the reservoir_state
+        if (active_[Water]) {
+            for (int c = 0; c < nc; ++c) {
+                reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
+            }
+        }
+
+        if (active_[Gas]) {
+            for (int c = 0; c < nc; ++c) {
+                reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
+            }
+        }
+
+        if (active_[ Oil ]) {
+            for (int c = 0; c < nc; ++c) {
+                reservoir_state.saturation()[c*np + pu.phase_pos[ Oil ]] = so[c];
+            }
+        }
+
         if (has_disgas_) {
             std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
         }
@@ -1361,8 +1366,7 @@ typedef Eigen::Array<double,
             std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
         }
 
-
-        asImpl().wellModel().updateWellState(dwells, dpMaxRel(), well_state);
+        asImpl().wellModel().updateWellState(dwells, dbhpMaxRel(), well_state);
 
         // Update phase conditions used for property calculations.
         updatePhaseCondFromPrimalVariable(reservoir_state);
diff --git a/opm/autodiff/BlackoilModelParameters.cpp b/opm/autodiff/BlackoilModelParameters.cpp
index 6d569399b..b3f1651df 100644
--- a/opm/autodiff/BlackoilModelParameters.cpp
+++ b/opm/autodiff/BlackoilModelParameters.cpp
@@ -19,8 +19,6 @@
 
 #include <opm/autodiff/BlackoilModelParameters.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
-#include <opm/parser/eclipse/Units/Units.hpp>
-
 
 namespace Opm
 {
@@ -44,6 +42,7 @@ namespace Opm
         dp_max_rel_  = param.getDefault("dp_max_rel", dp_max_rel_);
         ds_max_      = param.getDefault("ds_max", ds_max_);
         dr_max_rel_  = param.getDefault("dr_max_rel", dr_max_rel_);
+        dbhp_max_rel_=  param.getDefault("dbhp_max_rel", dbhp_max_rel_);
         max_residual_allowed_ = param.getDefault("max_residual_allowed", max_residual_allowed_);
         tolerance_mb_    = param.getDefault("tolerance_mb", tolerance_mb_);
         tolerance_cnv_   = param.getDefault("tolerance_cnv", tolerance_cnv_);
@@ -64,9 +63,10 @@ namespace Opm
     void BlackoilModelParameters::reset()
     {
         // default values for the solver parameters
-        dp_max_rel_      = 1.0e9;
+        dp_max_rel_      = 0.2;
         ds_max_          = 0.2;
         dr_max_rel_      = 1.0e9;
+        dbhp_max_rel_    = 1.0;
         max_residual_allowed_ = 1e7;
         tolerance_mb_    = 1.0e-5;
         tolerance_cnv_   = 1.0e-2;
diff --git a/opm/autodiff/BlackoilModelParameters.hpp b/opm/autodiff/BlackoilModelParameters.hpp
index c616b547d..037bfdfb6 100644
--- a/opm/autodiff/BlackoilModelParameters.hpp
+++ b/opm/autodiff/BlackoilModelParameters.hpp
@@ -36,6 +36,8 @@ namespace Opm
         double ds_max_;
         /// Max relative change in gas-oil or oil-gas ratio in single iteration.
         double dr_max_rel_;
+        /// Max relative change in bhp in single iteration.
+        double dbhp_max_rel_;
         /// Absolute max limit for residuals.
         double max_residual_allowed_;
         /// Relative mass balance tolerance (total mass balance error).