diff --git a/opm/autodiff/BlackoilPropsAd.cpp b/opm/autodiff/BlackoilPropsAd.cpp
index 22288afff..50025ff16 100644
--- a/opm/autodiff/BlackoilPropsAd.cpp
+++ b/opm/autodiff/BlackoilPropsAd.cpp
@@ -609,6 +609,23 @@ namespace Opm
         OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
     }
 
+    /// Bubble point curve for Rs as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    V BlackoilPropsAd::rsSat(const V& po,
+                             const V& so,
+                             const Cells& cells) const
+    {
+        // Suppress warning about "unused parameters".
+        static_cast<void>(po);
+        static_cast<void>(so);
+        static_cast<void>(cells);
+
+        OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
+    }
+
     /// Bubble point curve for Rs as function of oil pressure.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -623,6 +640,23 @@ namespace Opm
         OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
     }
 
+    /// Bubble point curve for Rs as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    ADB BlackoilPropsAd::rsSat(const ADB& po,
+                               const ADB& so,
+                               const Cells& cells) const
+    {
+        // Suppress warning about "unused parameters".
+        static_cast<void>(po);
+        static_cast<void>(so);
+        static_cast<void>(cells);
+
+        OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
+    }
+
     // ------ Rs bubble point curve ------
 
     /// Bubble point curve for Rs as function of oil pressure.
@@ -639,6 +673,23 @@ namespace Opm
         OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
     }
 
+    /// Bubble point curve for Rs as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    V BlackoilPropsAd::rvSat(const V& po,
+                             const V& so,
+                             const Cells& cells) const
+    {
+        // Suppress warning about "unused parameters".
+        static_cast<void>(po);
+        static_cast<void>(so);
+        static_cast<void>(cells);
+
+        OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
+    }
+
     /// Bubble point curve for Rs as function of oil pressure.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -653,6 +704,23 @@ namespace Opm
         OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
     }
 
+    /// Bubble point curve for Rs as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    ADB BlackoilPropsAd::rvSat(const ADB& po,
+                               const ADB& so,
+                               const Cells& cells) const
+    {
+        // Suppress warning about "unused parameters".
+        static_cast<void>(po);
+        static_cast<void>(so);
+        static_cast<void>(cells);
+
+        OPM_THROW(std::runtime_error, "Method rsMax() not implemented.");
+    }
+
     // ------ Relative permeability ------
 
     /// Relative permeabilities for all phases.
@@ -828,7 +896,12 @@ namespace Opm
     {
         OPM_THROW(std::logic_error, "BlackoilPropsAd class does not support hysteresis.");
     }
-
+    
+    /// Update for max oil saturation.
+    void BlackoilPropsAd::updateSatOilMax(const std::vector<double>& saturation)
+    {
+        OPM_THROW(std::logic_error, "BlackoilPropsAd class does not support this functionality.");
+    }
 
 } // namespace Opm
 
diff --git a/opm/autodiff/BlackoilPropsAd.hpp b/opm/autodiff/BlackoilPropsAd.hpp
index b972ab76c..e7656f6ad 100644
--- a/opm/autodiff/BlackoilPropsAd.hpp
+++ b/opm/autodiff/BlackoilPropsAd.hpp
@@ -254,6 +254,15 @@ namespace Opm
         V rsSat(const V& po,
                 const Cells& cells) const;
 
+        /// Solution gas/oil ratio and its derivatives at saturated condition as a function of p.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        V rsSat(const V& po,
+                const V& so,
+                const Cells& cells) const;
+
         /// Solution gas/oil ratio and its derivatives at saturated condition as a function of p.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -261,6 +270,15 @@ namespace Opm
         ADB rsSat(const ADB& po,
                   const Cells& cells) const;
 
+        /// Solution gas/oil ratio and its derivatives at saturated condition as a function of p.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        ADB rsSat(const ADB& po,
+                  const ADB& so,
+                  const Cells& cells) const;
+
         // ------ Rv condensation curve ------
 
         /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
