diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 54eb9daa..963a9a2a 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -172,6 +172,7 @@ list (APPEND TEST_SOURCE_FILES
 	tests/test_parallel_linearsolver.cpp
 	tests/test_param.cpp
 	tests/test_blackoilfluid.cpp
+	tests/test_satfunc.cpp
 	tests/test_shadow.cpp
 	tests/test_equil.cpp
 	tests/test_regionmapping.cpp
@@ -200,6 +201,12 @@ list (APPEND TEST_DATA_FILES
 	tests/equil_liveoil.DATA
 	tests/equil_rsvd_and_rvvd.DATA
 	tests/wetgas.DATA
+	tests/satfuncStandard.DATA
+	tests/satfuncEPSBase.DATA
+	tests/satfuncEPS_A.DATA
+	tests/satfuncEPS_B.DATA
+	tests/satfuncEPS_C.DATA
+	tests/satfuncEPS_D.DATA
 	tests/testBlackoilState1.DATA
 	tests/testBlackoilState2.DATA
 	tests/wells_manager_data.data
diff --git a/opm/core/props/satfunc/SatFuncGwseg.hpp b/opm/core/props/satfunc/SatFuncGwseg.hpp
index aaec8fa9..4a6d0fed 100644
--- a/opm/core/props/satfunc/SatFuncGwseg.hpp
+++ b/opm/core/props/satfunc/SatFuncGwseg.hpp
@@ -360,8 +360,8 @@ namespace Opm
 
             const double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
             const double krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
-            const double krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_);
-            const double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_);
+            const double krow = epst->watoil.scaleKr(ssow, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_);
+            const double krog = epst->gasoil.scaleKr(ssog, this->krog_(std::max(1.0-_ssog-_swco,eps)), this->krorg_);
 
             // xw and xg are the fractions occupied by water and gas zones.
             const double xw = (sw - swco) / d;
@@ -373,8 +373,8 @@ namespace Opm
             // Derivatives.
             double dkrww = _dsdsw*epst->wat.scaleKrDeriv(sw, this->krw_.derivative(_sw));
             double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(sg, this->krg_.derivative(_sg));
-            double dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(1.0-_ssow));
-            double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(1.0-_ssog-_swco));
+            double dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(std::max(1.0-_ssow,_swco)));
+            double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(std::max(1.0-_ssog-_swco,eps)));           
             dkrds[BlackoilPhases::Aqua   + BlackoilPhases::Aqua*np]   =  dkrww;
             dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np]   =  (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
             dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
@@ -469,21 +469,21 @@ namespace Opm
             if (ssow >= sat_hyst->sow_hyst) { // Drainage
                 double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
                 double _dsdssow = epst->watoil.scaleSatDeriv(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
-                krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_);
-                dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(1.0-_ssow));
+                krow = epst->watoil.scaleKr(ssow, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_);
+                dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(std::max(1.0-_ssow,_swco)));
             } else { // Imbibition
                 double ssow_shifted = ssow + sat_hyst->sow_shift;
                 double _ssow = epst_hyst->watoil.scaleSat(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
                 double _dsdssow = epst_hyst->watoil.scaleSatDeriv(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
-                krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(1.0-_ssow), this->krorw_);
-                dkrow = _dsdssow*epst_hyst->watoil.scaleKrDeriv(ssow_shifted, this->krow_.derivative(1.0-_ssow));
+                krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_);
+                dkrow = _dsdssow*epst_hyst->watoil.scaleKrDeriv(ssow_shifted, this->krow_.derivative(std::max(1.0-_ssow,_swco)));
             }
             
             // Oil in gas and connate water - use drainage curve only 
             double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
             double _dsdssog = epst->gasoil.scaleSatDeriv(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
-            double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_);
-            double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(1.0-_ssog-_swco));
+            double krog = epst->gasoil.scaleKr(ssog, this->krog_(std::max(1.0-_ssog-_swco,eps)), this->krorg_);
+            double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(std::max(1.0-_ssog-_swco,eps)));
 
             // xw and xg are the fractions occupied by water and gas zones.
             const double xw = (sw - swco) / d;
diff --git a/opm/core/props/satfunc/SatFuncSimple.cpp b/opm/core/props/satfunc/SatFuncSimple.cpp
index 91200053..dd5de753 100644
--- a/opm/core/props/satfunc/SatFuncSimple.cpp
+++ b/opm/core/props/satfunc/SatFuncSimple.cpp
@@ -159,7 +159,7 @@ namespace Opm
     double EPSTransforms::Transform::scaleKrDeriv(double s, double krDeriv)  const
     {
         if (doKrCrit) {
-            if (s <= scr) {
+            if (s < scr) {
                 return 0.0;
             } else if (s <= sr) {
                 return krDeriv*krSlopeCrit;
@@ -172,7 +172,7 @@ namespace Opm
                 return 0.0;
             }
         } else if (doKrMax) {
-            if (s <= scr) {
+            if (s < scr) {
                 return 0.0;
             } else if (s <= smax) {
                 return krDeriv*krSlopeMax;
@@ -180,7 +180,7 @@ namespace Opm
                 return 0.0;
             }
         } else {
-            if (s <= scr) {
+            if (s < scr) {
                 return 0.0;
             } else if (s <= smax) {
                 return krDeriv;
diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp
index 5b4a2fcb..a4667f8f 100644
--- a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp
+++ b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp
@@ -60,7 +60,7 @@ namespace Opm
         /// \param[in]  samples  Number of uniform sample points for saturation tables.
         /// NOTE: samples will only be used with the SatFuncSetUniform template argument.
         void init(Opm::DeckConstPtr deck,
-                  Opm::EclipseStateConstPtr eclState,
+                  Opm::EclipseStateConstPtr eclipseState,
                   const UnstructuredGrid& grid,
                   const int samples);
 
@@ -79,7 +79,7 @@ namespace Opm
         /// NOTE: samples will only be used with the SatFuncSetUniform template argument.
         template<class T>
         void init(Opm::DeckConstPtr deck,
-                  Opm::EclipseStateConstPtr eclState,
+                  Opm::EclipseStateConstPtr eclipseState,
                   int number_of_cells,
                   const int* global_cell,
                   const T& begin_cell_centroids,
@@ -165,14 +165,9 @@ namespace Opm
                      int number_of_cells,
                      const int* global_cell,
                      const T& begin_cell_centroids,
-                     int dimensions);
-        template<class T>
-        void initEPSHyst(Opm::DeckConstPtr deck,
-                         Opm::EclipseStateConstPtr eclipseState,
-                         int number_of_cells,
-                         const int* global_cell,
-                         const T& begin_cell_centroids,
-                         int dimensions);
+                     int dimensions,
+                     const std::vector<std::string>& eps_kw,
+                     std::vector<EPSTransforms>& eps_transf);
         template<class T>
         void initEPSKey(Opm::DeckConstPtr deck,
                         Opm::EclipseStateConstPtr eclipseState,
@@ -204,8 +199,8 @@ namespace Opm
 
         bool columnIsMasked_(Opm::DeckConstPtr deck,
                              const std::string& keywordName,
-                             int /* columnIdx */)
-        { return deck->getKeyword(keywordName)->getRecord(0)->getItem(0)->getSIDouble(0) != -1.0; }
+                             int columnIdx)
+        { return deck->getKeyword(keywordName)->getRecord(columnIdx)->getItem(0)->getSIDouble(0) != -1.0; }
     };
 
 
diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
index c65e1185..ed03ceb2 100644
--- a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
+++ b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
@@ -31,6 +31,7 @@
 #include <opm/parser/eclipse/Utility/ScalecrsWrapper.hpp>
 
 #include <iostream>
+#include <map>
 
 namespace Opm
 {
@@ -48,11 +49,11 @@ namespace Opm
     /// Initialize from deck.
     template <class SatFuncSet>
     void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr deck,
-                                                   Opm::EclipseStateConstPtr eclState,
+                                                   Opm::EclipseStateConstPtr eclipseState,
                                                    const UnstructuredGrid& grid,
                                                    const int samples)
     {
-        this->init(deck, eclState, grid.number_of_cells,
+        this->init(deck, eclipseState, grid.number_of_cells,
                    grid.global_cell, grid.cell_centroids,
                    grid.dimensions, samples);
     }
@@ -61,7 +62,7 @@ namespace Opm
     template <class SatFuncSet>
     template<class T>
     void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr deck,
-                                                   Opm::EclipseStateConstPtr eclState,
+                                                   Opm::EclipseStateConstPtr eclipseState,
                                                    int number_of_cells,
                                                    const int* global_cell,
                                                    const T& begin_cell_centroids,
@@ -128,7 +129,7 @@ namespace Opm
         // Initialize tables.
         satfuncset_.resize(num_tables);
         for (int table = 0; table < num_tables; ++table) {
-            satfuncset_[table].init(eclState, table, phase_usage_, samples);
+            satfuncset_[table].init(eclipseState, table, phase_usage_, samples);
         }
         
         // Check EHYSTR status
@@ -190,8 +191,11 @@ namespace Opm
             // TODO: ENPTVD/ENKRVD: Too few tables gives a cryptical message from parser, 
             //       superfluous tables are ignored by the parser without any warning ...
 
-            initEPS(deck, eclState, number_of_cells, global_cell, begin_cell_centroids,
-                    dimensions);
+            const std::vector<std::string> eps_kw{"SWL", "SWU", "SWCR", "SGL", "SGU", "SGCR", "SOWCR",
+                "SOGCR", "KRW", "KRG", "KRO", "KRWR", "KRGR", "KRORW", "KRORG", "PCW", "PCG"};
+            eps_transf_.resize(number_of_cells);
+            initEPS(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroids,
+                    dimensions, eps_kw, eps_transf_);
 
             if (do_hyst_) {
                 if (deck->hasKeyword("KRW")
@@ -230,15 +234,20 @@ namespace Opm
                     // TODO: Make actual use of IMBNUM.  For now we just consider the imbibition curve
                     //       to be a scaled version of the drainage curve (confer Norne model).
                 }
-                
-                initEPSHyst(deck, eclState, number_of_cells, global_cell, begin_cell_centroids,
-                            dimensions);
+
+                const std::vector<std::string> eps_i_kw{"ISWL", "ISWU", "ISWCR", "ISGL", "ISGU", "ISGCR", "ISOWCR",
+                    "ISOGCR", "IKRW", "IKRG", "IKRO", "IKRWR", "IKRGR", "IKRORW", "IKRORG", "IPCW", "IPCG"};
+                eps_transf_hyst_.resize(number_of_cells);
+                sat_hyst_.resize(number_of_cells);                
+                initEPS(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroids,
+                        dimensions, eps_i_kw, eps_transf_hyst_);
             }
         }
     }
 
 
 
+
     /// \return   P, the number of phases.
     template <class SatFuncSet>
     int SaturationPropsFromDeck<SatFuncSet>::numPhases() const
@@ -462,204 +471,17 @@ namespace Opm
                                                       int number_of_cells,
                                                       const int* global_cell,
                                                       const T& begin_cell_centroid,
-                                                      int dimensions)
+                                                      int dimensions,
+                                                      const std::vector<std::string>& eps_kw,
+                                                      std::vector<EPSTransforms>& eps_transf)
     {
-        std::vector<double> swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr;
-        std::vector<double> krw, krg, kro, krwr, krgr, krorw, krorg;
-        std::vector<double> pcw, pcg;
+        std::vector<std::vector<double> > eps_vec(eps_kw.size());
         const std::vector<double> dummy;
-        // Initialize saturation scaling parameter
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("SWL"),   swl);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("SWU"),   swu);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("SWCR"),  swcr);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("SGL"),   sgl);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("SGU"),   sgu);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("SGCR"),  sgcr);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("SOWCR"), sowcr);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("SOGCR"), sogcr);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("KRW"),   krw);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("KRG"),   krg);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("KRO"),   kro);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("KRWR"),  krwr);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("KRGR"),  krgr);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("KRORW"), krorw);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("KRORG"), krorg);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("PCW"), pcw);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("PCG"), pcg);
-
-        eps_transf_.resize(number_of_cells);
-
-        const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
-        const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
-        const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour];
-        const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
-        const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
-
-        for (int cell = 0; cell < number_of_cells; ++cell) {
-            if (oilWater) {
-                // ### krw
-                initEPSParam(cell, eps_transf_[cell].wat, false, 
-                             funcForCell(cell).smin_[wpos],
-                             funcForCell(cell).swcr_, 
-                             funcForCell(cell).smax_[wpos],
-                             funcForCell(cell).sowcr_, 
-                             -1.0, 
-                             funcForCell(cell).krwr_, 
-                             funcForCell(cell).krwmax_, 
-                             funcForCell(cell).pcwmax_, 
-                             swl, swcr, swu, sowcr, sgl, krwr, krw, pcw);
-                // ### krow
-                initEPSParam(cell, eps_transf_[cell].watoil, true, 
-                             0.0, 
-                             funcForCell(cell).sowcr_, 
-                             funcForCell(cell).smin_[wpos],
-                             funcForCell(cell).swcr_, 
-                             -1.0, 
-                             funcForCell(cell).krorw_, 
-                             funcForCell(cell).kromax_,
-                             0.0, 
-                             swl, sowcr, swl, swcr, sgl, krorw, kro, dummy);
-            } else if (oilGas) {
-                // ### krg
-                initEPSParam(cell, eps_transf_[cell].gas, false, 
-                             funcForCell(cell).smin_[gpos], 
-                             funcForCell(cell).sgcr_, 
-                             funcForCell(cell).smax_[gpos],
-                             funcForCell(cell).sogcr_, 
-                             -1.0, 
-                             funcForCell(cell).krgr_, 
-                             funcForCell(cell).krgmax_, 
-                             funcForCell(cell).pcgmax_, 
-                             sgl, sgcr, sgu, sogcr, swl, krgr, krg, pcg);
-                // ### krog
-                initEPSParam(cell, eps_transf_[cell].gasoil, true, 
-                             0.0, 
-                             funcForCell(cell).sogcr_, 
-                             funcForCell(cell).smin_[gpos],
-                             funcForCell(cell).sgcr_, 
-                             -1.0, 
-                             funcForCell(cell).krorg_, 
-                             funcForCell(cell).kromax_,
-                             0.0, 
-                             sgl, sogcr, sgl, sgcr, swl, krorg, kro, dummy);
-            } else if (threephase) {
-                // ### krw
-                initEPSParam(cell, eps_transf_[cell].wat, false, 
-                             funcForCell(cell).smin_[wpos], 
-                             funcForCell(cell).swcr_, 
-                             funcForCell(cell).smax_[wpos], 
-                             funcForCell(cell).sowcr_,
-                             funcForCell(cell).smin_[gpos], 
-                             funcForCell(cell).krwr_, 
-                             funcForCell(cell).krwmax_, 
-                             funcForCell(cell).pcwmax_, 
-                             swl, swcr, swu, sowcr, sgl, krwr, krw, pcw);
-                 // ### krow
-                initEPSParam(cell, eps_transf_[cell].watoil, true, 
-                            0.0, 
-                            funcForCell(cell).sowcr_, 
-                            funcForCell(cell).smin_[wpos], 
-                            funcForCell(cell).swcr_,
-                            funcForCell(cell).smin_[gpos], 
-                            funcForCell(cell).krorw_, 
-                            funcForCell(cell).kromax_,
-                            0.0, 
-                            swl, sowcr, swl, swcr, sgl, krorw, kro, dummy);
-                 // ### krg
-                initEPSParam(cell, eps_transf_[cell].gas, false, 
-                             funcForCell(cell).smin_[gpos], 
-                             funcForCell(cell).sgcr_, 
-                             funcForCell(cell).smax_[gpos], 
-                             funcForCell(cell).sogcr_,
-                             funcForCell(cell).smin_[wpos], 
-                             funcForCell(cell).krgr_, 
-                             funcForCell(cell).krgmax_, 
-                             funcForCell(cell).pcgmax_, 
-                             sgl, sgcr, sgu, sogcr, swl, krgr, krg, pcg);
-                 // ### krog
-                initEPSParam(cell, eps_transf_[cell].gasoil, true, 
-                             0.0, 
-                             funcForCell(cell).sogcr_, 
-                             funcForCell(cell).smin_[gpos], 
-                             funcForCell(cell).sgcr_,
-                             funcForCell(cell).smin_[wpos], 
-                             funcForCell(cell).krorg_, 
-                             funcForCell(cell).kromax_,
-                             0.0, 
-                             sgl, sogcr, sgl, sgcr, swl, krorg, kro, dummy);
-            }
+        
+        for (size_t i = 0; i < eps_kw.size(); ++i) {
+            initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                       eps_kw[i], eps_vec[i]);
         }
