diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
index e106cc7bf..747e5cc58 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
@@ -2052,29 +2052,35 @@ namespace {
         const PhaseUsage& pu = fluid_.phaseUsage();
         const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
 
+        // Water/Oil/Gas system
         assert (active_[ Gas ]);
 
         // reset the primary variables if RV and RS is not set Sg is used as primary variable.
         primalVariable_.resize(nc);
         std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg);
 
-        if (has_disgas_) {
-            // Oil/Gas or Water/Oil/Gas system
-            const V sg = s.col(pu.phase_pos[ Gas ]);
-            const V so = s.col(pu.phase_pos[ Oil ]);
+        const V sg = s.col(pu.phase_pos[ Gas ]);
+        const V so = s.col(pu.phase_pos[ Oil ]);
+        const V sw = s.col(pu.phase_pos[ Water ]);
 
+        const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
+        auto watOnly = sw >  (1 - epsilon);
+        auto hasOil = so > 0;
+        auto hasGas = sg > 0;
+
+        // For oil only cells Rs is used as primal variable. For cells almost full of water
+        // the default primal variable (Sg) is used.
+        if (has_disgas_) {         
             for (V::Index c = 0, e = sg.size(); c != e; ++c) {
-                if ( sg[c] <= 0 && so[c] > 0 ) {primalVariable_[c] = PrimalVariables::RS; }
+                if ( !watOnly[c] && hasOil[c] && !hasGas[c] ) {primalVariable_[c] = PrimalVariables::RS; }
             }
         }
 
+        // For gas only cells Rv is used as primal variable. For cells almost full of water
+        // the default primal variable (Sg) is used.
         if (has_vapoil_) {
-            // Oil/Gas or Water/Oil/Gas system
-            const V sg = s.col(pu.phase_pos[ Gas ]);
-            const V so = s.col(pu.phase_pos[ Oil ]);
-
             for (V::Index c = 0, e = so.size(); c != e; ++c) {
-                if (so[c] <= 0 && sg[c] > 0) {primalVariable_[c] = PrimalVariables::RV; }
+                if ( !watOnly[c] && hasGas[c] && !hasOil[c] ) {primalVariable_[c] = PrimalVariables::RV; }
             }
         }
         updatePhaseCondFromPrimalVariable();