@@ -270,6 +288,15 @@ namespace Opm
         V rvSat(const V& po,
                 const Cells& cells) const;
 
+        /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        V rvSat(const V& po,
+                const V& so,
+                const Cells& cells) const;
+
         /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -277,6 +304,15 @@ namespace Opm
         ADB rvSat(const ADB& po,
                   const Cells& cells) const;
 
+        /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        ADB rvSat(const ADB& po,
+                  const ADB& so,
+                  const Cells& cells) const;
+
         // ------ Relative permeability ------
 
         /// Relative permeabilities for all phases.
@@ -320,6 +356,11 @@ namespace Opm
         /// \param[in]  cells       Array of n cell indices to be associated with the saturation values.
         void updateSatHyst(const std::vector<double>& saturation,
                            const std::vector<int>& cells);
+
+
+        /// Update for max oil saturation.
+        void updateSatOilMax(const std::vector<double>& saturation);
+
     private:
         const BlackoilPropertiesInterface& props_;
         PhaseUsage pu_;
diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.cpp b/opm/autodiff/BlackoilPropsAdFromDeck.cpp
index ae5c9e186..d811346e5 100644
--- a/opm/autodiff/BlackoilPropsAdFromDeck.cpp
+++ b/opm/autodiff/BlackoilPropsAdFromDeck.cpp
@@ -188,6 +188,15 @@ namespace Opm
                 OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
             }
         }
+        // Oil vaporization controls (kw VAPPARS)
+        vap1_ = vap2_ = 0.0;
+        if (deck->hasKeyword("VAPPARS") && deck->hasKeyword("VAPOIL") && deck->hasKeyword("DISGAS")) {
+            vap1_ = deck->getKeyword("VAPPARS")->getRecord(0)->getItem(0)->getRawDouble(0);
+            vap2_ = deck->getKeyword("VAPPARS")->getRecord(0)->getItem(1)->getRawDouble(0);
+            satOilMax_.resize(number_of_cells, 0.0);
+        } else if (deck->hasKeyword("VAPPARS")) {
+            OPM_THROW(std::runtime_error, "Input has VAPPARS, but missing VAPOIL and/or DISGAS\n");
+        }
 
         SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
             = new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
@@ -199,6 +208,7 @@ namespace Opm
                       "Inconsistent number of phases in pvt data (" << phase_usage_.num_phases
                       << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
         }
+        vap_satmax_guard_ = 0.01;
     }
 
     ////////////////////////////
@@ -746,6 +756,20 @@ namespace Opm
         return rbub;
     }
 
+    /// Bubble point curve for Rs as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    V BlackoilPropsAdFromDeck::rsSat(const V& po,
+                                     const V& so,
+                                     const Cells& cells) const
+    {
+        V rs = rsSat(po, cells);
+        applyVap(rs, so, cells, vap2_);
+        return rs;
+    }
+
     /// Bubble point curve for Rs as function of oil pressure.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -770,6 +794,20 @@ namespace Opm
         return ADB::function(rbub, jacs);
     }
 
+    /// Bubble point curve for Rs as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    ADB BlackoilPropsAdFromDeck::rsSat(const ADB& po,
+                                       const ADB& so,
+                                       const Cells& cells) const
+    {
+        ADB rs = rsSat(po, cells);
+        applyVap(rs, so, cells, vap2_);
+        return rs;
+    }
+
     // ------ Condensation curve ------
 
     /// Condensation curve for Rv as function of oil pressure.
@@ -790,6 +828,20 @@ namespace Opm
         return rv;
     }
 
+    /// Condensation curve for Rv as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    V BlackoilPropsAdFromDeck::rvSat(const V& po,
+                                     const V& so,
+                                     const Cells& cells) const
+    {
+        V rv = rvSat(po, cells);
+        applyVap(rv, so, cells, vap1_);
+        return rv;
+    }
+
     /// Condensation curve for Rv as function of oil pressure.
     /// \param[in]  po     Array of n oil pressure values.
     /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -814,6 +866,20 @@ namespace Opm
         return ADB::function(rv, jacs);
     }
 