-    }
-
-    // Initialize hysteresis saturation scaling parameters
-    template <class SatFuncSet>
-    template<class T>
-    void SaturationPropsFromDeck<SatFuncSet>::initEPSHyst(Opm::DeckConstPtr deck,
-                                                      Opm::EclipseStateConstPtr eclipseState,
-                                                      int number_of_cells,
-                                                      const int* global_cell,
-                                                      const T& begin_cell_centroid,
-                                                      int dimensions)
-    {
-        std::vector<double> iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr;
-        std::vector<double> ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg;
-        std::vector<double> ipcw, ipcg;
-        const std::vector<double> dummy;
-        // Initialize hysteresis saturation scaling parameters
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("ISWL"),   iswl);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("ISWU"),   iswu);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("ISWCR"),  iswcr);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("ISGL"),   isgl);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("ISGU"),   isgu);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("ISGCR"),  isgcr);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("ISOWCR"), isowcr);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("ISOGCR"), isogcr);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("IKRW"),   ikrw);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("IKRG"),   ikrg);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("IKRO"),   ikro);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("IKRWR"),  ikrwr);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("IKRGR"),  ikrgr);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("IKRORW"), ikrorw);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("IKRORG"), ikrorg);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("IPCW"), ipcw);
-        initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
-                   std::string("IPCG"), ipcg);
-
-        eps_transf_hyst_.resize(number_of_cells);
-        sat_hyst_.resize(number_of_cells);
 
         const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
         const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