+    /// Condensation curve for Rv as function of oil pressure.
+    /// \param[in]  po     Array of n oil pressure values.
+    /// \param[in]  so     Array of n oil saturation values.
+    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+    /// \return            Array of n bubble point values for Rs.
+    ADB BlackoilPropsAdFromDeck::rvSat(const ADB& po,
+                                       const ADB& so,
+                                       const Cells& cells) const
+    {
+        ADB rv = rvSat(po, cells);
+        applyVap(rv, so, cells, vap1_);
+        return rv;
+    }
+
     // ------ Relative permeability ------
 
     /// Relative permeabilities for all phases.
@@ -987,6 +1053,76 @@ namespace Opm
         const int n = cells.size();
         satprops_->updateSatHyst(n, cells.data(), saturation.data());
     }
+    
+    /// Update for max oil saturation.
+    void BlackoilPropsAdFromDeck::updateSatOilMax(const std::vector<double>& saturation)
+    {
+        if (!satOilMax_.empty()) {
+            const int n = satOilMax_.size();
+            const int np = phase_usage_.num_phases;
+            const int posOil = phase_usage_.phase_pos[Oil];
+            const double* s = saturation.data();
+            for (int i=0; i<n; ++i) {
+                if (satOilMax_[i] < s[np*i+posOil]) {
+                    satOilMax_[i] = s[np*i+posOil];
+                }
+            }
+        }
+    }
+    
+    /// Apply correction to rs/rv according to kw VAPPARS
+    /// \param[in/out] r     Array of n rs/rv values.
+    /// \param[in]     so    Array of n oil saturation values.
+    /// \param[in]     cells Array of n cell indices to be associated with the r and so values.
+    /// \param[in]     vap   Correction parameter.
+    void BlackoilPropsAdFromDeck::applyVap(V& r,
+                                           const V& so,
+                                           const std::vector<int>& cells,
+                                           const double vap) const
+    {
+        if (!satOilMax_.empty() && vap > 0.0) { 
+            const int n = cells.size();
+            V factor = V::Ones(n, 1);
+            for (int i=0; i<n; ++i) {
+                if (satOilMax_[cells[i]] > vap_satmax_guard_ && so[i] < satOilMax_[cells[i]]) {
+                    factor[i] = std::pow(so[i]/satOilMax_[cells[i]], vap);
+                } 
+            }
+            r = factor*r;
+        }
+    }
+    
+    /// Apply correction to rs/rv according to kw VAPPARS
+    /// \param[in/out] r     Array of n rs/rv values.
+    /// \param[in]     so    Array of n oil saturation values.
+    /// \param[in]     cells Array of n cell indices to be associated with the r and so values.
+    /// \param[in]     vap   Correction parameter.
+    void BlackoilPropsAdFromDeck::applyVap(ADB& r,
+                                           const ADB& so,
+                                           const std::vector<int>& cells,
+                                           const double vap) const
+    {
+        if (!satOilMax_.empty() && vap > 0.0) { 
+            const int n = cells.size();
+            V factor = V::Ones(n, 1);
+            //V dfactor_dso = V::Zero(n, 1);  TODO: Consider effect of complete jacobian (including so-derivatives)
+            for (int i=0; i<n; ++i) {
+                if (satOilMax_[cells[i]] > vap_satmax_guard_ && so.value()[i] < satOilMax_[cells[i]]) {
+                    factor[i] = std::pow(so.value()[i]/satOilMax_[cells[i]], vap);
+                    //dfactor_dso[i] = vap*std::pow(so.value()[i]/satOilMax_[cells[i]], vap-1.0)/satOilMax_[cells[i]];
+                }
+            }
+            //ADB::M dfactor_dso_diag = spdiag(dfactor_dso);
+            //const int num_blocks = so.numBlocks();
+            //std::vector<ADB::M> jacs(num_blocks);
+            //for (int block = 0; block < num_blocks; ++block) {
+            //    jacs[block] = dfactor_dso_diag * so.derivative()[block];
+            //}
+            //r = ADB::function(factor, jacs)*r;
+            
+            r = factor*r;
+        }
+    }
 
 } // namespace Opm
 
diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.hpp b/opm/autodiff/BlackoilPropsAdFromDeck.hpp
index bced83fa9..8042c4d4d 100644
--- a/opm/autodiff/BlackoilPropsAdFromDeck.hpp
+++ b/opm/autodiff/BlackoilPropsAdFromDeck.hpp
@@ -277,6 +277,15 @@ namespace Opm
         V rsSat(const V& po,
                 const Cells& cells) const;
 
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        V rsSat(const V& po,
+                const V& so,
+                const Cells& cells) const;
+
         /// Bubble point curve for Rs as function of oil pressure.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -284,6 +293,15 @@ namespace Opm
         ADB rsSat(const ADB& po,
                   const Cells& cells) const;
 
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        ADB rsSat(const ADB& po,
+                  const ADB& so,
+                  const Cells& cells) const;
+
         // ------ Rv condensation curve ------
 
         /// Condensation curve for Rv as function of oil pressure.
@@ -293,6 +311,15 @@ namespace Opm
         V rvSat(const V& po,
                 const Cells& cells) const;
 
+        /// Condensation curve for Rv as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        V rvSat(const V& po,
+                const V& so,
+                const Cells& cells) const;
+
         /// Condensation curve for Rv as function of oil pressure.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -300,6 +327,15 @@ namespace Opm
         ADB rvSat(const ADB& po,
                   const Cells& cells) const;
 
+        /// Condensation curve for Rv as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        ADB rvSat(const ADB& po,
+                  const ADB& so,
+                  const Cells& cells) const;
+
         // ------ Relative permeability ------
 
         /// Relative permeabilities for all phases.
@@ -344,6 +380,9 @@ namespace Opm
         void updateSatHyst(const std::vector<double>& saturation,
                            const std::vector<int>& cells);
 
+        /// Update for max oil saturation.                  
+        void updateSatOilMax(const std::vector<double>& saturation);
+
     private:
         /// Initializes the properties.
         template <class CentroidIterator>
@@ -356,6 +395,17 @@ namespace Opm
                   int dimension,
                   const bool init_rock);
 
+        /// Correction to rs/rv according to kw VAPPARS
+        void applyVap(V& r,
+                      const V& so,
+                      const std::vector<int>& cells,
+                      const double vap) const;
+
+        void applyVap(ADB& r,
+                      const ADB& so,
+                      const std::vector<int>& cells,
+                      const double vap) const;
+
         RockFromDeck rock_;
         std::unique_ptr<SaturationPropsInterface> satprops_;
 
@@ -380,6 +430,13 @@ namespace Opm
         std::vector<int> pvtTableIdx_;
 
         std::vector<std::array<double, BlackoilPhases::MaxNumPhases> > densities_;
+        
+        // VAPPARS
+        double vap1_;
+        double vap2_;
+        std::vector<double> satOilMax_;
+        double vap_satmax_guard_;  //Threshold value to promote stability
+
     };
 } // namespace Opm
 
diff --git a/opm/autodiff/BlackoilPropsAdInterface.hpp b/opm/autodiff/BlackoilPropsAdInterface.hpp
index bde9f2b7f..f15b33ca1 100644
--- a/opm/autodiff/BlackoilPropsAdInterface.hpp
+++ b/opm/autodiff/BlackoilPropsAdInterface.hpp
@@ -247,6 +247,16 @@ namespace Opm
         V rsSat(const V& po,
                 const Cells& cells) const = 0;
 
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        virtual
+        V rsSat(const V& po,
+                const V& so,
+                const Cells& cells) const = 0;
+
         /// Bubble point curve for Rs as function of oil pressure.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -255,6 +265,16 @@ namespace Opm
         ADB rsSat(const ADB& po,
                   const Cells& cells) const = 0;
 
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        virtual
+        ADB rsSat(const ADB& po,
+                  const ADB& so,
+                  const Cells& cells) const = 0;
+
         // ------ Rs bubble point curve ------
 
         /// Bubble point curve for Rs as function of oil pressure.