@@ -668,97 +490,53 @@ namespace Opm
         const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
 
         for (int cell = 0; cell < number_of_cells; ++cell) {
-            if (oilWater) {
-                 // ### krw
-                 initEPSParam(cell, eps_transf_hyst_[cell].wat, false, 
-                              funcForCell(cell).smin_[wpos], 
-                              funcForCell(cell).swcr_, 
-                              funcForCell(cell).smax_[wpos],
-                              funcForCell(cell).sowcr_, 
-                              -1.0, 
-                              funcForCell(cell).krwr_, 
-                              funcForCell(cell).krwmax_, 
-                              funcForCell(cell).pcwmax_, 
-                              iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw, ipcw);
-                 // ### krow
-                 initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 
-                              0.0, 
-                              funcForCell(cell).sowcr_, 
-                              funcForCell(cell).smin_[wpos],
-                              funcForCell(cell).swcr_, 
-                              -1.0, 
-                              funcForCell(cell).krorw_, 
-                              funcForCell(cell).kromax_,
-                              0.0, 
-                              iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro, dummy);
-             } else if (oilGas) {
-                 // ### krg
-                 initEPSParam(cell, eps_transf_hyst_[cell].gas, false, 
-                              funcForCell(cell).smin_[gpos], 
-                              funcForCell(cell).sgcr_, 
-                              funcForCell(cell).smax_[gpos],
-                              funcForCell(cell).sogcr_, 
-                              -1.0, 
-                              funcForCell(cell).krgr_, 
-                              funcForCell(cell).krgmax_, 
-                              funcForCell(cell).pcgmax_, 
-                              isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg, ipcg);
-                 // ### krog
-                 initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 
-                              0.0, 
-                              funcForCell(cell).sogcr_, 
-                              funcForCell(cell).smin_[gpos],
-                              funcForCell(cell).sgcr_, 
-                              -1.0, 
-                              funcForCell(cell).krorg_, 
-                              funcForCell(cell).kromax_,
-                              0.0, 
-                              isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro, dummy);
-             } else if (threephase) {
-                 // ### krw
-                 initEPSParam(cell, eps_transf_hyst_[cell].wat, false, 
-                              funcForCell(cell).smin_[wpos], 
-                              funcForCell(cell).swcr_, 
-                              funcForCell(cell).smax_[wpos], 
-                              funcForCell(cell).sowcr_,
-                              funcForCell(cell).smin_[gpos], 
-                              funcForCell(cell).krwr_, 
-                              funcForCell(cell).krwmax_, 
-                              funcForCell(cell).pcwmax_, 
-                              iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw, ipcw);
-                 // ### krow
-                 initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 
-                              0.0, 
-                              funcForCell(cell).sowcr_, 
-                              funcForCell(cell).smin_[wpos], 
-                              funcForCell(cell).swcr_,
-                              funcForCell(cell).smin_[gpos], 
-                              funcForCell(cell).krorw_, 
-                              funcForCell(cell).kromax_,
-                              0.0, 
-                              iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro, dummy);
-                 // ### krg
-                 initEPSParam(cell, eps_transf_hyst_[cell].gas, false, 
-                              funcForCell(cell).smin_[gpos], 
-                              funcForCell(cell).sgcr_, 
-                              funcForCell(cell).smax_[gpos], 
-                              funcForCell(cell).sogcr_,
-                              funcForCell(cell).smin_[wpos], 
-                              funcForCell(cell).krgr_, 
-                              funcForCell(cell).krgmax_, 
-                              funcForCell(cell).pcgmax_, 
-                              isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg, ipcg);
-                 // ### krog
-                 initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 
-                              0.0, 
-                              funcForCell(cell).sogcr_, 
-                              funcForCell(cell).smin_[gpos], 
-                              funcForCell(cell).sgcr_, 
-                              funcForCell(cell).smin_[wpos], 
-                              funcForCell(cell).krorg_, 
-                              funcForCell(cell).kromax_,
-                              0.0, 
-                              isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro, dummy);
+            if (threephase || oilWater) {
+                // ### krw
+                initEPSParam(cell, eps_transf[cell].wat, false,
+                             funcForCell(cell).smin_[wpos],
+                             funcForCell(cell).swcr_,
+                             funcForCell(cell).smax_[wpos],
+                             funcForCell(cell).sowcr_,
+                             oilWater ? -1.0 : funcForCell(cell).smin_[gpos],
+                             funcForCell(cell).krwr_,
+                             funcForCell(cell).krwmax_,
+                             funcForCell(cell).pcwmax_,
+                             eps_vec[0], eps_vec[2], eps_vec[1], eps_vec[6], eps_vec[3], eps_vec[11], eps_vec[8], eps_vec[15]);
+                // ### krow
+                initEPSParam(cell, eps_transf[cell].watoil, true,
+                             0.0,
+                             funcForCell(cell).sowcr_,
+                             funcForCell(cell).smin_[wpos],
+                             funcForCell(cell).swcr_,
+                             oilWater ? -1.0 : funcForCell(cell).smin_[gpos],
+                             funcForCell(cell).krorw_,
+                             funcForCell(cell).kromax_,
+                             0.0,
+                             eps_vec[0], eps_vec[6], eps_vec[0], eps_vec[2], eps_vec[3], eps_vec[13], eps_vec[10], dummy);
+            }
+            if (threephase || oilGas) {
+                // ### krg
+                initEPSParam(cell, eps_transf[cell].gas, false,
+                             funcForCell(cell).smin_[gpos],
+                             funcForCell(cell).sgcr_,
+                             funcForCell(cell).smax_[gpos],
+                             funcForCell(cell).sogcr_,
+                             oilGas ? -1.0 : funcForCell(cell).smin_[wpos],
+                             funcForCell(cell).krgr_,
+                             funcForCell(cell).krgmax_,
+                             funcForCell(cell).pcgmax_,
+                             eps_vec[3], eps_vec[5], eps_vec[4], eps_vec[7], eps_vec[0], eps_vec[12], eps_vec[9], eps_vec[16]);
+                // ### krog
+                initEPSParam(cell, eps_transf[cell].gasoil, true,
+                             0.0,
+                             funcForCell(cell).sogcr_,
+                             funcForCell(cell).smin_[gpos],
+                             funcForCell(cell).sgcr_,
+                             oilGas ? -1.0 : funcForCell(cell).smin_[wpos],
+                             funcForCell(cell).krorg_,
+                             funcForCell(cell).kromax_,
+                             0.0,
+                             eps_vec[3], eps_vec[7], eps_vec[3], eps_vec[5], eps_vec[0], eps_vec[14], eps_vec[10], dummy);
             }
         }
     }
@@ -779,6 +557,12 @@ namespace Opm
         const bool useLiquid = phase_usage_.phase_used[Liquid];
         const bool useVapour = phase_usage_.phase_used[Vapour];
         bool useKeyword = deck->hasKeyword(keyword);
+        bool useStateKeyword = eclipseState->hasDoubleGridProperty(keyword);
+        const std::map<std::string, int> kw2tab = {
+            {"SWL", 1}, {"SWCR", 2}, {"SWU", 3}, {"SGL", 4},
+            {"SGCR", 5}, {"SGU", 6}, {"SOWCR", 7}, {"SOGCR", 8},
+            {"ISWL", 1}, {"ISWCR", 2}, {"ISWU", 3}, {"ISGL", 4},
+            {"ISGCR", 5}, {"ISGU", 6}, {"ISOWCR", 7}, {"ISOGCR", 8}};
         bool hasENPTVD = deck->hasKeyword("ENPTVD");
         bool hasENKRVD = deck->hasKeyword("ENKRVD");
         int itab = 0;
@@ -787,67 +571,10 @@ namespace Opm
         std::vector<std::string> col_names;
 
         // Active keyword assigned default values for each cell (in case of possible box-wise assignment)
-        int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua];
-        int phase_pos_vapour = phase_usage_.phase_pos[BlackoilPhases::Vapour];
-        if ((keyword[0] == 'S' && (useKeyword || hasENPTVD)) || (keyword[1] == 'S' && useKeyword) ) {
-            if (keyword == std::string("SWL") || keyword == std::string("ISWL") ) {
-                if (useAqua && (useKeyword || columnIsMasked_(deck, "ENPTVD", 0))) {
-                    itab = 1;
-                    scaleparam.resize(number_of_cells);
-                    for (int i=0; i<number_of_cells; ++i)
-                        scaleparam[i] = funcForCell(i).smin_[phase_pos_aqua];
-                }
-            } else if (keyword == std::string("SWCR") || keyword == std::string("ISWCR") ) {
-                if (useAqua && (useKeyword || columnIsMasked_(deck, "ENPTVD", 1))) {
-                    itab = 2;
-                    scaleparam.resize(number_of_cells);
-                    for (int i=0; i<number_of_cells; ++i)
-                        scaleparam[i] = funcForCell(i).swcr_;
-                }
-            } else if (keyword == std::string("SWU") || keyword == std::string("ISWU") ) {
-                if (useAqua && (useKeyword || columnIsMasked_(deck, "ENPTVD", 2))) {
-                    itab = 3;
-                    scaleparam.resize(number_of_cells);
-                    for (int i=0; i<number_of_cells; ++i)
-                        scaleparam[i] = funcForCell(i).smax_[phase_pos_aqua];
-                }
-            } else if (keyword == std::string("SGL") || keyword == std::string("ISGL") ) {
-                if (useVapour && (useKeyword || columnIsMasked_(deck, "ENPTVD", 3))) {
-                    itab = 4;
-                    scaleparam.resize(number_of_cells);
-                    for (int i=0; i<number_of_cells; ++i)
-                        scaleparam[i] = funcForCell(i).smin_[phase_pos_vapour];
-                }
-            } else if (keyword == std::string("SGCR") || keyword == std::string("ISGCR") ) {
-                if (useVapour && (useKeyword || columnIsMasked_(deck, "ENPTVD", 4))) {
-                    itab = 5;
-                    scaleparam.resize(number_of_cells);
-                    for (int i=0; i<number_of_cells; ++i)
-                        scaleparam[i] = funcForCell(i).sgcr_;
-                }
-            } else if (keyword == std::string("SGU") || keyword == std::string("ISGU") ) {
-                if (useVapour && (useKeyword || columnIsMasked_(deck, "ENPTVD", 5))) {
-                    itab = 6;
-                    scaleparam.resize(number_of_cells);
-                    for (int i=0; i<number_of_cells; ++i)
-                        scaleparam[i] = funcForCell(i).smax_[phase_pos_vapour];
-                }
-            } else if (keyword == std::string("SOWCR") || keyword == std::string("ISOWCR") ) {
-                if (useAqua && (useKeyword || columnIsMasked_(deck, "ENPTVD", 6))) {
-                    itab = 7;
-                    scaleparam.resize(number_of_cells);
-                    for (int i=0; i<number_of_cells; ++i)
-                        scaleparam[i] = funcForCell(i).sowcr_;
-                }
-            } else if (keyword == std::string("SOGCR") || keyword == std::string("ISOGCR") ) {
-                if (useVapour && (useKeyword || columnIsMasked_(deck, "ENPTVD", 7))) {
-                    itab = 8;
-                    scaleparam.resize(number_of_cells);
-                    for (int i=0; i<number_of_cells; ++i)
-                        scaleparam[i] = funcForCell(i).sogcr_;
-                }
-            } else {
-                OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
+        if ((keyword[0] == 'S' && (useStateKeyword || hasENPTVD)) || (keyword[1] == 'S' && useStateKeyword) ) {
+            if (useAqua && (useStateKeyword || columnIsMasked_(deck, "ENPTVD", kw2tab.find(keyword)->second-1))) {
+                itab = kw2tab.find(keyword)->second;
+                scaleparam.resize(number_of_cells);
             }
             if (!useKeyword && itab > 0) {
                 const auto& enptvdTables = eclipseState->getEnptvdTables();
@@ -926,7 +653,7 @@ namespace Opm
                     param_col[table_num] = enkrvdTable.getColumn(itab); // itab=[1-7]: krw krg kro krwr krgr krorw krorg
                 }
             }
-         } else if (useKeyword && (keyword[0] == 'P' || keyword[1] == 'P') ) {
+        } else if (useKeyword && (keyword[0] == 'P' || keyword[1] == 'P') ) {
              if (useAqua && (keyword == std::string("PCW") || keyword == std::string("IPCW")) ) {
                  scaleparam.resize(number_of_cells);
                  for (int i=0; i<number_of_cells; ++i)
@@ -940,16 +667,25 @@ namespace Opm
 
         if (scaleparam.empty()) {
             return;
-        } else if (useKeyword) {
+        }
+
+        if (useKeyword || useStateKeyword) {
             // Keyword values from deck
             std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl;
             const int* gc = global_cell;
-            const std::vector<double>& val = deck->getKeyword(keyword)->getSIDoubleData();
+            std::vector<double> val;
+            if (keyword[0] == 'S' || keyword[1] == 'S') { // Saturation from EclipseState
+                val = eclipseState->getDoubleGridProperty(keyword)->getData();
+            } else {
+                val = deck->getKeyword(keyword)->getSIDoubleData(); //KR and PC directly from deck.
+            }
             for (int c = 0; c < int(scaleparam.size()); ++c) {
                 const int deck_pos = (gc == NULL) ? c : gc[c];
                 scaleparam[c] = val[deck_pos];
             }
-        } else {
+        }
+
+        if (itab > 0) {
             const int dim = dimensions;
             std::vector<int> endnum;
             if ( deck->hasKeyword("ENDNUM")) {
@@ -966,17 +702,39 @@ namespace Opm
                 // Default deck value is one
                 endnum.assign(number_of_cells, 0);
             }
-            for (int cell = 0; cell < number_of_cells; ++cell) {
-                if (endnum[cell] >= 0 && param_col[endnum[cell]][0] >= 0.0) {
-                    double zc = UgGridHelpers
-                        ::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim),
-                                       dim-1);
-                    if (zc >= depth_col[endnum[cell]].front() && zc <= depth_col[endnum[cell]].back()) { //don't want extrap outside depth interval
-                        scaleparam[cell] = linearInterpolation(depth_col[endnum[cell]], param_col[endnum[cell]], zc);
+            if (keyword[0] == 'S' || keyword[1] == 'S') { // From EclipseState
+                for (int cell = 0; cell < number_of_cells; ++cell) {
+                    if (!std::isfinite(scaleparam[cell]) && endnum[cell] >= 0 && param_col[endnum[cell]][0] >= 0.0) {
+                        double zc = UgGridHelpers
+                            ::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim),
+                                           dim-1);
+                        if (zc >= depth_col[endnum[cell]].front() && zc <= depth_col[endnum[cell]].back()) { //don't want extrap outside depth interval
+                            scaleparam[cell] = linearInterpolation(depth_col[endnum[cell]], param_col[endnum[cell]], zc);
+                        }
+                    } else if (!std::isfinite(scaleparam[cell]) && endnum[cell] >= 0) {
+                        // As of 1/9-2014:  Reflects remaining work on opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp ...
+                        OPM_THROW(std::runtime_error, " -- Inconsistent EclipseState: '" << keyword << "' (ENPTVD)");
+                    }
+                }
+            } else { //KR and PC from deck.
+                for (int cell = 0; cell < number_of_cells; ++cell) {
+                    if (endnum[cell] >= 0 && param_col[endnum[cell]][0] >= 0.0) {
+                        double zc = UgGridHelpers
+                            ::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim),
+                                           dim-1);
+                        if (zc >= depth_col[endnum[cell]].front() && zc <= depth_col[endnum[cell]].back()) { //don't want extrap outside depth interval
+                            scaleparam[cell] = linearInterpolation(depth_col[endnum[cell]], param_col[endnum[cell]], zc);
+                        }
                     }
                 }
             }
         }