@@ -265,6 +285,16 @@ namespace Opm
         V rvSat(const V& po,
                 const Cells& cells) const = 0;
 
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        virtual
+        V rvSat(const V& po,
+                const V& so,
+                const Cells& cells) const = 0;
+
         /// Bubble point curve for Rs as function of oil pressure.
         /// \param[in]  po     Array of n oil pressure values.
         /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
@@ -273,6 +303,16 @@ namespace Opm
         ADB rvSat(const ADB& po,
                   const Cells& cells) const = 0;
 
+        /// Bubble point curve for Rs as function of oil pressure.
+        /// \param[in]  po     Array of n oil pressure values.
+        /// \param[in]  so     Array of n oil saturation values.
+        /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
+        /// \return            Array of n bubble point values for Rs.
+        virtual
+        ADB rvSat(const ADB& po,
+                  const ADB& so,
+                  const Cells& cells) const = 0;
+
         // ------ Relative permeability ------
 
         /// Relative permeabilities for all phases.
@@ -321,6 +361,10 @@ namespace Opm
         virtual
         void updateSatHyst(const std::vector<double>& saturation,
                            const std::vector<int>& cells) = 0;
+
+        /// Update for max oil saturation.
+        virtual                  
+        void updateSatOilMax(const std::vector<double>& saturation) = 0;
     };
 
 } // namespace Opm
diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp
index 6cb85380f..f46413b42 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp
@@ -254,18 +254,22 @@ namespace Opm {
 
         V
         fluidRsSat(const V&                p,
+                   const V&                so,
                    const std::vector<int>& cells) const;
 
         ADB
         fluidRsSat(const ADB&              p,
+                   const ADB&              so,
                    const std::vector<int>& cells) const;
 
         V
         fluidRvSat(const V&                p,
+                   const V&                so,
                    const std::vector<int>& cells) const;
 
         ADB
         fluidRvSat(const ADB&              p,
+                   const ADB&              so,
                    const std::vector<int>& cells) const;
 
         ADB
diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
index 621c8c001..5a204de3b 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
@@ -476,12 +476,12 @@ namespace {
 
             if (active_[ Gas]) {
                 // Define Sg Rs and Rv in terms of xvar.
-                const ADB& rsSat = fluidRsSat(state.pressure, cells_);
-                const ADB& rvSat = fluidRvSat(state.pressure, cells_);
                 const ADB& xvar = vars[ nextvar++ ];
                 const ADB& sg = isSg*xvar + isRv* so;
                 state.saturation[ pu.phase_pos[ Gas ] ] = sg;
                 so = so - sg;
+                const ADB rsSat = fluidRsSat(state.pressure, so, cells_);
+                const ADB rvSat = fluidRvSat(state.pressure, so, cells_);
 
                 if (has_disgas_) {
                     state.rs = (1-isRs) * rsSat + isRs*xvar;
@@ -583,18 +583,20 @@ namespace {
             const ADB bw = fluid_.bWat(perf_press, well_cells);
             b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw.value();
         }
+        assert(active_[Oil]);
+        const ADB perf_so =  subset(state.saturation[pu.phase_pos[Oil]], well_cells);
         if (pu.phase_used[BlackoilPhases::Liquid]) {
             const ADB perf_rs = subset(state.rs, well_cells);
             const ADB bo = fluid_.bOil(perf_press, perf_rs, perf_cond, well_cells);
             b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo.value();
-            const V rssat = fluidRsSat(perf_press.value(), well_cells);
+            const V rssat = fluidRsSat(perf_press.value(), perf_so.value(), well_cells);
             rssat_perf.assign(rssat.data(), rssat.data() + nperf);
         }
         if (pu.phase_used[BlackoilPhases::Vapour]) {
             const ADB perf_rv = subset(state.rv, well_cells);
             const ADB bg = fluid_.bGas(perf_press, perf_rv, perf_cond, well_cells);
             b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg.value();
-            const V rvsat = fluidRvSat(perf_press.value(), well_cells);
+            const V rvsat = fluidRvSat(perf_press.value(), perf_so.value(), well_cells);
             rvsat_perf.assign(rvsat.data(), rvsat.data() + nperf);
         }
         // b is row major, so can just copy data.
@@ -1273,8 +1275,8 @@ namespace {
 
 
         // phase translation sg <-> rs
-        const V rsSat0 = fluidRsSat(p_old, cells_);
-        const V rsSat = fluidRsSat(p, cells_);
+        const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
+        const V rsSat = fluidRsSat(p, so, cells_);
 
         std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg);
 
@@ -1298,8 +1300,8 @@ namespace {
         }
 
         // phase transitions so <-> rv
-        const V rvSat0 = fluidRvSat(p_old, cells_);
-        const V rvSat = fluidRvSat(p, cells_);
+        const V rvSat0 = fluidRvSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
+        const V rvSat = fluidRvSat(p, so, cells_);
 
         if (has_vapoil_) {
             // The obvious case
@@ -1851,9 +1853,10 @@ namespace {
     template<class T>
     V
     FullyImplicitBlackoilSolver<T>::fluidRsSat(const V&                p,
+                                            const V&                satOil,
                                             const std::vector<int>& cells) const
     {
-        return fluid_.rsSat(p, cells);
+        return fluid_.rsSat(p, satOil, cells);
     }
 
 
@@ -1863,17 +1866,19 @@ namespace {
     template<class T>
     ADB
     FullyImplicitBlackoilSolver<T>::fluidRsSat(const ADB&              p,
+                                            const ADB&              satOil,
                                             const std::vector<int>& cells) const
     {
-        return fluid_.rsSat(p, cells);
+        return fluid_.rsSat(p, satOil, cells);
     }
 
     template<class T>
     V
     FullyImplicitBlackoilSolver<T>::fluidRvSat(const V&                p,
+                                            const V&              satOil,
                                             const std::vector<int>& cells) const
     {
-        return fluid_.rvSat(p, cells);
+        return fluid_.rvSat(p, satOil, cells);
     }
 
 
@@ -1883,9 +1888,10 @@ namespace {
     template<class T>
     ADB
     FullyImplicitBlackoilSolver<T>::fluidRvSat(const ADB&              p,
+                                            const ADB&              satOil,
                                             const std::vector<int>& cells) const
     {
-        return fluid_.rvSat(p, cells);
+        return fluid_.rvSat(p, satOil, cells);
     }
 
 
diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
index cbac54486..ab11cbc0d 100644
--- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
+++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
@@ -302,10 +302,14 @@ namespace Opm
             SimulatorReport sreport;
 
             FullyImplicitBlackoilSolver<T> solver(param_, grid_, props_, geo_, rock_comp_props_, *wells, solver_, has_disgas_, has_vapoil_);
+            // Max oil saturation
+            props_.updateSatOilMax(state.saturation());
+
+            // Hysteresis
+            props_.updateSatHyst(state.saturation(), allcells_);
 
             // Run solver.
             solver_timer.start();
-            std::vector<double> initial_pressure = state.pressure();
             solver.step(timer.currentStepLength(), state, well_state);
 
             // Stop timer and report.
@@ -316,9 +320,6 @@ namespace Opm
             stime += st;
             sreport.pressure_time = st;
 
-            // Hysteresis
-            props_.updateSatHyst(state.saturation(), allcells_);
-
             sreport.total_time =  step_timer.secsSinceStart();
             if (output_) {
                 sreport.reportParam(tstep_os);