+
+//        std::cout << keyword << ":" << std::endl;
+//        for (int c = 0; c < int(scaleparam.size()); ++c) {
+//                std::cout << c << "    " << scaleparam[c] << std::endl;
+//        }
+
     }
 
     // Saturation scaling
diff --git a/tests/satfuncEPSBase.DATA b/tests/satfuncEPSBase.DATA
new file mode 100644
index 00000000..5f47737c
--- /dev/null
+++ b/tests/satfuncEPSBase.DATA
@@ -0,0 +1,168 @@
+NOECHO
+
+RUNSPEC   ======
+
+WATER
+OIL
+GAS
+DISGAS
+VAPOIL
+
+TABDIMS
+  1    1   40   20    1   20  /
+
+DIMENS
+1 1 20
+/
+
+WELLDIMS
+   30   10    2   30 /
+   
+ENDSCALE
+--DIR      REV      NTENDP    NSENDP
+'NODIR'  'REVERS'    1          20   /
+/ 
+
+START
+   1 'JAN' 1990  /
+
+NSTACK
+   25 /
+
+EQLDIMS
+-- NTEQUL
+     1 / 
+     
+
+FMTOUT
+FMTIN
+
+GRID      ======
+
+DXV
+1.0
+/
+
+DYV
+1.0
+/
+
+DZV
+20*5.0
+/
+
+
+PORO
+20*0.2
+/
+
+
+PERMZ
+  20*1.0
+/
+
+PERMY
+20*100.0
+/
+
+PERMX
+20*100.0
+/
+
+BOX
+ 1 1 1 1 1 1 /
+
+TOPS
+0.0
+/
+
+PROPS     ======
+
+PVTO
+--     Rs       Pbub       Bo        Vo
+         0          1.    1.0000     1.20  /
+        20         40.    1.0120     1.17  /
+        40         80.    1.0255     1.14  /
+        60        120.    1.0380     1.11  /
+        80        160.    1.0510     1.08  /
+       100        200.    1.0630     1.06  /
+       120        240.    1.0750     1.03  /
+       140        280.    1.0870     1.00  /
+       160        320.    1.0985      .98  /
+       180        360.    1.1100      .95  /
+       200        400.    1.1200      .94
+                  500.    1.1189      .94  /
+/
+
+PVTG
+--  Pg     Rv        Bg       Vg
+   100   0.0001       0.010      0.1
+         0.0          0.0104     0.1 /
+   200   0.0004       0.005      0.2
+         0.0          0.0054     0.2 /
+/
+
+SCALECRS
+--  YES /
+ NO/
+
+SWOF
+0.1 0.0 1.0 0.9
+0.2 0.0 0.8 0.8
+0.3 0.1 0.6 0.7
+0.4 0.2 0.4 0.6
+0.7 0.5 0.1 0.3
+0.8 0.6 0.0 0.2
+0.9 0.7 0.0 0.1
+/
+
+SGOF
+0.0 0.0 1.0 0.2
+0.1 0.0 0.7 0.4
+0.2 0.1 0.6 0.6
+0.8 0.7 0.0 2.0
+0.9 1.0 0.0 2.1
+/
+
+PVTW
+--RefPres  Bw      Comp   Vw    Cv
+   1.      1.0   4.0E-5  0.96  0.0 /
+   
+
+ROCK
+--RefPres  Comp
+   1.   5.0E-5 /
+
+DENSITY
+700 1000 1
+/
+
+SOLUTION  ======
+
+EQUIL
+45 150 50 0.25 45 0.35 1 1 0
+/
+
+RSVD
+ 0  0.0
+ 100 100. /
+ 
+RVVD
+   0.  0.
+ 100.  0.0001 /
+
+RPTSOL
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
+
+SUMMARY   ======
+RUNSUM
+
+SEPARATE
+
+SCHEDULE  ======
+
+RPTSCHED
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
+
+
+END
diff --git a/tests/satfuncEPS_A.DATA b/tests/satfuncEPS_A.DATA
new file mode 100644
index 00000000..47e9612a
--- /dev/null
+++ b/tests/satfuncEPS_A.DATA
@@ -0,0 +1,177 @@
+NOECHO
+
+RUNSPEC   ======
+
+WATER
+OIL
+GAS
+DISGAS
+VAPOIL
+
+TABDIMS
+  1    1   40   20    1   20  /
+
+DIMENS
+1 1 10
+/
+
+WELLDIMS
+   30   10    2   30 /
+   
+ENDSCALE
+--DIR      REV      NTENDP    NSENDP
+'NODIR'  'REVERS'    1          20   /
+/ 
+
+START
+   1 'JAN' 1990  /
+
+NSTACK
+   25 /
+
+EQLDIMS
+-- NTEQUL
+     1 / 
+     
+
+FMTOUT
+FMTIN
+
+GRID      ======
+
+DXV
+1.0
+/
+
+DYV
+1.0
+/
+
+DZV
+10*5.0
+/
+
+
+PORO
+10*0.2
+/
+
+
+PERMZ
+  10*1.0
+/
+
+PERMY
+10*100.0
+/
+
+PERMX
+10*100.0
+/
+
+BOX
+ 1 1 1 1 1 1 /
+
+TOPS
+0.0
+/
+
+PROPS     ======
+
+PVTO
+--     Rs       Pbub       Bo        Vo
+         0          1.    1.0000     1.20  /
+        20         40.    1.0120     1.17  /
+        40         80.    1.0255     1.14  /
+        60        120.    1.0380     1.11  /
+        80        160.    1.0510     1.08  /
+       100        200.    1.0630     1.06  /
+       120        240.    1.0750     1.03  /
+       140        280.    1.0870     1.00  /
+       160        320.    1.0985      .98  /
+       180        360.    1.1100      .95  /
+       200        400.    1.1200      .94
+                  500.    1.1189      .94  /
+/
+
+PVTG
+--  Pg     Rv        Bg       Vg
+   100   0.0001       0.010      0.1
+         0.0          0.0104     0.1 /
+   200   0.0004       0.005      0.2
+         0.0          0.0054     0.2 /
+/
+
+SCALECRS
+--  YES /
+ NO/
+ 
+SWL
+4*0.1 4*0.2 2*0.1/
+ 
+SWCR
+0.2 0.2 0.4 0.4 0.2 0.2 0.4 0.4 2*0.2 /
+
+SWU
+0.9 0.7 0.9 0.7 0.9 0.7 0.9 0.7 2*0.9 /
+
+SWOF
+0.1 0.0 1.0 0.9
+0.2 0.0 0.8 0.8
+0.3 0.1 0.6 0.7
+0.4 0.2 0.4 0.6
+0.7 0.5 0.1 0.3
+0.8 0.6 0.0 0.2
+0.9 0.7 0.0 0.1
+/
+
+SGOF
+0.0 0.0 1.0 0.2
+0.1 0.0 0.7 0.4
+0.2 0.1 0.6 0.6
+0.8 0.7 0.0 2.0
+0.9 1.0 0.0 2.1
+/
+
+PVTW
+--RefPres  Bw      Comp   Vw    Cv
+   1.      1.0   4.0E-5  0.96  0.0 /
+   
+
+ROCK
+--RefPres  Comp
+   1.   5.0E-5 /
+
+DENSITY
+700 1000 1
+/
+
+SOLUTION  ======
+
+EQUIL
+45 150 50 0.25 45 0.35 1 1 0
+/
+
+RSVD
+ 0  0.0
+ 100 100. /
+ 
+RVVD
+   0.  0.
+ 100.  0.0001 /
+
+RPTSOL
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
+
+SUMMARY   ======
+RUNSUM
+
+SEPARATE
+
+SCHEDULE  ======
+
+RPTSCHED
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
+
+
+END
diff --git a/tests/satfuncEPS_B.DATA b/tests/satfuncEPS_B.DATA
new file mode 100644
index 00000000..a1d8a317
--- /dev/null
+++ b/tests/satfuncEPS_B.DATA
@@ -0,0 +1,190 @@
+NOECHO
+
+RUNSPEC   ======
+
+WATER
+OIL
+GAS
+DISGAS
+VAPOIL
+
+TABDIMS
+  1    1   40   20    1   20  /
+
+DIMENS
+1 1 10
+/
+
+WELLDIMS
+   30   10    2   30 /
+   
+ENDSCALE
+--DIR      REV      NTENDP    NSENDP
+'NODIR'  'REVERS'    1          20   /
+/ 
+
+START
+   1 'JAN' 1990  /
+
+NSTACK
+   25 /
+
+EQLDIMS
+-- NTEQUL
+     1 / 
+     
+
+FMTOUT
+FMTIN
+
+GRID      ======
+
+DXV
+1.0
+/
+
+DYV
+1.0
+/
+
+DZV
+10*5.0
+/
+
+
+PORO
+10*0.2
+/
+
+
+PERMZ
+  10*1.0
+/
+
+PERMY
+10*100.0
+/
+
+PERMX
+10*100.0
+/
+
+BOX
+ 1 1 1 1 1 1 /
+
+TOPS
+0.0
+/
+
+PROPS     ======
+
+PVTO
+--     Rs       Pbub       Bo        Vo
+         0          1.    1.0000     1.20  /
+        20         40.    1.0120     1.17  /
+        40         80.    1.0255     1.14  /
+        60        120.    1.0380     1.11  /
+        80        160.    1.0510     1.08  /
+       100        200.    1.0630     1.06  /
+       120        240.    1.0750     1.03  /
+       140        280.    1.0870     1.00  /
+       160        320.    1.0985      .98  /
+       180        360.    1.1100      .95  /
+       200        400.    1.1200      .94
+                  500.    1.1189      .94  /
+/
+
+PVTG
+--  Pg     Rv        Bg       Vg
+   100   0.0001       0.010      0.1
+         0.0          0.0104     0.1 /
+   200   0.0004       0.005      0.2
+         0.0          0.0054     0.2 /
+/
+
+SCALECRS
+--  YES /
+ NO/
+
+ENPTVD 
+-- Depth   swl   swcr   swu   sgl   sgcr   sgu  sowcr  sogcr 
+-- Sandstone 
+   2.5     .1    0.2    0.9   .0     0.1   0.9   0.2    0.2
+   7.5     .1    0.2    0.7   .0     0.1   0.9   0.2    0.2 
+  12.5     .1    0.4    0.9   .0     0.1   0.9   0.2    0.2
+  17.5     .1    0.4    0.7   .0     0.1   0.9   0.2    0.2
+  22.5     .2    0.2    0.9   .0     0.1   0.9   0.2    0.2
+  27.5     .2    0.2    0.7   .0     0.1   0.9   0.2    0.2
+  32.5     .2    0.4    0.9   .0     0.1   0.9   0.2    0.2
+  37.5     .2    0.4    0.7   .0     0.1   0.9   0.2    0.2
+  42.5     .1    0.2    0.9   .0     0.1   0.9   0.2    0.2
+  47.5     .1    0.2    0.9   .0     0.1   0.9   0.2    0.2
+  50.0     .1    0.2    0.9   .0     0.1   0.9   0.2    0.2 /
+/ 
+
+SWOF
+0.1 0.0 1.0 0.9
+0.2 0.0 0.8 0.8
+0.3 0.1 0.6 0.7
+0.4 0.2 0.4 0.6
+0.7 0.5 0.1 0.3
+0.8 0.6 0.0 0.2
+0.9 0.7 0.0 0.1
+/
+
+SGOF
+0.0 0.0 1.0 0.2
+0.1 0.0 0.7 0.4
+0.2 0.1 0.6 0.6
+0.8 0.7 0.0 2.0
+0.9 1.0 0.0 2.1
+/
+
+PVTW
+--RefPres  Bw      Comp   Vw    Cv
+   1.      1.0   4.0E-5  0.96  0.0 /
+   
+
+ROCK
+--RefPres  Comp
+   1.   5.0E-5 /
+
+DENSITY
+700 1000 1
+/
+
+REGIONS   ======
+
+ENDNUM
+10*1
+/
+
+SOLUTION  ======
+
+EQUIL
+45 150 50 0.25 45 0.35 1 1 0
+/
+
+RSVD
+ 0  0.0
+ 100 100. /
+ 
+RVVD
+   0.  0.
+ 100.  0.0001 /
+
+RPTSOL
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
+
+SUMMARY   ======
+RUNSUM
+
+SEPARATE
+
+SCHEDULE  ======
+
+RPTSCHED
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
+
+
+END
diff --git a/tests/satfuncEPS_C.DATA b/tests/satfuncEPS_C.DATA
new file mode 100644
index 00000000..74a89c02
--- /dev/null
+++ b/tests/satfuncEPS_C.DATA
@@ -0,0 +1,225 @@
+NOECHO
+
+RUNSPEC   ======
+
+WATER
+OIL
+GAS
+DISGAS
+VAPOIL
+
+TABDIMS
+  1    1   40   20    1   20  /
+
+DIMENS
+1 1 10
+/
+
+WELLDIMS
+   30   10    2   30 /
+   
+ENDSCALE
+--DIR      REV      NTENDP    NSENDP
+'NODIR'  'REVERS'    1          20   /
+/ 
+
+START
+   1 'JAN' 1990  /
+
+NSTACK
+   25 /
+
+EQLDIMS
+-- NTEQUL
+     1 / 
+     
+
+FMTOUT
+FMTIN
+
+GRID      ======
+
+DXV
+1.0
+/
+
+DYV
+1.0
+/
+
+DZV
+10*5.0
+/
+
+
+PORO
+10*0.2
+/
+
+
+PERMZ
+  10*1.0
+/
+
+PERMY
+10*100.0
+/
+
+PERMX
+10*100.0
+/
+
+BOX
+ 1 1 1 1 1 1 /
+
+TOPS
+0.0
+/
+
+PROPS     ======
+
+PVTO
+--     Rs       Pbub       Bo        Vo
+         0          1.    1.0000     1.20  /
+        20         40.    1.0120     1.17  /
+        40         80.    1.0255     1.14  /
+        60        120.    1.0380     1.11  /
+        80        160.    1.0510     1.08  /
+       100        200.    1.0630     1.06  /
+       120        240.    1.0750     1.03  /
+       140        280.    1.0870     1.00  /
+       160        320.    1.0985      .98  /
+       180        360.    1.1100      .95  /
+       200        400.    1.1200      .94
+                  500.    1.1189      .94  /
+/
+
+PVTG
+--  Pg     Rv        Bg       Vg
+   100   0.0001       0.010      0.1
+         0.0          0.0104     0.1 /
+   200   0.0004       0.005      0.2
+         0.0          0.0054     0.2 /
+/
+
+SCALECRS
+--  YES /
+ NO/
+ 
+--SWL
+--4*0.1 4*0.2 2*0.1/
+
+EQUALS
+SWL 0.1 1 1 1 1 1 1  /
+SWL 0.1 1 1 1 1 2 2  /
+SWL 0.1 1 1 1 1 3 3  /
+SWL 0.1 1 1 1 1 4 4  /
+SWL 0.2 1 1 1 1 5 5  /
+SWL 0.2 1 1 1 1 6 6  /
+SWL 0.2 1 1 1 1 7 7  /
+SWL 0.2 1 1 1 1 8 8  /
+SWL 0.1 1 1 1 1 9 9  /
+SWL 0.1 1 1 1 1 10 10  /
+/
+ 
+--SWCR
+--0.2 0.2 0.4 0.4 0.2 0.2 0.4 0.4 2*0.2 /
+
+COPY
+  SWL SWCR /
+  SWL SWU  /
+/
+
+ADD
+  SWCR 0.1 1 1 1 1 1 2  /
+  SWCR 0.3 1 1 1 1 3 4  /
+  SWCR 0.2 1 1 1 1 7 8  /
+  SWCR 0.1 1 1 1 1 9 10 /
+/
+
+--SWU
+--0.9 0.7 0.9 0.7 0.9 0.7 0.9 0.7 2*0.9 /
+
+MULTIPLY
+  SWU -1.0 1 1 1 1 1 10 /
+/
+
+ADD
+  SWU  1.0 1 1 1 1 1 10 /
+/
+
+ADD
+  SWU -0.2 1 1 1 1 2 2 /
+  SWU -0.2 1 1 1 1 4 4 /
+  SWU  0.1 1 1 1 1 5 5 /
+  SWU -0.1 1 1 1 1 6 6 /
+  SWU  0.1 1 1 1 1 7 7 /
+  SWU -0.1 1 1 1 1 8 8 /
+/
+
+SWOF
+0.1 0.0 1.0 0.9
+0.2 0.0 0.8 0.8
+0.3 0.1 0.6 0.7
+0.4 0.2 0.4 0.6
+0.7 0.5 0.1 0.3
+0.8 0.6 0.0 0.2
+0.9 0.7 0.0 0.1
+/
+
+SGOF
+0.0 0.0 1.0 0.2
+0.1 0.0 0.7 0.4
+0.2 0.1 0.6 0.6
+0.8 0.7 0.0 2.0
+0.9 1.0 0.0 2.1
+/
+
+PVTW
+--RefPres  Bw      Comp   Vw    Cv
+   1.      1.0   4.0E-5  0.96  0.0 /
+   
+
+ROCK
+--RefPres  Comp
+   1.   5.0E-5 /
+
+DENSITY
+700 1000 1
+/
+
+REGIONS   ======
+
+ENDNUM
+10*1
+/
+
+SOLUTION  ======
+
+EQUIL
+45 150 50 0.25 45 0.35 1 1 0
+/
+
+RSVD
+ 0  0.0
+ 100 100. /
+ 
+RVVD
+   0.  0.
+ 100.  0.0001 /
+
+RPTSOL
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
+
+SUMMARY   ======
+RUNSUM
+
+SEPARATE
+
+SCHEDULE  ======
+
+RPTSCHED
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
+
+
+END
diff --git a/tests/satfuncEPS_D.DATA b/tests/satfuncEPS_D.DATA
new file mode 100644
index 00000000..1077121e
--- /dev/null
+++ b/tests/satfuncEPS_D.DATA
@@ -0,0 +1,181 @@
+NOECHO
+
+RUNSPEC   ======
+
+WATER
+OIL
+GAS
+DISGAS
+VAPOIL
+
+TABDIMS
+  1    1   40   20    1   20  /
+
+DIMENS
+1 1 10
+/
+
+WELLDIMS
+   30   10    2   30 /
+   
+ENDSCALE
+--DIR      REV      NTENDP    NSENDP
+'NODIR'  'REVERS'    1          20   /
+/
+
+SATOPTS
+HYSTER  / 
+
+START
+   1 'JAN' 1990  /
+
+NSTACK
+   25 /
+
+EQLDIMS
+-- NTEQUL
+     1 / 
+     
+
+FMTOUT
+FMTIN
+
+GRID      ======
+
+DXV
+1.0
+/
+
+DYV
+1.0
+/
+
+DZV
+10*5.0
+/
+
+
+PORO
+10*0.2
+/
+
+
+PERMZ
+  10*1.0
+/
+
+PERMY
+10*100.0
+/
+
+PERMX
+10*100.0
+/
+
+BOX
+ 1 1 1 1 1 1 /
+
+TOPS
+0.0
+/
+
+PROPS     ======
+
+-- Hysteresis input
+EHYSTR
+   0.1   0  0.1 1* KR /
+
+PVTO
+--     Rs       Pbub       Bo        Vo
+         0          1.    1.0000     1.20  /
+        20         40.    1.0120     1.17  /
+        40         80.    1.0255     1.14  /
+        60        120.    1.0380     1.11  /
+        80        160.    1.0510     1.08  /
+       100        200.    1.0630     1.06  /
+       120        240.    1.0750     1.03  /
+       140        280.    1.0870     1.00  /
+       160        320.    1.0985      .98  /
+       180        360.    1.1100      .95  /
+       200        400.    1.1200      .94
+                  500.    1.1189      .94  /
+/
+
+PVTG
+--  Pg     Rv        Bg       Vg
+   100   0.0001       0.010      0.1
+         0.0          0.0104     0.1 /
+   200   0.0004       0.005      0.2
+         0.0          0.0054     0.2 /
+/
+
+SCALECRS
+--  YES /
+ NO/
+
+SWOF
+0.1 0.0 1.0 0.9
+0.2 0.0 0.8 0.8
+0.3 0.1 0.6 0.7
+0.4 0.2 0.4 0.6
+0.7 0.5 0.1 0.3
+0.8 0.6 0.0 0.2
+0.9 0.7 0.0 0.1
+/
+
+SGOF
+0.0 0.0 1.0 0.2
+0.1 0.0 0.7 0.4
+0.2 0.1 0.6 0.6
+0.8 0.7 0.0 2.0
+0.9 1.0 0.0 2.1
+/
+
+PVTW
+--RefPres  Bw      Comp   Vw    Cv
+   1.      1.0   4.0E-5  0.96  0.0 /
+   
+
+ROCK
+--RefPres  Comp
+   1.   5.0E-5 /
+
+DENSITY
+700 1000 1
+/
+
+REGIONS   ======
+
+ENDNUM
+10*1
+/
+
+SOLUTION  ======
+
+EQUIL
+45 150 50 0.25 45 0.35 1 1 0
+/
+
+RSVD
+ 0  0.0
+ 100 100. /
+ 
+RVVD
+   0.  0.
+ 100.  0.0001 /
+
+RPTSOL
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
+
+SUMMARY   ======
+RUNSUM
+
+SEPARATE
+
+SCHEDULE  ======
+
+RPTSCHED
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
+
+
+END
diff --git a/tests/satfuncStandard.DATA b/tests/satfuncStandard.DATA
new file mode 100644
index 00000000..d9173131
--- /dev/null
+++ b/tests/satfuncStandard.DATA
@@ -0,0 +1,168 @@
+NOECHO
+
+RUNSPEC   ======
+
+WATER
+OIL
+GAS
+DISGAS
+VAPOIL
+
+TABDIMS
+  1    1   40   20    1   20  /
+
+DIMENS
+1 1 10
+/
+
+WELLDIMS
+   30   10    2   30 /
+   
+--ENDSCALE
+--DIR      REV      NTENDP    NSENDP
+--'NODIR'  'REVERS'    1          20   /
+--/ 
+
+START
+   1 'JAN' 1990  /
+
+NSTACK
+   25 /
+
+EQLDIMS
+-- NTEQUL
+     1 / 
+     
+
+FMTOUT
+FMTIN
+
+GRID      ======
+
+DXV
+1.0
+/
+
+DYV
+1.0
+/
+
+DZV
+10*5.0
+/
+
+
+PORO
+10*0.2
+/
+
+
+PERMZ
+  10*1.0
+/
+
+PERMY
+10*100.0
+/
+
+PERMX
+10*100.0
+/
+
+BOX
+ 1 1 1 1 1 1 /
+
+TOPS
+0.0
+/
+
+PROPS     ======
+
+PVTO
+--     Rs       Pbub       Bo        Vo
+         0          1.    1.0000     1.20  /
+        20         40.    1.0120     1.17  /
+        40         80.    1.0255     1.14  /
+        60        120.    1.0380     1.11  /
+        80        160.    1.0510     1.08  /
+       100        200.    1.0630     1.06  /
+       120        240.    1.0750     1.03  /
+       140        280.    1.0870     1.00  /
+       160        320.    1.0985      .98  /
+       180        360.    1.1100      .95  /
+       200        400.    1.1200      .94
+                  500.    1.1189      .94  /
+/
+
+PVTG
+--  Pg     Rv        Bg       Vg
+   100   0.0001       0.010      0.1
+         0.0          0.0104     0.1 /
+   200   0.0004       0.005      0.2
+         0.0          0.0054     0.2 /
+/
+
+--SCALECRS
+--  YES /
+-- NO/
+
+SWOF
+0.1 0.0 1.0 0.9
+0.2 0.0 0.8 0.8
+0.3 0.1 0.6 0.7
+0.4 0.2 0.4 0.6
+0.7 0.5 0.1 0.3
+0.8 0.6 0.0 0.2
+0.9 0.7 0.0 0.1
+/
+
+SGOF
+0.0 0.0 1.0 0.2
+0.1 0.0 0.7 0.4
+0.2 0.1 0.6 0.6
+0.8 0.7 0.0 2.0
+0.9 1.0 0.0 2.1
+/
+
+PVTW
+--RefPres  Bw      Comp   Vw    Cv
+   1.      1.0   4.0E-5  0.96  0.0 /
+   
+
+ROCK
+--RefPres  Comp
+   1.   5.0E-5 /
+
+DENSITY
+700 1000 1
+/
+
+SOLUTION  ======
+
+EQUIL
+45 150 50 0.25 45 0.35 1 1 0
+/
+
+RSVD
+ 0  0.0
+ 100 100. /
+ 
+RVVD
+   0.  0.
+ 100.  0.0001 /
+
+RPTSOL
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
+
+SUMMARY   ======
+RUNSUM
+
+SEPARATE
+
+SCHEDULE  ======
+
+RPTSCHED
+'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
+
+
+END
diff --git a/tests/test_satfunc.cpp b/tests/test_satfunc.cpp
new file mode 100644
index 00000000..406ca81f
--- /dev/null
+++ b/tests/test_satfunc.cpp
@@ -0,0 +1,648 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+*/
+
+#include "config.h"
+
+/* --- Boost.Test boilerplate --- */
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE  // Suppress own messages when throw()ing
+
+#define BOOST_TEST_MODULE UnitsTest
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+/* --- our own headers --- */
+
+#include <opm/core/simulator/initStateEquil.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/cart_grid.h>
+#include <opm/core/grid/GridManager.hpp>
+
+#include <opm/core/props/BlackoilPropertiesBasic.hpp>
+#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
+#include <opm/core/pressure/msmfem/partition.h>
+
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/utility/Units.hpp>
+
+#include <array>
+#include <iostream>
+#include <limits>
+#include <memory>
+#include <numeric>
+#include <sstream>
+#include <string>
+#include <vector>
+
+BOOST_AUTO_TEST_SUITE ()
+
+BOOST_AUTO_TEST_CASE (GwsegStandard)
+{
+    // This is the basic (no eps and hysteris) version of 
+    // the Gwseg model.
+    
+    //std::cout << "==================================== GwsegStandard ====================================" << std::endl;
+
+    Opm::parameter::ParameterGroup param;
+    
+    Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0);
+    const UnstructuredGrid& grid = *(gm.c_grid());
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("satfuncStandard.DATA");
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
+    Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
+    
+    const int np = props.numPhases();
+    const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua];
+    const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid];
+    const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour];
+    
+    BOOST_REQUIRE(np == 3);
+    BOOST_REQUIRE(wpos == 0);
+    BOOST_REQUIRE(opos == 1);
+    BOOST_REQUIRE(gpos == 2);
+    
+    const int n=11;
+    double s[n*np];
+    int cells[n];
+    double kr[n*np];
+    double dkrds[n*np*np];
+    
+    for (int i=0; i<n; ++i) {
+      cells[i] = 0;
+      s[i*np+wpos] = i*0.1;
+      s[i*np+opos] = 1.0-s[i*np+wpos];
+      s[i*np+gpos] = 0.0;
+    }
+    
+    props.relperm(n, s, cells, kr, dkrds);
+    
+    double krw[11] = {0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7};
+    double kro[11] = {1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0};
+    double DkrwDsw[11] = {0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0};
+    double DkroDsw[11] = {-2.0, -2.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0};
+    double DkroDsg[11] = {-5.0, -5.0, -3.0, -2.0, -0.66666666666666741, -0.75, -0.8, 
+                          -0.83333333333333237, 0.14285714285714296, 0.0, 0.0};
+    
+    const double reltol = 1.0e-6;
+    for (int i=0; i<n; ++i) {
+      BOOST_CHECK_CLOSE(kr[i*np+wpos], krw[i], reltol);
+      BOOST_CHECK_CLOSE(kr[i*np+opos], kro[i], reltol);
+      BOOST_CHECK_CLOSE(dkrds[i*np*np+wpos], DkrwDsw[i], reltol);
+      BOOST_CHECK_CLOSE(dkrds[i*np*np+opos], DkroDsw[i], reltol);
+      BOOST_CHECK_CLOSE(dkrds[i*np*np+np*gpos+opos], DkroDsg[i], reltol);
+    }
+
+/*    
+    std::cout << std::setw(12) << "sw";
+    std::cout << std::setw(12) << "so";
+    std::cout << std::setw(12) << "sg";
+    std::cout << std::setw(12) << "krw";
+    std::cout << std::setw(12) << "kro";
+    std::cout << std::setw(12) << "krg";
+    std::cout << std::setw(12) << "DkrwDsw";
+    std::cout << std::setw(12) << "DkroDsw";
+    std::cout << std::setw(12) << "DkrgDsw";
+    std::cout << std::setw(12) << "DkrwDso";
+    std::cout << std::setw(12) << "DkroDso";
+    std::cout << std::setw(12) << "DkrgDso";
+    std::cout << std::setw(12) << "DkrwDsg";
+    std::cout << std::setw(12) << "DkroDsg";
+    std::cout << std::setw(12) << "DkrgDsg";
+    std::cout << std::endl;
+    for (int i=0; i<n; ++i) {
+      std::cout << std::setw(12) << s[i*np+wpos] << std::setw(12) << s[i*np+opos] << std::setw(12) << s[i*np+gpos]
+                << std::setw(12) << kr[i*np+wpos] << std::setw(12) << kr[i*np+opos] << std::setw(12) << kr[i*np+gpos];
+      for (int j=0; j<np*np; ++j) {
+        std::cout << std::setw(12) << dkrds[i*np*np+j];
+      }
+      std::cout << std::endl;
+    }
+*/
+}
+
+BOOST_AUTO_TEST_CASE (GwsegEPSBase)
+{
+    // This is the eps (but no hysteris) version of the Gwseg model.
+    // However, only default scaling parameters, i.e no scaling.
+    
+    //std::cout << "==================================== GwsegEPSBase ====================================" << std::endl;
+    
+    Opm::parameter::ParameterGroup param;
+    
+    Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0);
+    const UnstructuredGrid& grid = *(gm.c_grid());
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("satfuncEPSBase.DATA");
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
+    Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
+    
+    const int np = props.numPhases();
+    const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua];
+    const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid];
+    const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour];
+    
+    BOOST_REQUIRE(np == 3);
+    BOOST_REQUIRE(wpos == 0);
+    BOOST_REQUIRE(opos == 1);
+    BOOST_REQUIRE(gpos == 2);
+    
+    const int n=11;
+    double s[n*np];
+    int cells[n];
+    double kr[n*np];
+    double dkrds[n*np*np];
+    
+    for (int i=0; i<n; ++i) {
+      cells[i] = 0;
+      s[i*np+wpos] = i*0.1;
+      s[i*np+opos] = 1.0-s[i*np+wpos];
+      s[i*np+gpos] = 0.0;
+    }
+    
+    props.relperm(n, s, cells, kr, dkrds);
+ 
+    double krw[11] = {0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7};
+    double kro[11] = {1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0};
+    double DkrwDsw[11] = {0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0};
+    double DkroDsw[11] = {-2.0, -2.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0};
+    double DkroDsg[11] = {-5.0, -5.0, -3.0, -2.0,-0.66666666666666741, -0.75, -0.8, -0.83333333333333237, 0.14285714285714296, 0.0, 0.0};
+   
+    const double reltol = 1.0e-6;
+    for (int i=0; i<n; ++i) {
+      BOOST_CHECK_CLOSE(kr[i*np+wpos], krw[i], reltol);
+      BOOST_CHECK_CLOSE(kr[i*np+opos], kro[i], reltol);
+      BOOST_CHECK_CLOSE(dkrds[i*np*np+wpos], DkrwDsw[i], reltol);
+      BOOST_CHECK_CLOSE(dkrds[i*np*np+opos], DkroDsw[i], reltol);
+      BOOST_CHECK_CLOSE(dkrds[i*np*np+np*gpos+opos], DkroDsg[i], reltol);
+    }
+ 
+/*
+    std::cout << std::setw(12) << "sw";
+    std::cout << std::setw(12) << "so";
+    std::cout << std::setw(12) << "sg";
+    std::cout << std::setw(12) << "krw";
+    std::cout << std::setw(12) << "kro";
+    std::cout << std::setw(12) << "krg";
+    std::cout << std::setw(12) << "DkrwDsw";
+    std::cout << std::setw(12) << "DkroDsw";
+    std::cout << std::setw(12) << "DkrgDsw";
+    std::cout << std::setw(12) << "DkrwDso";
+    std::cout << std::setw(12) << "DkroDso";
+    std::cout << std::setw(12) << "DkrgDso";
+    std::cout << std::setw(12) << "DkrwDsg";
+    std::cout << std::setw(12) << "DkroDsg";
+    std::cout << std::setw(12) << "DkrgDsg";
+    std::cout << std::endl;
+    for (int i=0; i<n; ++i) {
+      std::cout << std::setw(12) << s[i*np+wpos] << std::setw(12) << s[i*np+opos] << std::setw(12) << s[i*np+gpos]
+                << std::setw(12) << kr[i*np+wpos] << std::setw(12) << kr[i*np+opos] << std::setw(12) << kr[i*np+gpos];
+      for (int j=0; j<np*np; ++j) {
+        std::cout << std::setw(12) << dkrds[i*np*np+j];
+      }
+      std::cout << std::endl;
+    }
+*/
+}
+
+
+BOOST_AUTO_TEST_CASE (GwsegEPS_A)    
+{
+    // This is the eps (but no hysteris) version of the Gwseg model.
+    // Scaling parameters from keyword SWL etc.
+    
+    //std::cout << "==================================== GwsegEPS_A ====================================" << std::endl;
+    
+    Opm::parameter::ParameterGroup param;
+    
+    Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0);
+    const UnstructuredGrid& grid = *(gm.c_grid());
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("satfuncEPS_A.DATA");
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
+    Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
+    
+    const int np = props.numPhases();
+    const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua];
+    const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid];
+    const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour];
+    
+    BOOST_REQUIRE(np == 3);
+    BOOST_REQUIRE(wpos == 0);
+    BOOST_REQUIRE(opos == 1);
+    BOOST_REQUIRE(gpos == 2);
+    
+    const int n=11;
+    double s[n*np];
+    int cells[n];
+    double kr[n*np];
+    double dkrds[n*np*np];
+    
+    const int ncell = 8;
+ 
+    double krw[ncell][n] = {{0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7},
+                            {0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7, 0.7, 0.7},
+                            {0, 0, 0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7},
+                            {0, 0, 0, 0, 0, 0.233333, 0.466667, 0.7, 0.7, 0.7, 0.7},
+                            {0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7},
+                            {0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7, 0.7, 0.7},
+                            {0, 0, 0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7},
+                            {0, 0, 0, 0, 0, 0.233333, 0.466667, 0.7, 0.7, 0.7, 0.7}};
+    double kro[ncell][n] = {{1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0},
+                            {1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
+                            {1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
+                            {1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
+                            {1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
+                            {1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
+                            {1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
+                            {1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0}};
+    double DkrwDsw[ncell][n] = {{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
+                                {0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
+                                {0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
+                                {0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0},
+                                {0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
+                                {0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
+                                {0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
+                                {0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0}};
+    double DkroDsw[ncell][n] = {{0.0, 0.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0},
+                                {0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
+                                {0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
+                                {0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
+                                {0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
+                                {0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
+                                {0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
+                                {0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0}};
+    double DkroDsg[ncell][n] = {{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
+                                {-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
+                                {-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
+                                {-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
+                                {-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
+                                {-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
+                                {-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
+                                {-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0}}; 
+    
+    for (int icell=0; icell<ncell; ++icell) {
+      for (int i=0; i<n; ++i) {
+        cells[i] = icell;
+        s[i*np+wpos] = i*0.1;
+        s[i*np+opos] = 1.0-s[i*np+wpos];
+        s[i*np+gpos] = 0.0;
+      }
+      
+      props.relperm(n, s, cells, kr, dkrds);
+      
+      const double reltol = 1.0e-3;
+      for (int i=0; i<n; ++i) {
+        BOOST_CHECK_CLOSE(kr[i*np+wpos], krw[icell][i], reltol);
+        BOOST_CHECK_CLOSE(kr[i*np+opos], kro[icell][i], reltol);
+        BOOST_CHECK_CLOSE(dkrds[i*np*np+wpos], DkrwDsw[icell][i], reltol);
+        BOOST_CHECK_CLOSE(dkrds[i*np*np+opos], DkroDsw[icell][i], reltol);
+        BOOST_CHECK_CLOSE(dkrds[i*np*np+np*gpos+opos], DkroDsg[icell][i], reltol);
+      }
+      
+/*      
+      std::cout << std::setw(12) << "sw: ";
+      for (int i=0; i<n; ++i) std::cout << s[i*np+wpos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "so: ";
+      for (int i=0; i<n; ++i) std::cout << s[i*np+opos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "sg: ";
+      for (int i=0; i<n; ++i) std::cout << s[i*np+gpos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "krw: ";
+      for (int i=0; i<n; ++i) std::cout << kr[i*np+wpos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "kro: ";
+      for (int i=0; i<n; ++i) std::cout << kr[i*np+opos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "krg: ";
+      for (int i=0; i<n; ++i) std::cout << kr[i*np+gpos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "DkrwDsw: ";
+      for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+wpos*np+wpos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "DkroDsw: ";
+      for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+wpos*np+opos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "DkrgDsw: ";
+      for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+wpos*np+gpos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "DkrwDso: ";
+      for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+opos*np+wpos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "DkroDso: ";
+      for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+opos*np+opos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "DkrgDso: ";
+      for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+opos*np+gpos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "DkrwDsg: ";
+      for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+gpos*np+wpos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "DkroDsg: ";
+      for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+gpos*np+opos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::setw(12) << "DkrgDsg: ";
+      for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+gpos*np+gpos] << ", ";
+      std::cout << std::endl;
+      std::cout << std::endl;
+*/
+    
+    }
+
+}
+
+BOOST_AUTO_TEST_CASE (GwsegEPS_B)    
+{
+    // This is the eps (but no hysteris) version of the Gwseg model.
+    // Scaling parameters from ENPTVD table.
+    
+    //std::cout << "==================================== GwsegEPS_B ====================================" << std::endl;
+ /*   
+    Opm::parameter::ParameterGroup param;
+    
+    Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0);
+    const UnstructuredGrid& grid = *(gm.c_grid());
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("satfuncEPS_B.DATA");
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
+    Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
+    
+    const int np = props.numPhases();
+    const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua];
+    const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid];
+    const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour];
+    
+    BOOST_REQUIRE(np == 3);
+    BOOST_REQUIRE(wpos == 0);
+    BOOST_REQUIRE(opos == 1);
+    BOOST_REQUIRE(gpos == 2);
+    
+    const int n=11;
+    double s[n*np];
+    int cells[n];
+    double kr[n*np];
+    double dkrds[n*np*np];
+    
+    const int ncell = 8;
+ 
+    double krw[ncell][n] = {{0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7},
+                            {0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7, 0.7, 0.7},
+                            {0, 0, 0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7},
+                            {0, 0, 0, 0, 0, 0.233333, 0.466667, 0.7, 0.7, 0.7, 0.7},
+                            {0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7},
+                            {0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7, 0.7, 0.7},
+                            {0, 0, 0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7},
+                            {0, 0, 0, 0, 0, 0.233333, 0.466667, 0.7, 0.7, 0.7, 0.7}};
+    double kro[ncell][n] = {{1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0},
+                            {1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
+                            {1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
+                            {1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
+                            {1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
+                            {1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
+                            {1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
+                            {1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0}};
+    double DkrwDsw[ncell][n] = {{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
+                                {0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
+                                {0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
+                                {0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0},
+                                {0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
+                                {0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
+                                {0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
+                                {0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0}};
+    double DkroDsw[ncell][n] = {{0.0, 0.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0},
+                                {0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
+                                {0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
+                                {0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
+                                {0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
+                                {0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
+                                {0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
+                                {0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0}};
+    double DkroDsg[ncell][n] = {{-2.32831e-10, -2.32831e-10, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
+                                {-2.32831e-10, -2.32831e-10, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
+                                {-2.32831e-10, -2.32831e-10, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
+                                {-2.32831e-10, -2.32831e-10, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
+                                {-2.32831e-10, -2.32831e-10, -2.32831e-10, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
+                                {-2.32831e-10, -2.32831e-10, -2.32831e-10, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
+                                {-2.32831e-10, -2.32831e-10, -2.32831e-10, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
+                                {-2.32831e-10, -2.32831e-10, -2.32831e-10, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0}}; 
+    
+    for (int icell=0; icell<ncell; ++icell) {
+      for (int i=0; i<n; ++i) {
+        cells[i] = icell;
+        s[i*np+wpos] = i*0.1;
+        s[i*np+opos] = 1.0-s[i*np+wpos];
+        s[i*np+gpos] = 0.0;
+      }
+      
+      props.relperm(n, s, cells, kr, dkrds);
+      
+      const double reltol = 1.0e-3;
+      for (int i=0; i<n; ++i) {
+        BOOST_CHECK_CLOSE(kr[i*np+wpos], krw[icell][i], reltol);
+        BOOST_CHECK_CLOSE(kr[i*np+opos], kro[icell][i], reltol);
+        BOOST_CHECK_CLOSE(dkrds[i*np*np+wpos], DkrwDsw[icell][i], reltol);
+        BOOST_CHECK_CLOSE(dkrds[i*np*np+opos], DkroDsw[icell][i], reltol);
+        BOOST_CHECK_CLOSE(dkrds[i*np*np+np*gpos+opos], DkroDsg[icell][i], reltol);
+      }
+  
+    }
+    */
+
+}
+
+BOOST_AUTO_TEST_CASE (GwsegEPS_C)    
+{
+    // This is the eps (but no hysteris) version of the Gwseg model.
+    // Scaling parameters given the "norne-way", i.e EQUALS, COPY, ADD, MULTIPLY.
+    
+    //std::cout << "==================================== GwsegEPS_C ====================================" << std::endl;
+    
+    Opm::parameter::ParameterGroup param;
+    
+    Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0);
+    const UnstructuredGrid& grid = *(gm.c_grid());
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("satfuncEPS_C.DATA");
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
+    Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
+    
+    const int np = props.numPhases();
+    const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua];
+    const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid];
+    const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour];
+    
+    BOOST_REQUIRE(np == 3);
+    BOOST_REQUIRE(wpos == 0);
+    BOOST_REQUIRE(opos == 1);
+    BOOST_REQUIRE(gpos == 2);
+    
+    const int n=11;
+    double s[n*np];
+    int cells[n];
+    double kr[n*np];
+    double dkrds[n*np*np];
+    
+    const int ncell = 8;
+ 
+    double krw[ncell][n] = {{0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7},
+                            {0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7, 0.7, 0.7},
+                            {0, 0, 0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7},
+                            {0, 0, 0, 0, 0, 0.233333, 0.466667, 0.7, 0.7, 0.7, 0.7},
+                            {0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7},
+                            {0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7, 0.7, 0.7},
+                            {0, 0, 0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7},
+                            {0, 0, 0, 0, 0, 0.233333, 0.466667, 0.7, 0.7, 0.7, 0.7}};
+    double kro[ncell][n] = {{1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0},
+                            {1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
+                            {1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
+                            {1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
+                            {1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
+                            {1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
+                            {1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
+                            {1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0}};
+    double DkrwDsw[ncell][n] = {{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
+                                {0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
+                                {0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
+                                {0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0},
+                                {0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
+                                {0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
+                                {0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
+                                {0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0}};
+    double DkroDsw[ncell][n] = {{0.0, 0.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0},
+                                {0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
+                                {0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
+                                {0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
+                                {0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
+                                {0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
+                                {0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
+                                {0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0}};
+    double DkroDsg[ncell][n] = {{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
+                                {-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
+                                {-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
+                                {-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
+                                {-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
+                                {-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
+                                {-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
+                                {-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0}}; 
+    
+    for (int icell=0; icell<ncell; ++icell) {
+      for (int i=0; i<n; ++i) {
+        cells[i] = icell;
+        s[i*np+wpos] = i*0.1;
+        s[i*np+opos] = 1.0-s[i*np+wpos];
+        s[i*np+gpos] = 0.0;
+      }
+      
+      props.relperm(n, s, cells, kr, dkrds);
+      
+      const double reltol = 1.0e-3;
+      for (int i=0; i<n; ++i) {
+        BOOST_CHECK_CLOSE(kr[i*np+wpos], krw[icell][i], reltol);
+        BOOST_CHECK_CLOSE(kr[i*np+opos], kro[icell][i], reltol);
+        BOOST_CHECK_CLOSE(dkrds[i*np*np+wpos], DkrwDsw[icell][i], reltol);
+        BOOST_CHECK_CLOSE(dkrds[i*np*np+opos], DkroDsw[icell][i], reltol);
+        BOOST_CHECK_CLOSE(dkrds[i*np*np+np*gpos+opos], DkroDsg[icell][i], reltol);
+      }
+    
+    }
+
+}
+
+BOOST_AUTO_TEST_CASE (GwsegEPS_D)
+{
+    // This is the eps and hysteris version of the Gwseg model.
+    // However, only default scaling parameters, i.e no scaling.
+    
+    //std::cout << "==================================== GwsegEPS_D ====================================" << std::endl;
+
+    Opm::parameter::ParameterGroup param;
+    
+    Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0);
+    const UnstructuredGrid& grid = *(gm.c_grid());
+    Opm::ParserPtr parser(new Opm::Parser() );
+    Opm::DeckConstPtr deck = parser->parseFile("satfuncEPS_D.DATA");
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
+    Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
+    
+    const int np = props.numPhases();
+    const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua];
+    const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid];
+    const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour];
+    
+    BOOST_REQUIRE(np == 3);
+    BOOST_REQUIRE(wpos == 0);
+    BOOST_REQUIRE(opos == 1);
+    BOOST_REQUIRE(gpos == 2);
+    
+    const int n=11;
+    double s[n*np];
+    int cells[n];
+    double kr[n*np];
+    double dkrds[n*np*np];
+    
+    for (int i=0; i<n; ++i) {
+      cells[i] = 0;
+      s[i*np+wpos] = i*0.1;
+      s[i*np+opos] = 1.0-s[i*np+wpos];
+      s[i*np+gpos] = 0.0;
+    }
+    
+    props.relperm(n, s, cells, kr, dkrds);
+    
+    double krw[11] = {0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7};
+    double kro[11] = {1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0};
+    double DkrwDsw[11] = {0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0};
+    double DkroDsw[11] = {-2.0, -2.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0};
+    double DkroDsg[11] = {-5.0, -5.0, -3.0, -2.0, -0.66666666666666741, -0.75, -0.8, 
+                          -0.83333333333333237, 0.14285714285714296, 0.0, 0.0};
+    
+    const double reltol = 1.0e-6;
+    for (int i=0; i<n; ++i) {
+      BOOST_CHECK_CLOSE(kr[i*np+wpos], krw[i], reltol);
+      BOOST_CHECK_CLOSE(kr[i*np+opos], kro[i], reltol);
+      BOOST_CHECK_CLOSE(dkrds[i*np*np+wpos], DkrwDsw[i], reltol);
+      BOOST_CHECK_CLOSE(dkrds[i*np*np+opos], DkroDsw[i], reltol);
+      BOOST_CHECK_CLOSE(dkrds[i*np*np+np*gpos+opos], DkroDsg[i], reltol);
+    }
+
+/*    
+    std::cout << std::setw(12) << "sw";
+    std::cout << std::setw(12) << "so";
+    std::cout << std::setw(12) << "sg";
+    std::cout << std::setw(12) << "krw";
+    std::cout << std::setw(12) << "kro";
+    std::cout << std::setw(12) << "krg";
+    std::cout << std::setw(12) << "DkrwDsw";
+    std::cout << std::setw(12) << "DkroDsw";
+    std::cout << std::setw(12) << "DkrgDsw";
+    std::cout << std::setw(12) << "DkrwDso";
+    std::cout << std::setw(12) << "DkroDso";
+    std::cout << std::setw(12) << "DkrgDso";
+    std::cout << std::setw(12) << "DkrwDsg";
+    std::cout << std::setw(12) << "DkroDsg";
+    std::cout << std::setw(12) << "DkrgDsg";
+    std::cout << std::endl;
+    for (int i=0; i<n; ++i) {
+      std::cout << std::setw(12) << s[i*np+wpos] << std::setw(12) << s[i*np+opos] << std::setw(12) << s[i*np+gpos]
+                << std::setw(12) << kr[i*np+wpos] << std::setw(12) << kr[i*np+opos] << std::setw(12) << kr[i*np+gpos];
+      for (int j=0; j<np*np; ++j) {
+        std::cout << std::setw(12) << dkrds[i*np*np+j];
+      }
+      std::cout << std::endl;
+    }
+*/
+}
+
+BOOST_AUTO_TEST_SUITE_END()