diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index f86e19cc..3e5db299 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -84,11 +84,11 @@ list (APPEND MAIN_SOURCE_FILES
 	opm/core/props/pvt/BlackoilPvtProperties.cpp
 	opm/core/props/pvt/PvtPropertiesBasic.cpp
 	opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
-	opm/core/props/pvt/SinglePvtDead.cpp
-	opm/core/props/pvt/SinglePvtDeadSpline.cpp
-	opm/core/props/pvt/SinglePvtInterface.cpp
-	opm/core/props/pvt/SinglePvtLiveGas.cpp
-	opm/core/props/pvt/SinglePvtLiveOil.cpp
+	opm/core/props/pvt/PvtDead.cpp
+	opm/core/props/pvt/PvtDeadSpline.cpp
+	opm/core/props/pvt/PvtInterface.cpp
+	opm/core/props/pvt/PvtLiveGas.cpp
+	opm/core/props/pvt/PvtLiveOil.cpp
 	opm/core/props/rock/RockBasic.cpp
 	opm/core/props/rock/RockCompressibility.cpp
 	opm/core/props/rock/RockFromDeck.cpp
@@ -311,12 +311,12 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/core/props/pvt/BlackoilPvtProperties.hpp
 	opm/core/props/pvt/PvtPropertiesBasic.hpp
 	opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp
-	opm/core/props/pvt/SinglePvtConstCompr.hpp
-	opm/core/props/pvt/SinglePvtDead.hpp
-	opm/core/props/pvt/SinglePvtDeadSpline.hpp
-	opm/core/props/pvt/SinglePvtInterface.hpp
-	opm/core/props/pvt/SinglePvtLiveGas.hpp
-	opm/core/props/pvt/SinglePvtLiveOil.hpp
+	opm/core/props/pvt/PvtConstCompr.hpp
+	opm/core/props/pvt/PvtDead.hpp
+	opm/core/props/pvt/PvtDeadSpline.hpp
+	opm/core/props/pvt/PvtInterface.hpp
+	opm/core/props/pvt/PvtLiveGas.hpp
+	opm/core/props/pvt/PvtLiveOil.hpp
 	opm/core/props/rock/RockBasic.hpp
 	opm/core/props/rock/RockCompressibility.hpp
 	opm/core/props/rock/RockFromDeck.hpp
diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp
index 606e04ba..24ac6438 100644
--- a/opm/core/pressure/CompressibleTpfa.cpp
+++ b/opm/core/pressure/CompressibleTpfa.cpp
@@ -238,7 +238,7 @@ namespace Opm
                 const int cell = wells_->well_cells[j];
                 const double cell_depth = grid_.cell_centroids[dim * cell + dim - 1];
                 props_.matrix(1, &state.pressure()[cell], &state.surfacevol()[np*cell], &cell, &A[0], 0);
-                props_.density(1, &A[0], &rho[0]);
+                props_.density(1, &A[0], &cell, &rho[0]);
                 for (int phase = 0; phase < np; ++phase) {
                     const double s_phase = state.saturation()[np*cell + phase];
                     wellperf_wdp_[j] += s_phase*rho[phase]*grav*(cell_depth - ref_depth);
@@ -380,7 +380,7 @@ namespace Opm
                     // Gravity contribution, gravcontrib = rho*(face_z - cell_z) [per phase].
                     if (grav != 0.0) {
                         const double depth_diff = face_depth - grid_.cell_centroids[c[j]*dim + dim - 1];
-                        props_.density(1, &cell_A_[np*np*c[j]], &gravcontrib[j][0]);
+                        props_.density(1, &cell_A_[np*np*c[j]], &c[j], &gravcontrib[j][0]);
                         for (int p = 0; p < np; ++p) {
                             gravcontrib[j][p] *= depth_diff*grav;
                         }
diff --git a/opm/core/props/BlackoilPropertiesBasic.cpp b/opm/core/props/BlackoilPropertiesBasic.cpp
index 8700af11..ab9b238b 100644
--- a/opm/core/props/BlackoilPropertiesBasic.cpp
+++ b/opm/core/props/BlackoilPropertiesBasic.cpp
@@ -157,9 +157,11 @@ namespace Opm
     ///                    matrix A = RB^{-1} which relates z to u by z = Au. The matrices
     ///                    are assumed to be in Fortran order, and are typically the result
     ///                    of a call to the method matrix().
+    /// \param[in]  cells  The index of the grid cell of each data point.
     /// \param[out] rho    Array of nP density values, array must be valid before calling.
     void BlackoilPropertiesBasic::density(const int n,
                                           const double* A,
+                                          const int* /*cells*/,
                                           double* rho) const
     {
         const int np = numPhases();
@@ -177,7 +179,7 @@ namespace Opm
 
     /// Densities of stock components at surface conditions.
     /// \return Array of P density values.
-    const double* BlackoilPropertiesBasic::surfaceDensity() const
+    const double* BlackoilPropertiesBasic::surfaceDensity(int /*cellIdx*/) const
     {
         return pvt_.surfaceDensities();
     }
diff --git a/opm/core/props/BlackoilPropertiesBasic.hpp b/opm/core/props/BlackoilPropertiesBasic.hpp
index ad84b1a3..b85991ed 100644
--- a/opm/core/props/BlackoilPropertiesBasic.hpp
+++ b/opm/core/props/BlackoilPropertiesBasic.hpp
@@ -59,6 +59,11 @@ namespace Opm
         /// \return   N, the number of cells.
         virtual int numCells() const;
 
+        /// Return an array containing the PVT table index for each
+        /// grid cell
+        virtual const int* cellPvtRegionIndex() const
+        { return NULL; }
+
         /// \return   Array of N porosity values.
         virtual const double* porosity() const;
 
@@ -114,14 +119,17 @@ namespace Opm
         ///                    matrix A = RB^{-1} which relates z to u by z = Au. The matrices
         ///                    are assumed to be in Fortran order, and are typically the result
         ///                    of a call to the method matrix().
+        /// \param[in]  cells  The index of the grid cell of each data point.
         /// \param[out] rho    Array of nP density values, array must be valid before calling.
         virtual void density(const int n,
                              const double* A,
+                             const int* cells,
                              double* rho) const;
 
         /// Densities of stock components at surface conditions.
+        /// \param[in]  cellIdx  The index of the cell for which the surface density ought to be calculated
         /// \return Array of P density values.
-        virtual const double* surfaceDensity() const;
+        virtual const double* surfaceDensity(int cellIdx = 0) const;
 
         /// \param[in]  n      Number of data points.
         /// \param[in]  s      Array of nP saturation values.
diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp
index 5a781031..2a2a020a 100644
--- a/opm/core/props/BlackoilPropertiesFromDeck.cpp
+++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp
@@ -96,14 +96,20 @@ namespace Opm
     void BlackoilPropertiesFromDeck::viscosity(const int n,
                                                const double* p,
                                                const double* z,
-                                               const int* /*cells*/,
+                                               const int* cells,
                                                double* mu,
                                                double* dmudp) const
     {
         if (dmudp) {
             OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::viscosity()  --  derivatives of viscosity not yet implemented.");
         } else {
-            pvt_.mu(n, p, z, mu);
+            const int *cellPvtTableIdx = cellPvtRegionIndex();
+            assert(cellPvtTableIdx != 0);
+            std::vector<int> pvtTableIdx(n);
+            for (int i = 0; i < n; ++ i)
+                pvtTableIdx[i] = cellPvtTableIdx[cells[i]];
+
+            pvt_.mu(n, &pvtTableIdx[0], p, z, mu);
         }
     }
 
@@ -120,21 +126,27 @@ namespace Opm
     void BlackoilPropertiesFromDeck::matrix(const int n,
                                             const double* p,
                                             const double* z,
-                                            const int* /*cells*/,
+                                            const int* cells,
                                             double* A,
                                             double* dAdp) const
     {
         const int np = numPhases();
+
+        const int *cellPvtTableIdx = cellPvtRegionIndex();
+        std::vector<int> pvtTableIdx(n);
+        for (int i = 0; i < n; ++ i)
+            pvtTableIdx[i] = cellPvtTableIdx[cells[i]];
+
         B_.resize(n*np);
         R_.resize(n*np);
         if (dAdp) {
             dB_.resize(n*np);
             dR_.resize(n*np);
-            pvt_.dBdp(n, p, z, &B_[0], &dB_[0]);
-            pvt_.dRdp(n, p, z, &R_[0], &dR_[0]);
+            pvt_.dBdp(n, &pvtTableIdx[0], p, z, &B_[0], &dB_[0]);
+            pvt_.dRdp(n, &pvtTableIdx[0], p, z, &R_[0], &dR_[0]);
         } else {
-            pvt_.B(n, p, z, &B_[0]);
-            pvt_.R(n, p, z, &R_[0]);
+            pvt_.B(n, &pvtTableIdx[0], p, z, &B_[0]);
+            pvt_.R(n, &pvtTableIdx[0], p, z, &R_[0]);
         }
         const int* phase_pos = pvt_.phasePosition();
         bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] &&
@@ -207,15 +219,19 @@ namespace Opm
     ///                    matrix A = RB^{-1} which relates z to u by z = Au. The matrices
     ///                    are assumed to be in Fortran order, and are typically the result
     ///                    of a call to the method matrix().
+    /// \param[in]  cells  The index of the grid cell of each data point.
     /// \param[out] rho    Array of nP density values, array must be valid before calling.
     void BlackoilPropertiesFromDeck::density(const int n,
                                              const double* A,
+                                             const int* cells,
                                              double* rho) const
     {
         const int np = numPhases();
-        const double* sdens = pvt_.surfaceDensities();
 // #pragma omp parallel for
         for (int i = 0; i < n; ++i) {
+            int cellIdx = cells?cells[i]:i;
+            int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx);
+            const double* sdens = pvt_.surfaceDensities(pvtRegionIdx);
             for (int phase = 0; phase < np; ++phase) {
                 rho[np*i + phase] = 0.0;
                 for (int comp = 0; comp < np; ++comp) {
@@ -227,9 +243,10 @@ namespace Opm
 
     /// Densities of stock components at surface conditions.
     /// \return Array of P density values.
-    const double* BlackoilPropertiesFromDeck::surfaceDensity() const
+    const double* BlackoilPropertiesFromDeck::surfaceDensity(int cellIdx) const
     {
-        return pvt_.surfaceDensities();
+        int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx);
+        return pvt_.surfaceDensities(pvtRegionIdx);
     }
 
     /// \param[in]  n      Number of data points.
diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp
index e2919424..448c5cb5 100644
--- a/opm/core/props/BlackoilPropertiesFromDeck.hpp
+++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp
@@ -97,6 +97,11 @@ namespace Opm
         /// \return   N, the number of cells.
         virtual int numCells() const;
 
+        /// Return an array containing the PVT table index for each
+        /// grid cell
+        virtual const int* cellPvtRegionIndex() const
+        { return &cellPvtRegionIdx_[0]; }
+
         /// \return   Array of N porosity values.
         virtual const double* porosity() const;
 
@@ -152,14 +157,16 @@ namespace Opm
         ///                    matrix A = RB^{-1} which relates z to u by z = Au. The matrices
         ///                    are assumed to be in Fortran order, and are typically the result
         ///                    of a call to the method matrix().
+        /// \param[in]  cells  The index of the grid cell of each data point.
         /// \param[out] rho    Array of nP density values, array must be valid before calling.
         virtual void density(const int n,
                              const double* A,
+                             const int* cells,
                              double* rho) const;
 
         /// Densities of stock components at surface conditions.
         /// \return Array of P density values.
-        virtual const double* surfaceDensity() const;
+        virtual const double* surfaceDensity(int cellIdx = 0) const;
 
         /// \param[in]  n      Number of data points.
         /// \param[in]  s      Array of nP saturation values.
@@ -206,6 +213,13 @@ namespace Opm
                               double* smax) const;
 
     private:
+        int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
+        {
+            if (!pvtTableIdx)
+                return 0;
+            return pvtTableIdx[cellIdx];
+        }
+
         template<class CentroidIterator>
         void init(Opm::DeckConstPtr deck,
                   int number_of_cells,
@@ -224,6 +238,7 @@ namespace Opm
                   const parameter::ParameterGroup& param,
                   bool init_rock);
         RockFromDeck rock_;
+        std::vector<int> cellPvtRegionIdx_;
         BlackoilPvtProperties pvt_;
         std::unique_ptr<SaturationPropsInterface> satprops_;
         mutable std::vector<double> B_;
diff --git a/opm/core/props/BlackoilPropertiesFromDeck_impl.hpp b/opm/core/props/BlackoilPropertiesFromDeck_impl.hpp
index 9a9b8bad..586021c4 100644
--- a/opm/core/props/BlackoilPropertiesFromDeck_impl.hpp
+++ b/opm/core/props/BlackoilPropertiesFromDeck_impl.hpp
@@ -23,8 +23,14 @@ namespace Opm
                                                            const parameter::ParameterGroup& param,
                                                            bool init_rock)
     {
-        init(deck, number_of_cells, global_cell, cart_dims, begin_cell_centroids, dimension,
-             param, init_rock);
+        init(deck,
+             number_of_cells,
+             global_cell,
+             cart_dims,
+             begin_cell_centroids,
+             dimension,
+             param,
+             init_rock);
     }
 
     template<class CentroidIterator>
@@ -36,6 +42,10 @@ namespace Opm
                                                  int dimension,
                                                  bool init_rock)
     {
+        // retrieve the cell specific PVT table index from the deck
+        // and using the grid...
+        extractPvtTableIndex(cellPvtRegionIdx_, deck, number_of_cells, global_cell);
+
         if (init_rock){
            rock_.init(deck, number_of_cells, global_cell, cart_dims);
         }
@@ -62,6 +72,9 @@ namespace Opm
                                                  const parameter::ParameterGroup& param,
                                                  bool init_rock)
     {
+        // retrieve the cell specific PVT table index from the deck
+        // and using the grid...
+        extractPvtTableIndex(cellPvtRegionIdx_, deck, number_of_cells, global_cell);
 
         if(init_rock){
             rock_.init(deck, number_of_cells, global_cell, cart_dims);
diff --git a/opm/core/props/BlackoilPropertiesInterface.hpp b/opm/core/props/BlackoilPropertiesInterface.hpp
index 80e3ef7b..08edc8ae 100644
--- a/opm/core/props/BlackoilPropertiesInterface.hpp
+++ b/opm/core/props/BlackoilPropertiesInterface.hpp
@@ -47,6 +47,10 @@ namespace Opm
         /// \return   N, the number of cells.
         virtual int numCells() const = 0;
 
+        /// Return an array containing the PVT table index for each
+        /// grid cell
+        virtual const int* cellPvtRegionIndex() const = 0;
+
         /// \return   Array of N porosity values.
         virtual const double* porosity() const = 0;
 
@@ -102,14 +106,16 @@ namespace Opm
         ///                    matrix A = RB^{-1} which relates z to u by z = Au. The matrices
         ///                    are assumed to be in Fortran order, and are typically the result
         ///                    of a call to the method matrix().
+        /// \param[in]  cells  The index of the grid cell of each data point.
         /// \param[out] rho    Array of nP density values, array must be valid before calling.
         virtual void density(const int n,
                              const double* A,
+                             const int* cells,
                              double* rho) const = 0;
 
         /// Densities of stock components at surface conditions.
         /// \return Array of P density values.
-        virtual const double* surfaceDensity() const = 0;
+        virtual const double* surfaceDensity(int regionIdx = 0) const = 0;
 
         /// \param[in]  n      Number of data points.
         /// \param[in]  s      Array of nP saturation values.
diff --git a/opm/core/props/pvt/BlackoilPvtProperties.cpp b/opm/core/props/pvt/BlackoilPvtProperties.cpp
index 04e56ad1..d1026c00 100644
--- a/opm/core/props/pvt/BlackoilPvtProperties.cpp
+++ b/opm/core/props/pvt/BlackoilPvtProperties.cpp
@@ -21,11 +21,11 @@
 
 #include "config.h"
 #include <opm/core/props/pvt/BlackoilPvtProperties.hpp>
-#include <opm/core/props/pvt/SinglePvtDead.hpp>
-#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
-#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
-#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
-#include <opm/core/props/pvt/SinglePvtConstCompr.hpp>
+#include <opm/core/props/pvt/PvtDead.hpp>
+#include <opm/core/props/pvt/PvtDeadSpline.hpp>
+#include <opm/core/props/pvt/PvtLiveOil.hpp>
+#include <opm/core/props/pvt/PvtLiveGas.hpp>
+#include <opm/core/props/pvt/PvtConstCompr.hpp>
 #include <opm/core/props/phaseUsageFromDeck.hpp>
 #include <opm/core/utility/Units.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
@@ -42,93 +42,97 @@ namespace Opm
     {
     }
 
-    void BlackoilPvtProperties::init(Opm::DeckConstPtr deck, int samples)
+    void BlackoilPvtProperties::init(Opm::DeckConstPtr deck,
+                                     int numSamples)
     {
-        // If we need multiple regions, this class and the SinglePvt* classes must change.
-        region_number_ = 0;
-
         phase_usage_ = phaseUsageFromDeck(deck);
 
         // Surface densities. Accounting for different orders in eclipse and our code.
-        if (deck->hasKeyword("DENSITY")) {
-            Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY");
+        Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY");
+        int numRegions = densityKeyword->size();
+
+        densities_.resize(numRegions);
+        for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
             if (phase_usage_.phase_used[Liquid]) {
-                densities_[phase_usage_.phase_pos[Liquid]]
-                    = densityKeyword->getRecord(region_number_)->getItem("OIL")->getSIDouble(0);
+                densities_[regionIdx][phase_usage_.phase_pos[Liquid]]
+                    = densityKeyword->getRecord(regionIdx)->getItem("OIL")->getSIDouble(0);
             }
             if (phase_usage_.phase_used[Aqua]) {
-                densities_[phase_usage_.phase_pos[Aqua]]
-                    = densityKeyword->getRecord(region_number_)->getItem("WATER")->getSIDouble(0);
+                densities_[regionIdx][phase_usage_.phase_pos[Aqua]]
+                    = densityKeyword->getRecord(regionIdx)->getItem("WATER")->getSIDouble(0);
             }
             if (phase_usage_.phase_used[Vapour]) {
-                densities_[phase_usage_.phase_pos[Vapour]]
-                    = densityKeyword->getRecord(region_number_)->getItem("GAS")->getSIDouble(0);
+                densities_[regionIdx][phase_usage_.phase_pos[Vapour]]
+                    = densityKeyword->getRecord(regionIdx)->getItem("GAS")->getSIDouble(0);
             }
-        } else {
-            OPM_THROW(std::runtime_error, "Input is missing DENSITY\n");
         }
 
-        // Set the properties.
+        // Resize the property objects container
         props_.resize(phase_usage_.num_phases);
+
         // Water PVT
         if (phase_usage_.phase_used[Aqua]) {
-            if (deck->hasKeyword("PVTW")) {
-                Opm::PvtwTable pvtwTable(deck->getKeyword("PVTW"));
+            // if water is used, we require the presence of the "PVTW"
+            // keyword for now...
+            std::shared_ptr<PvtConstCompr> pvtw(new PvtConstCompr);
+            pvtw->initFromWater(deck->getKeyword("PVTW"));
 
-                props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable));
-            } else {
-                // Eclipse 100 default.
-                props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
-            }
+            props_[phase_usage_.phase_pos[Aqua]] = pvtw;
         }
         // Oil PVT
         if (phase_usage_.phase_used[Liquid]) {
+            // for oil, we support the "PVDO", "PVTO" and "PVCDO"
+            // keywords...
             if (deck->hasKeyword("PVDO")) {
-                Opm::PvdoTable pvdoTable(deck->getKeyword("PVDO"), region_number_);
-                if (samples > 0) {
-                    props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples));
+                Opm::DeckKeywordConstPtr pvdoKeyword = deck->getKeyword("PVDO");
+                if (numSamples > 0) {
+                    auto splinePvt = std::shared_ptr<PvtDeadSpline>(new PvtDeadSpline);
+                    splinePvt->initFromOil(pvdoKeyword, numSamples);
+                    props_[phase_usage_.phase_pos[Liquid]] = splinePvt;
                 } else {
-                    props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(pvdoTable));
+                    auto deadPvt = std::shared_ptr<PvtDead>(new PvtDead);
+                    deadPvt->initFromOil(pvdoKeyword);
+                    props_[phase_usage_.phase_pos[Liquid]] = deadPvt;
                 }
             } else if (deck->hasKeyword("PVTO")) {
-                Opm::PvtoTable pvtoTable(deck->getKeyword("PVTO"), /*tableIdx=*/0);
-
-                props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable));
+                props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO")));
             } else if (deck->hasKeyword("PVCDO")) {
-                Opm::PvcdoTable pvcdoTable(deck->getKeyword("PVCDO"));
+                std::shared_ptr<PvtConstCompr> pvcdo(new PvtConstCompr);
+                pvcdo->initFromOil(deck->getKeyword("PVCDO"));
 
-                props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable));
+                props_[phase_usage_.phase_pos[Liquid]] = pvcdo;
             } else {
-                OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n");
+                OPM_THROW(std::runtime_error, "Input is missing PVDO, PVCDO or PVTO\n");
             }
         }
         // Gas PVT
         if (phase_usage_.phase_used[Vapour]) {
+            // gas can be specified using the "PVDG" or "PVTG" keywords...
             if (deck->hasKeyword("PVDG")) {
-                Opm::PvdgTable pvdgTable(deck->getKeyword("PVDG"), region_number_);
-                if (samples > 0) {
-                    props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples));
+                Opm::DeckKeywordConstPtr pvdgKeyword = deck->getKeyword("PVDG");
+
+                if (numSamples > 0) {
+                    std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
+                    splinePvt->initFromGas(pvdgKeyword, numSamples);
+
+                    props_[phase_usage_.phase_pos[Vapour]] = splinePvt;
                 } else {
-                    props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(pvdgTable));
+                    std::shared_ptr<PvtDead> deadPvt(new PvtDead);
+                    deadPvt->initFromGas(pvdgKeyword);
+
+                    props_[phase_usage_.phase_pos[Vapour]] = deadPvt;
                 }
             } else if (deck->hasKeyword("PVTG")) {
-                Opm::PvtgTable pvtgTable(deck->getKeyword("PVTG"), /*tableIdx=*/0);
-
-                props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable));
+                props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG")));
             } else {
                 OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
             }
         }
-
-        // Must inform pvt property objects of phase structure.
-        for (int i = 0; i < phase_usage_.num_phases; ++i) {
-            props_[i]->setPhaseConfiguration(phase_usage_.num_phases, phase_usage_.phase_pos);
-        }
     }
 
-    const double* BlackoilPvtProperties::surfaceDensities() const
+    const double* BlackoilPvtProperties::surfaceDensities(int regionIdx) const
     {
-        return densities_;
+        return &densities_[regionIdx][0];
     }
 
 
@@ -154,13 +158,14 @@ namespace Opm
 
 
     void BlackoilPvtProperties::mu(const int n,
+                                   const int* pvtTableIdx,
                                    const double* p,
                                    const double* z,
                                    double* output_mu) const
     {
         data1_.resize(n);
         for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
-            props_[phase]->mu(n, p, z, &data1_[0]);
+            props_[phase]->mu(n, pvtTableIdx, p, z, &data1_[0]);
 // #pragma omp parallel for
             for (int i = 0; i < n; ++i) {
                 output_mu[phase_usage_.num_phases*i + phase] = data1_[i];
@@ -169,13 +174,14 @@ namespace Opm
     }
 
     void BlackoilPvtProperties::B(const int n,
+                                  const int* pvtTableIdx,
                                   const double* p,
                                   const double* z,
                                   double* output_B) const
     {
         data1_.resize(n);
         for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
-            props_[phase]->B(n, p, z, &data1_[0]);
+            props_[phase]->B(n, pvtTableIdx, p, z, &data1_[0]);
 // #pragma omp parallel for
             for (int i = 0; i < n; ++i) {
                 output_B[phase_usage_.num_phases*i + phase] = data1_[i];
@@ -184,6 +190,7 @@ namespace Opm
     }
 
     void BlackoilPvtProperties::dBdp(const int n,
+                                     const int* pvtTableIdx,
                                      const double* p,
                                      const double* z,
                                      double* output_B,
@@ -192,7 +199,7 @@ namespace Opm
         data1_.resize(n);
         data2_.resize(n);
         for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
-            props_[phase]->dBdp(n, p, z, &data1_[0], &data2_[0]);
+            props_[phase]->dBdp(n, pvtTableIdx, p, z, &data1_[0], &data2_[0]);
 // #pragma omp parallel for
             for (int i = 0; i < n; ++i) {
                 output_B[phase_usage_.num_phases*i + phase] = data1_[i];
@@ -203,13 +210,14 @@ namespace Opm
 
 
     void BlackoilPvtProperties::R(const int n,
+                                  const int* pvtTableIdx,
                                   const double* p,
                                   const double* z,
                                   double* output_R) const
     {
         data1_.resize(n);
         for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
-            props_[phase]->R(n, p, z, &data1_[0]);
+            props_[phase]->R(n, pvtTableIdx, p, z, &data1_[0]);
 // #pragma omp parallel for
             for (int i = 0; i < n; ++i) {
                 output_R[phase_usage_.num_phases*i + phase] = data1_[i];
@@ -218,6 +226,7 @@ namespace Opm
     }
 
     void BlackoilPvtProperties::dRdp(const int n,
+                                     const int* pvtTableIdx,
                                      const double* p,
                                      const double* z,
                                      double* output_R,
@@ -226,7 +235,7 @@ namespace Opm
         data1_.resize(n);
         data2_.resize(n);
         for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
-            props_[phase]->dRdp(n, p, z, &data1_[0], &data2_[0]);
+            props_[phase]->dRdp(n, pvtTableIdx, p, z, &data1_[0], &data2_[0]);
 // #pragma omp parallel for
             for (int i = 0; i < n; ++i) {
                 output_R[phase_usage_.num_phases*i + phase] = data1_[i];
@@ -234,5 +243,4 @@ namespace Opm
             }
         }
     }
-
 } // namespace Opm
diff --git a/opm/core/props/pvt/BlackoilPvtProperties.hpp b/opm/core/props/pvt/BlackoilPvtProperties.hpp
index b8bff837..a4a44adc 100644
--- a/opm/core/props/pvt/BlackoilPvtProperties.hpp
+++ b/opm/core/props/pvt/BlackoilPvtProperties.hpp
@@ -20,26 +20,30 @@
 #ifndef OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED
 #define OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED
 
-#include <opm/core/props/pvt/SinglePvtInterface.hpp>
+#include <opm/core/props/pvt/PvtInterface.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
 
 #include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Utility/PvtoTable.hpp>
 
 #include <string>
 #include <memory>
+#include <array>
 
 namespace Opm
 {
 
     /// Class collecting the pvt properties for all active phases.
-    /// For all the methods, the following apply: p and z
-    /// are expected to be of size n and n*num_phases, respectively.
-    /// Output arrays shall be of size n*num_phases, and must be valid
-    /// before calling the method.
+    /// For all the methods, the following apply:
+    /// - p and z are expected to be of size n and n*num_phases, respectively.
+    /// - pvtTableIdx specifies the PVT table to be used for each data
+    ///               point and is thus expected to be an array of size n
+    /// - Output arrays shall be of size n*num_phases, and must be valid
+    ///   before calling the method.
     /// NOTE: The difference between this interface and the one defined
-    /// by SinglePvtInterface is that this collects all phases' properties,
+    /// by PvtInterface is that this collects all phases' properties,
     /// and therefore the output arrays are of size n*num_phases as opposed
-    /// to size n in SinglePvtInterface.
+    /// to size n in PvtInterface.
     class BlackoilPvtProperties : public BlackoilPhases
     {
     public:
@@ -47,8 +51,10 @@ namespace Opm
         BlackoilPvtProperties();
 
         /// Initialize from deck.
+        ///
         /// \param deck     An input deck from the opm-parser module.
-        void init(Opm::DeckConstPtr deck, int samples);
+        void init(Opm::DeckConstPtr deck,
+                  int samples);
 
         /// \return   Object describing the active phases.
         PhaseUsage phaseUsage() const;
@@ -68,22 +74,25 @@ namespace Opm
 
         /// Densities of stock components at surface conditions.
         /// \return  Array of size numPhases().
-        const double* surfaceDensities() const;
+        const double* surfaceDensities(int regionIdx = 0) const;
 
         /// Viscosity as a function of p and z.
         void mu(const int n,
+                const int *pvtTableIdx,
                 const double* p,
                 const double* z,
                 double* output_mu) const;
 
         /// Formation volume factor as a function of p and z.
         void B(const int n,
+               const int *pvtTableIdx,
                const double* p,
                const double* z,
                double* output_B) const;
 
         /// Formation volume factor and p-derivative as functions of p and z.
         void dBdp(const int n,
+                  const int *pvtTableIdx,
                   const double* p,
                   const double* z,
                   double* output_B,
@@ -91,12 +100,14 @@ namespace Opm
 
         /// Solution factor as a function of p and z.
         void R(const int n,
+               const int *pvtTableIdx,
                const double* p,
                const double* z,
                double* output_R) const;
 
         /// Solution factor and p-derivative as functions of p and z.
         void dRdp(const int n,
+                  const int *pvtTableIdx,
                   const double* p,
                   const double* z,
                   double* output_R,
@@ -109,11 +120,11 @@ namespace Opm
 
         PhaseUsage phase_usage_;
 
-        int region_number_;
+        // The PVT properties. We need to store one object per PVT
+        // region per active fluid phase.
+        std::vector<std::shared_ptr<PvtInterface> > props_;
+        std::vector<std::array<double, MaxNumPhases> > densities_;
 
-        std::vector<std::shared_ptr<SinglePvtInterface> > props_;
-
-        double densities_[MaxNumPhases];
         mutable std::vector<double> data1_;
         mutable std::vector<double> data2_;
     };
diff --git a/opm/core/props/pvt/PvtConstCompr.hpp b/opm/core/props/pvt/PvtConstCompr.hpp
new file mode 100644
index 00000000..fff6b83d
--- /dev/null
+++ b/opm/core/props/pvt/PvtConstCompr.hpp
@@ -0,0 +1,307 @@
+/*
+  Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_PVTCONSTCOMPR_HEADER_INCLUDED
+#define OPM_PVTCONSTCOMPR_HEADER_INCLUDED
+
+#include <opm/core/props/pvt/PvtInterface.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+
+#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
+#include <opm/parser/eclipse/Utility/PvcdoTable.hpp>
+
+#include <vector>
+#include <algorithm>
+
+
+namespace Opm
+{
+
+    /// Class for constant compressible phases (PVTW or PVCDO).
+    /// The PVT properties can either be given as a function of
+    /// pressure (p) and surface volume (z) or pressure (p) and gas
+    /// resolution factor (r). Also, since this class supports
+    /// multiple PVT regions, the concrete table to be used for each
+    /// data point needs to be specified via the pvtTableIdx argument
+    /// of the respective method. For all the virtual methods, the
+    /// following apply: pvtTableIdx, p, r and z are expected to be of
+    /// size n, size n, size n and n*num_phases, respectively.  Output
+    /// arrays shall be of size n, and must be valid before calling
+    /// the method.
+    class PvtConstCompr : public PvtInterface
+    {
+    public:
+        PvtConstCompr()
+        {}
+
+        void initFromWater(Opm::DeckKeywordConstPtr pvtwKeyword)
+        {
+            int numRegions = pvtwKeyword->size();
+
+            ref_press_.resize(numRegions);
+            ref_B_.resize(numRegions);
+            comp_.resize(numRegions);
+            viscosity_.resize(numRegions);
+            visc_comp_.resize(numRegions);
+
+            for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
+                Opm::DeckRecordConstPtr pvtwRecord = pvtwKeyword->getRecord(regionIdx);
+
+                ref_press_[regionIdx] = pvtwRecord->getItem("P_REF")->getSIDouble(0);
+                ref_B_[regionIdx]     = pvtwRecord->getItem("WATER_VOL_FACTOR")->getSIDouble(0);
+                comp_[regionIdx]      = pvtwRecord->getItem("WATER_COMPRESSIBILITY")->getSIDouble(0);
+                viscosity_[regionIdx] = pvtwRecord->getItem("WATER_VISCOSITY")->getSIDouble(0);
+                visc_comp_[regionIdx] = pvtwRecord->getItem("WATER_VISCOSIBILITY")->getSIDouble(0);
+            }
+        }
+
+        void initFromOil(Opm::DeckKeywordConstPtr pvcdoKeyword)
+        {
+            int numRegions = pvcdoKeyword->size();
+
+            ref_press_.resize(numRegions);
+            ref_B_.resize(numRegions);
+            comp_.resize(numRegions);
+            viscosity_.resize(numRegions);
+            visc_comp_.resize(numRegions);
+
+            for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
+                Opm::DeckRecordConstPtr pvcdoRecord = pvcdoKeyword->getRecord(regionIdx);
+
+                ref_press_[regionIdx] = pvcdoRecord->getItem("P_REF")->getSIDouble(0);
+                ref_B_[regionIdx]     = pvcdoRecord->getItem("OIL_VOL_FACTOR")->getSIDouble(0);
+                comp_[regionIdx]      = pvcdoRecord->getItem("OIL_COMPRESSIBILITY")->getSIDouble(0);
+                viscosity_[regionIdx] = pvcdoRecord->getItem("OIL_VISCOSITY")->getSIDouble(0);
+                visc_comp_[regionIdx] = pvcdoRecord->getItem("OIL_VISCOSIBILITY")->getSIDouble(0);
+            }
+        }
+
+        /*!
+         * \brief Create a PVT object with a given viscosity that
+         *        assumes all fluid phases to be incompressible.
+         */
+        explicit PvtConstCompr(double visc)
+            : ref_press_(1, 0.0),
+              ref_B_(1, 1.0),
+              comp_(1, 0.0),
+              viscosity_(1, visc),
+              visc_comp_(1, 0.0)
+        {
+        }
+
+        virtual ~PvtConstCompr()
+        {
+        }
+
+        virtual void mu(const int n,
+                        const int* pvtRegionIdx,
+                        const double* p,
+                        const double* /*z*/,
+                        double* output_mu) const
+        {
+// #pragma omp parallel for
+            for (int i = 0; i < n; ++i) {
+                // Computing a polynomial approximation to the exponential.
+                int tableIdx = getTableIndex_(pvtRegionIdx, i);
+                double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
+                output_mu[i] = viscosity_[tableIdx]/(1.0 + x + 0.5*x*x);
+            }
+        }
+
+        virtual void mu(const int n,
+                        const int* pvtRegionIdx,
+                        const double* p,
+                        const double* /*r*/,
+                        double* output_mu,
+                        double* output_dmudp,
+                        double* output_dmudr) const
+        {
+            // #pragma omp parallel for
+            for (int i = 0; i < n; ++i) {
+                // Computing a polynomial approximation to the exponential.
+                int tableIdx = getTableIndex_(pvtRegionIdx, i);
+                double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
+                double d = (1.0 + x + 0.5*x*x);
+                output_mu[i] = viscosity_[tableIdx]/d;
+                output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx];
+            }
+            std::fill(output_dmudr, output_dmudr + n, 0.0);
+        }
+
+        virtual void mu(const int n,
+                        const int* pvtRegionIdx,
+                        const double* p,
+                        const double* /*r*/,
+                        const PhasePresence* /*cond*/,
+                        double* output_mu,
+                        double* output_dmudp,
+                        double* output_dmudr) const
+        {
+            // #pragma omp parallel for
+            for (int i = 0; i < n; ++i) {
+                // Computing a polynomial approximation to the exponential.
+                int tableIdx = getTableIndex_(pvtRegionIdx, i);
+                double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
+                double d = (1.0 + x + 0.5*x*x);
+                output_mu[i] = viscosity_[tableIdx]/d;
+                output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx];
+            }
+            std::fill(output_dmudr, output_dmudr + n, 0.0);
+        }
+
+        virtual void B(const int n,
+                       const int* pvtRegionIdx,
+                       const double* p,
+                       const double* /*z*/,
+                       double* output_B) const
+        {
+// #pragma omp parallel for
+            for (int i = 0; i < n; ++i) {
+                // Computing a polynomial approximation to the exponential.
+                int tableIdx = getTableIndex_(pvtRegionIdx, i);
+                double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
+                output_B[i] = ref_B_[tableIdx]/(1.0 + x + 0.5*x*x);
+            }
+        }
+
+        virtual void dBdp(const int n,
+                          const int* pvtRegionIdx,
+                          const double* p,
+                          const double* /*z*/,
+                          double* output_B,
+                          double* output_dBdp) const
+        {
+// #pragma omp parallel for
+            for (int i = 0; i < n; ++i) {
+                int tableIdx = getTableIndex_(pvtRegionIdx, i);
+                double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
+                double d = (1.0 + x + 0.5*x*x);
+                output_B[i] = ref_B_[tableIdx]/d;
+                output_dBdp[i] = (-ref_B_[tableIdx]/(d*d))*(1 + x) * comp_[tableIdx];
+            }
+        }
+
+        virtual void b(const int n,
+                       const int* pvtRegionIdx,
+                       const double* p,
+                       const double* /*r*/,
+                       double* output_b,
+                       double* output_dbdp,
+                       double* output_dbdr) const
+        {
+        // #pragma omp parallel for
+            for (int i = 0; i < n; ++i) {
+                // Computing a polynomial approximation to the exponential.
+                int tableIdx = getTableIndex_(pvtRegionIdx, i);
+                double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
+                double d = (1.0 + x + 0.5*x*x);
+
+                // b = 1/B = d/ref_B_B;
+                output_b[i] = d/ref_B_[tableIdx];
+                output_dbdp[i] = (1 + x) * comp_[tableIdx]/ref_B_[tableIdx];
+            }
+
+            std::fill(output_dbdr, output_dbdr + n, 0.0);
+        }
+
+        virtual void b(const int n,
+                       const int* pvtRegionIdx,
+                       const double* p,
+                       const double* /*r*/,
+                       const PhasePresence* /*cond*/,
+                       double* output_b,
+                       double* output_dbdp,
+                       double* output_dbdr) const
+        {
+        // #pragma omp parallel for
+            for (int i = 0; i < n; ++i) {
+                // Computing a polynomial approximation to the exponential.
+                int tableIdx = getTableIndex_(pvtRegionIdx, i);
+                double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
+                double d = (1.0 + x + 0.5*x*x);
+
+                // b = 1/B = d/ref_B_[tableIdx]B;
+                output_b[i] = d/ref_B_[tableIdx];
+                output_dbdp[i] = (1 + x) * comp_[tableIdx]/ref_B_[tableIdx];
+            }
+            std::fill(output_dbdr, output_dbdr + n, 0.0);
+        }
+
+        virtual void rsSat(const int n,
+                           const int* /*pvtRegionIdx*/,
+                           const double* /*p*/,
+                           double* output_rsSat,
+                           double* output_drsSatdp) const
+        {
+            std::fill(output_rsSat, output_rsSat + n, 0.0);
+            std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
+        }
+
+        virtual void rvSat(const int n,
+                           const int* /*pvtRegionIdx*/,
+                           const double* /*p*/,
+                           double* output_rvSat,
+                           double* output_drvSatdp) const
+        {
+            std::fill(output_rvSat, output_rvSat + n, 0.0);
+            std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
+        }
+
+        virtual void R(const int n,
+                       const int* /*pvtRegionIdx*/,
+                       const double* /*p*/,
+                       const double* /*z*/,
+                       double* output_R) const
+        {
+            std::fill(output_R, output_R + n, 0.0);
+        }
+
+        virtual void dRdp(const int n,
+                          const int* /*pvtRegionIdx*/,
+                          const double* /*p*/,
+                          const double* /*z*/,
+                          double* output_R,
+                          double* output_dRdp) const
+        {
+            std::fill(output_R, output_R + n, 0.0);
+            std::fill(output_dRdp, output_dRdp + n, 0.0);
+        }
+
+    private:
+        int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
+        {
+            if (!pvtTableIdx)
+                return 0;
+            return pvtTableIdx[cellIdx];
+        }
+
+        // The PVT properties. We need to store one value per PVT
+        // region.
+        std::vector<double> ref_press_;
+        std::vector<double> ref_B_;
+        std::vector<double> comp_;
+        std::vector<double> viscosity_;
+        std::vector<double> visc_comp_;
+    };
+
+}
+
+
+#endif // OPM_PVTCONSTCOMPR_HEADER_INCLUDED
+
diff --git a/opm/core/props/pvt/SinglePvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp
similarity index 55%
rename from opm/core/props/pvt/SinglePvtDead.cpp
rename to opm/core/props/pvt/PvtDead.cpp
index 68449de2..09dc32d3 100644
--- a/opm/core/props/pvt/SinglePvtDead.cpp
+++ b/opm/core/props/pvt/PvtDead.cpp
@@ -19,7 +19,7 @@
 
 
 #include "config.h"
-#include <opm/core/props/pvt/SinglePvtDead.hpp>
+#include <opm/core/props/pvt/PvtDead.hpp>
 #include <algorithm>
 
 // Extra includes for debug dumping of tables.
@@ -34,72 +34,80 @@ namespace Opm
     // Member functions
     //-------------------------------------------------------------------------
     /// Constructor
-    SinglePvtDead::SinglePvtDead(const Opm::PvdoTable& pvdoTable)
+    void PvtDead::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword)
     {
-        // Copy data
-        const std::vector<double>& press = pvdoTable.getPressureColumn();
-        const std::vector<double>& b = pvdoTable.getFormationFactorColumn();
-        const std::vector<double>& visc = pvdoTable.getViscosityColumn();
+        int numRegions = Opm::PvdoTable::numTables(pvdoKeyword);
 
-        const int sz = b.size();
-        std::vector<double> bInv(sz);
-        for (int i = 0; i < sz; ++i) {
-            bInv[i] = 1.0 / b[i];
+        // resize the attributes of the object
+        b_.resize(numRegions);
+        viscosity_.resize(numRegions);
+
+        for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
+            Opm::PvdoTable pvdoTable(pvdoKeyword, regionIdx);
+            // Copy data
+            const std::vector<double>& press = pvdoTable.getPressureColumn();
+            const std::vector<double>& b = pvdoTable.getFormationFactorColumn();
+            const std::vector<double>& visc = pvdoTable.getViscosityColumn();
+
+            const int sz = b.size();
+            std::vector<double> bInv(sz);
+            for (int i = 0; i < sz; ++i) {
+                bInv[i] = 1.0 / b[i];
+            }
+            b_[regionIdx] = NonuniformTableLinear<double>(press, bInv);
+            viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
         }
-        b_ = NonuniformTableLinear<double>(press, bInv);
-        viscosity_ = NonuniformTableLinear<double>(press, visc);
-
-        // Dumping the created tables.
-//         static int count = 0;
-//         std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
-//         os.precision(15);
-//         os << "1/B\n\n" << one_over_B_
-//            << "\n\nvisc\n\n" << viscosity_ << std::endl;
     }
 
-    /// Constructor
-    SinglePvtDead::SinglePvtDead(const Opm::PvdgTable& pvdgTable)
+
+    void PvtDead::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword)
     {
-        // Copy data
-        const std::vector<double>& press = pvdgTable.getPressureColumn();
-        const std::vector<double>& b = pvdgTable.getFormationFactorColumn();
-        const std::vector<double>& visc = pvdgTable.getViscosityColumn();
+        int numRegions = Opm::PvdgTable::numTables(pvdgKeyword);
 
-        const int sz = b.size();
-        std::vector<double> bInv(sz);
-        for (int i = 0; i < sz; ++i) {
-            bInv[i] = 1.0 / b[i];
+        // resize the attributes of the object
+        b_.resize(numRegions);
+        viscosity_.resize(numRegions);
+
+        for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
+            Opm::PvdgTable pvdgTable(pvdgKeyword, regionIdx);
+
+            // Copy data
+            const std::vector<double>& press = pvdgTable.getPressureColumn();
+            const std::vector<double>& b = pvdgTable.getFormationFactorColumn();
+            const std::vector<double>& visc = pvdgTable.getViscosityColumn();
+
+            const int sz = b.size();
+            std::vector<double> bInv(sz);
+            for (int i = 0; i < sz; ++i) {
+                bInv[i] = 1.0 / b[i];
+            }
+            b_[regionIdx] = NonuniformTableLinear<double>(press, bInv);
+            viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
         }
-        b_ = NonuniformTableLinear<double>(press, bInv);
-        viscosity_ = NonuniformTableLinear<double>(press, visc);
-
-        // Dumping the created tables.
-//         static int count = 0;
-//         std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
-//         os.precision(15);
-//         os << "1/B\n\n" << one_over_B_
-//            << "\n\nvisc\n\n" << viscosity_ << std::endl;
     }
 
     // Destructor
-    SinglePvtDead::~SinglePvtDead()
+    PvtDead::~PvtDead()
     {
     }
 
 
 
-    void SinglePvtDead::mu(const int n,
+    void PvtDead::mu(const int n,
+                     const int* pvtTableIdx,
                            const double* p,
                            const double* /*z*/,
                            double* output_mu) const
     {
 // #pragma omp parallel for
         for (int i = 0; i < n; ++i) {
-            output_mu[i] = viscosity_(p[i]);
+            int regionIdx = getTableIndex_(pvtTableIdx, i);
+            output_mu[i] = viscosity_[regionIdx](p[i]);
         }
     }
 
-    void SinglePvtDead::mu(const int n,
+    void PvtDead::mu(const int n,
+                     const int* pvtTableIdx,
                                const double* p,
                                const double* /*r*/,
                                double* output_mu,
@@ -108,14 +116,16 @@ namespace Opm
         {
     // #pragma omp parallel for
             for (int i = 0; i < n; ++i) {
-                output_mu[i] = viscosity_(p[i]);
-                output_dmudp[i] = viscosity_.derivative(p[i]);
+                int regionIdx = getTableIndex_(pvtTableIdx, i);
+                output_mu[i] = viscosity_[regionIdx](p[i]);
+                output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
             }
             std::fill(output_dmudr, output_dmudr + n, 0.0);
 
         }
 
-    void SinglePvtDead::mu(const int n,
+    void PvtDead::mu(const int n,
+                     const int* pvtTableIdx,
                                const double* p,
                                const double* /*r*/,
                                const PhasePresence* /*cond*/,
@@ -125,14 +135,16 @@ namespace Opm
         {
     // #pragma omp parallel for
             for (int i = 0; i < n; ++i) {
-                output_mu[i] = viscosity_(p[i]);
-                output_dmudp[i] = viscosity_.derivative(p[i]);
+                int regionIdx = getTableIndex_(pvtTableIdx, i);
+                output_mu[i] = viscosity_[regionIdx](p[i]);
+                output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
             }
             std::fill(output_dmudr, output_dmudr + n, 0.0);
 
         }
 
-    void SinglePvtDead::B(const int n,
+    void PvtDead::B(const int n,
+                    const int* pvtTableIdx,
                           const double* p,
                           const double* /*z*/,
                           double* output_B) const
@@ -140,25 +152,29 @@ namespace Opm
 // #pragma omp parallel for
         // B = 1/b
         for (int i = 0; i < n; ++i) {
-            output_B[i] = 1.0/b_(p[i]);
+            int regionIdx = getTableIndex_(pvtTableIdx, i);
+            output_B[i] = 1.0/b_[regionIdx](p[i]);
         }
     }
 
-    void SinglePvtDead::dBdp(const int n,
+    void PvtDead::dBdp(const int n,
+                       const int* pvtTableIdx,
                              const double* p,
                              const double* /*z*/,
                              double* output_B,
                              double* output_dBdp) const
     {
-        B(n, p, 0, output_B);
+        B(n, pvtTableIdx, p, 0, output_B);
 // #pragma omp parallel for
         for (int i = 0; i < n; ++i) {
+            int regionIdx = getTableIndex_(pvtTableIdx, i);
             double Bg = output_B[i];
-            output_dBdp[i] = -Bg*Bg*b_.derivative(p[i]);
+            output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]);
         }
     }
 
-    void SinglePvtDead::b(const int n,
+    void PvtDead::b(const int n,
+                    const int* pvtTableIdx,
                               const double* p,
                               const double* /*r*/,
                               double* output_b,
@@ -168,15 +184,18 @@ namespace Opm
         {
     // #pragma omp parallel for
             for (int i = 0; i < n; ++i) {
-                output_b[i] = b_(p[i]);
-                output_dbdp[i] = b_.derivative(p[i]);
+                int regionIdx = getTableIndex_(pvtTableIdx, i);
+
+                output_b[i] = b_[regionIdx](p[i]);
+                output_dbdp[i] = b_[regionIdx].derivative(p[i]);
 
             }
             std::fill(output_dbdr, output_dbdr + n, 0.0);
 
         }
 
-    void SinglePvtDead::b(const int n,
+    void PvtDead::b(const int n,
+                    const int* pvtTableIdx,
                               const double* p,
                               const double* /*r*/,
                               const PhasePresence* /*cond*/,
@@ -187,15 +206,18 @@ namespace Opm
         {
     // #pragma omp parallel for
             for (int i = 0; i < n; ++i) {
-                output_b[i] = b_(p[i]);
-                output_dbdp[i] = b_.derivative(p[i]);
+                int regionIdx = getTableIndex_(pvtTableIdx, i);
+
+                output_b[i] = b_[regionIdx](p[i]);
+                output_dbdp[i] = b_[regionIdx].derivative(p[i]);
 
             }
             std::fill(output_dbdr, output_dbdr + n, 0.0);
 
         }
 
-    void SinglePvtDead::rsSat(const int n,
+    void PvtDead::rsSat(const int n,
+                        const int* pvtTableIdx,
                              const double* /*p*/,
                              double* output_rsSat,
                              double* output_drsSatdp) const
@@ -204,7 +226,8 @@ namespace Opm
         std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
     }
 
-    void SinglePvtDead::rvSat(const int n,
+    void PvtDead::rvSat(const int n,
+                        const int* pvtTableIdx,
                              const double* /*p*/,
                              double* output_rvSat,
                              double* output_drvSatdp) const
@@ -213,7 +236,8 @@ namespace Opm
         std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
     }
 
-    void SinglePvtDead::R(const int n,
+    void PvtDead::R(const int n,
+                    const int* pvtTableIdx,
                           const double* /*p*/,
                           const double* /*z*/,
                           double* output_R) const
@@ -221,7 +245,8 @@ namespace Opm
         std::fill(output_R, output_R + n, 0.0);
     }
 
-    void SinglePvtDead::dRdp(const int n,
+    void PvtDead::dRdp(const int n,
+                       const int* pvtTableIdx,
                              const double* /*p*/,
                              const double* /*z*/,
                              double* output_R,
diff --git a/opm/core/props/pvt/SinglePvtDead.hpp b/opm/core/props/pvt/PvtDead.hpp
similarity index 69%
rename from opm/core/props/pvt/SinglePvtDead.hpp
rename to opm/core/props/pvt/PvtDead.hpp
index e7d28872..3d01eb2a 100644
--- a/opm/core/props/pvt/SinglePvtDead.hpp
+++ b/opm/core/props/pvt/PvtDead.hpp
@@ -17,11 +17,11 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#ifndef OPM_SINGLEPVTDEAD_HEADER_INCLUDED
-#define OPM_SINGLEPVTDEAD_HEADER_INCLUDED
+#ifndef OPM_PVTDEAD_HEADER_INCLUDED
+#define OPM_PVTDEAD_HEADER_INCLUDED
 
 
-#include <opm/core/props/pvt/SinglePvtInterface.hpp>
+#include <opm/core/props/pvt/PvtInterface.hpp>
 #include <opm/core/utility/NonuniformTableLinear.hpp>
 
 #include <opm/parser/eclipse/Utility/PvdoTable.hpp>
@@ -33,21 +33,28 @@ namespace Opm
 {
 
     /// Class for immiscible dead oil and dry gas.
-    /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
-    /// or pressure (p) and gas resolution factor (r).
-    /// For all the virtual methods, the following apply: p, r and z
-    /// are expected to be of size n, size n and n*num_phases, respectively.
-    /// Output arrays shall be of size n, and must be valid before
-    /// calling the method.
-    class SinglePvtDead : public SinglePvtInterface
+    /// The PVT properties can either be given as a function of
+    /// pressure (p) and surface volume (z) or pressure (p) and gas
+    /// resolution factor (r). Also, since this class supports
+    /// multiple PVT regions, the concrete table to be used for each
+    /// data point needs to be specified via the pvtTableIdx argument
+    /// of the respective method. For all the virtual methods, the
+    /// following apply: pvtTableIdx, p, r and z are expected to be of
+    /// size n, size n, size n and n*num_phases, respectively.  Output
+    /// arrays shall be of size n, and must be valid before calling
+    /// the method.
+    class PvtDead : public PvtInterface
     {
     public:
-        SinglePvtDead(const Opm::PvdoTable& pvdoTable);
-        SinglePvtDead(const Opm::PvdgTable& pvdgTable);
-        virtual ~SinglePvtDead();
+        PvtDead() {};
+
+        void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword);
+        void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword);
+        virtual ~PvtDead();
 
         /// Viscosity as a function of p and z.
         virtual void mu(const int n,
+                        const int* pvtTableIdx,
                         const double* p,
                         const double* z,
                         double* output_mu) const;
@@ -55,6 +62,7 @@ namespace Opm
         /// Viscosity and its derivatives as a function of p and r.
         /// The fluid is considered saturated if r >= rsSat(p).
         virtual void mu(const int n,
+                        const int* pvtTableIdx,
                         const double* p,
                         const double* r,
                         double* output_mu,
@@ -64,6 +72,7 @@ namespace Opm
         /// Viscosity as a function of p and r.
         /// State condition determined by 'cond'.
         virtual void mu(const int n,
+                        const int* pvtTableIdx,
                         const double* p,
                         const double* r,
                         const PhasePresence* cond,
@@ -73,12 +82,14 @@ namespace Opm
 
         /// Formation volume factor as a function of p and z.
         virtual void B(const int n,
+                       const int* pvtTableIdx,
                        const double* p,
                        const double* z,
                        double* output_B) const;
 
         /// Formation volume factor and p-derivative as functions of p and z.
         virtual void dBdp(const int n,
+                          const int* pvtTableIdx,
                           const double* p,
                           const double* z,
                           double* output_B,
@@ -87,6 +98,7 @@ namespace Opm
         /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
         /// The fluid is considered saturated if r >= rsSat(p).
         virtual void b(const int n,
+                       const int* pvtTableIdx,
                        const double* p,
                        const double* r,
                        double* output_b,
@@ -96,6 +108,7 @@ namespace Opm
         /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
         /// State condition determined by 'cond'.
         virtual void b(const int n,
+                       const int* pvtTableIdx,
                        const double* p,
                        const double* r,
                        const PhasePresence* cond,
@@ -106,35 +119,46 @@ namespace Opm
 
         /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
         virtual void rsSat(const int n,
+                           const int* pvtTableIdx,
                           const double* p,
                           double* output_rsSat,
                           double* output_drsSatdp) const;
 
         /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
         virtual void rvSat(const int n,
+                           const int* pvtTableIdx,
                           const double* p,
                           double* output_rvSat,
                           double* output_drvSatdp) const;
 
         /// Solution factor as a function of p and z.
         virtual void R(const int n,
+                       const int* pvtTableIdx,
                        const double* p,
                        const double* z,
                        double* output_R) const;
 
         /// Solution factor and p-derivative as functions of p and z.
         virtual void dRdp(const int n,
+                          const int* pvtTableIdx,
                           const double* p,
                           const double* z,
                           double* output_R,
                           double* output_dRdp) const;
     private:
-        // PVT properties of dry gas or dead oil
-        NonuniformTableLinear<double> b_;
-        NonuniformTableLinear<double> viscosity_;
-    };
+        int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
+        {
+            if (!pvtTableIdx)
+                return 0;
+            return pvtTableIdx[cellIdx];
+        }
 
+        // PVT properties of dry gas or dead oil. We need to store one
+        // table per PVT region.
+        std::vector<NonuniformTableLinear<double> > b_;
+        std::vector<NonuniformTableLinear<double> > viscosity_;
+    };
 }
 
 
-#endif // OPM_SINGLEPVTDEAD_HEADER_INCLUDED
+#endif // OPM_PVTDEAD_HEADER_INCLUDED
diff --git a/opm/core/props/pvt/PvtDeadSpline.cpp b/opm/core/props/pvt/PvtDeadSpline.cpp
new file mode 100644
index 00000000..69f749b7
--- /dev/null
+++ b/opm/core/props/pvt/PvtDeadSpline.cpp
@@ -0,0 +1,260 @@
+/*
+  Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "config.h"
+#include <opm/core/props/pvt/PvtDeadSpline.hpp>
+#include <opm/core/utility/buildUniformMonotoneTable.hpp>
+
+#include <opm/parser/eclipse/Utility/SingleRecordTable.hpp>
+
+#include <algorithm>
+
+// Extra includes for debug dumping of tables.
+// #include <boost/lexical_cast.hpp>
+// #include <string>
+// #include <fstream>
+
+namespace Opm
+{
+    //------------------------------------------------------------------------
+    // Member functions
+    //-------------------------------------------------------------------------
+
+    PvtDeadSpline::PvtDeadSpline()
+    {}
+
+    void PvtDeadSpline::initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword,
+                                    int numSamples)
+    {
+        int numRegions = Opm::PvdoTable::numTables(pvdoKeyword);
+
+        // resize the attributes of the object
+        b_.resize(numRegions);
+        viscosity_.resize(numRegions);
+
+        for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
+            Opm::PvdoTable pvdoTable(pvdoKeyword, regionIdx);
+
+            int numRows = pvdoTable.numRows();
+
+            // Copy data
+            const std::vector<double> &press = pvdoTable.getPressureColumn();
+            const std::vector<double> &B = pvdoTable.getFormationFactorColumn();
+            const std::vector<double> &visc = pvdoTable.getViscosityColumn();
+
+            std::vector<double> B_inv(numRows);
+            for (int i = 0; i < numRows; ++i) {
+                B_inv[i] = 1.0 / B[i];
+            }
+
+            buildUniformMonotoneTable(press, B_inv, numSamples, b_[regionIdx]);
+            buildUniformMonotoneTable(press, visc, numSamples, viscosity_[regionIdx]);
+        }
+    }
+
+    void PvtDeadSpline::initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword,
+                                    int numSamples)
+    {
+        int numRegions = Opm::PvdgTable::numTables(pvdgKeyword);
+
+        // resize the attributes of the object
+        b_.resize(numRegions);
+        viscosity_.resize(numRegions);
+
+        for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
+            Opm::PvdgTable pvdgTable(pvdgKeyword, regionIdx);
+
+            int numRows = pvdgTable.numRows();
+
+            // Copy data
+            const std::vector<double> &press = pvdgTable.getPressureColumn();
+            const std::vector<double> &B = pvdgTable.getFormationFactorColumn();
+            const std::vector<double> &visc = pvdgTable.getViscosityColumn();
+
+            std::vector<double> B_inv(numRows);
+            for (int i = 0; i < numRows; ++i) {
+                B_inv[i] = 1.0 / B[i];
+            }
+
+            buildUniformMonotoneTable(press, B_inv, numSamples, b_[regionIdx]);
+            buildUniformMonotoneTable(press, visc, numSamples, viscosity_[regionIdx]);
+        }
+    }
+
+    // Destructor
+    PvtDeadSpline::~PvtDeadSpline()
+    {
+    }
+
+
+
+    void PvtDeadSpline::mu(const int n,
+                           const int* pvtTableIdx,
+                                 const double* p,
+                                 const double* /*z*/,
+                                 double* output_mu) const
+    {
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            int regionIdx = getTableIndex_(pvtTableIdx, i);
+            output_mu[i] = viscosity_[regionIdx](p[i]);
+        }
+    }
+
+    void PvtDeadSpline::mu(const int n,
+                           const int* pvtTableIdx,
+                           const double* p,
+                           const double* /*r*/,
+                           double* output_mu,
+                           double* output_dmudp,
+                           double* output_dmudr) const
+    {
+        // #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            int regionIdx = getTableIndex_(pvtTableIdx, i);
+            output_mu[i] = viscosity_[regionIdx](p[i]);
+            output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
+        }
+        std::fill(output_dmudr, output_dmudr + n, 0.0);
+    }
+
+    void PvtDeadSpline::mu(const int n,
+                           const int* pvtTableIdx,
+                           const double* p,
+                           const double* /*r*/,
+                           const PhasePresence* /*cond*/,
+                           double* output_mu,
+                           double* output_dmudp,
+                           double* output_dmudr) const
+    {
+    // #pragma omp parallel for
+
+        for (int i = 0; i < n; ++i) {
+            int regionIdx = getTableIndex_(pvtTableIdx, i);
+            output_mu[i] = viscosity_[regionIdx](p[i]);
+            output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
+        }
+        std::fill(output_dmudr, output_dmudr + n, 0.0);
+    }
+
+    void PvtDeadSpline::B(const int n,
+                          const int* pvtTableIdx,
+                                const double* p,
+                                const double* /*z*/,
+                                double* output_B) const
+    {
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            int regionIdx = getTableIndex_(pvtTableIdx, i);
+            output_B[i] = 1.0/b_[regionIdx](p[i]);
+        }
+    }
+
+    void PvtDeadSpline::dBdp(const int n,
+                             const int* pvtTableIdx,
+                             const double* p,
+                             const double* /*z*/,
+                             double* output_B,
+                             double* output_dBdp) const
+    {
+        B(n, pvtTableIdx, p, 0, output_B);
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            int regionIdx = getTableIndex_(pvtTableIdx, i);
+            double Bg = output_B[i];
+            output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]);
+        }
+    }
+
+    void PvtDeadSpline::b(const int n,
+                          const int* pvtTableIdx,
+                              const double* p,
+                              const double* /*r*/,
+                              double* output_b,
+                              double* output_dbdp,
+                              double* output_dbdr) const
+    {
+    // #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            int regionIdx = getTableIndex_(pvtTableIdx, i);
+            output_b[i] = b_[regionIdx](p[i]);
+            output_dbdp[i] = b_[regionIdx].derivative(p[i]);
+        }
+        std::fill(output_dbdr, output_dbdr + n, 0.0);
+    }
+
+    void PvtDeadSpline::b(const int n,
+                          const int* pvtTableIdx,
+                          const double* p,
+                          const double* /*r*/,
+                          const PhasePresence* /*cond*/,
+                          double* output_b,
+                          double* output_dbdp,
+                          double* output_dbdr) const
+    {
+    // #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            int regionIdx = getTableIndex_(pvtTableIdx, i);
+            output_b[i] = b_[regionIdx](p[i]);
+            output_dbdp[i] = b_[regionIdx].derivative(p[i]);
+        }
+        std::fill(output_dbdr, output_dbdr + n, 0.0);
+    }
+
+    void PvtDeadSpline::rsSat(const int n,
+                              const int* /*pvtTableIdx*/,
+                             const double* /*p*/,
+                             double* output_rsSat,
+                             double* output_drsSatdp) const
+    {
+        std::fill(output_rsSat, output_rsSat + n, 0.0);
+        std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
+    }
+
+    void PvtDeadSpline::rvSat(const int n,
+                              const int* /*pvtTableIdx*/,
+                             const double* /*p*/,
+                             double* output_rvSat,
+                             double* output_drvSatdp) const
+    {
+        std::fill(output_rvSat, output_rvSat + n, 0.0);
+        std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
+    }
+
+    void PvtDeadSpline::R(const int n,
+                          const int* /*pvtTableIdx*/,
+                          const double* /*p*/,
+                          const double* /*z*/,
+                          double* output_R) const
+    {
+        std::fill(output_R, output_R + n, 0.0);
+    }
+
+    void PvtDeadSpline::dRdp(const int n,
+                             const int* /*pvtTableIdx*/,
+                             const double* /*p*/,
+                             const double* /*z*/,
+                             double* output_R,
+                             double* output_dRdp) const
+    {
+        std::fill(output_R, output_R + n, 0.0);
+        std::fill(output_dRdp, output_dRdp + n, 0.0);
+    }
+
+}
diff --git a/opm/core/props/pvt/SinglePvtDeadSpline.hpp b/opm/core/props/pvt/PvtDeadSpline.hpp
similarity index 77%
rename from opm/core/props/pvt/SinglePvtDeadSpline.hpp
rename to opm/core/props/pvt/PvtDeadSpline.hpp
index 88320ed2..b7c8beba 100644
--- a/opm/core/props/pvt/SinglePvtDeadSpline.hpp
+++ b/opm/core/props/pvt/PvtDeadSpline.hpp
@@ -17,10 +17,10 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#ifndef OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
-#define OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
+#ifndef OPM_PVTDEADSPLINE_HEADER_INCLUDED
+#define OPM_PVTDEADSPLINE_HEADER_INCLUDED
 
-#include <opm/core/props/pvt/SinglePvtInterface.hpp>
+#include <opm/core/props/pvt/PvtInterface.hpp>
 #include <opm/core/utility/UniformTableLinear.hpp>
 
 #include <opm/parser/eclipse/Utility/PvdoTable.hpp>
@@ -38,15 +38,21 @@ namespace Opm
     /// are expected to be of size n, size n and n*num_phases, respectively.
     /// Output arrays shall be of size n, and must be valid before
     /// calling the method.
-    class SinglePvtDeadSpline : public SinglePvtInterface
+    class PvtDeadSpline : public PvtInterface
     {
     public:
-        SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples);
-        SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples);
-        virtual ~SinglePvtDeadSpline();
+        PvtDeadSpline();
+
+        void initFromOil(Opm::DeckKeywordConstPtr pvdoKeyword,
+                         int numSamples);
+        void initFromGas(Opm::DeckKeywordConstPtr pvdgKeyword,
+                         int numSamples);
+
+        virtual ~PvtDeadSpline();
 
         /// Viscosity as a function of p and z.
         virtual void mu(const int n,
+                        const int* pvtTableIdx,
                         const double* p,
                         const double* z,
                         double* output_mu) const;
@@ -54,6 +60,7 @@ namespace Opm
         /// Viscosity and its derivatives as a function of p and r.
         /// The fluid is considered saturated if r >= rsSat(p).
         virtual void mu(const int n,
+                        const int* pvtTableIdx,
                         const double* p,
                         const double* r,
                         double* output_mu,
@@ -63,6 +70,7 @@ namespace Opm
         /// Viscosity as a function of p and r.
         /// State condition determined by 'cond'.
         virtual void mu(const int n,
+                        const int* pvtTableIdx,
                         const double* p,
                         const double* r,
                         const PhasePresence* cond,
@@ -72,12 +80,14 @@ namespace Opm
 
         /// Formation volume factor as a function of p and z.
         virtual void B(const int n,
+                       const int* pvtTableIdx,
                        const double* p,
                        const double* z,
                        double* output_B) const;
 
         /// Formation volume factor and p-derivative as functions of p and z.
         virtual void dBdp(const int n,
+                          const int* pvtTableIdx,
                           const double* p,
                           const double* z,
                           double* output_B,
@@ -86,6 +96,7 @@ namespace Opm
         /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
         /// The fluid is considered saturated if r >= rsSat(p).
         virtual void b(const int n,
+                       const int* pvtTableIdx,
                        const double* p,
                        const double* r,
                        double* output_b,
@@ -95,6 +106,7 @@ namespace Opm
         /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
         /// State condition determined by 'cond'.
         virtual void b(const int n,
+                       const int* pvtTableIdx,
                        const double* p,
                        const double* r,
                        const PhasePresence* cond,
@@ -104,35 +116,47 @@ namespace Opm
 
         /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
         virtual void rsSat(const int n,
+                           const int* pvtTableIdx,
                           const double* p,
                           double* output_rsSat,
                           double* output_drsSatdp) const;
 
         /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
         virtual void rvSat(const int n,
+                           const int* pvtTableIdx,
                           const double* p,
                           double* output_rvSat,
                           double* output_drvSatdp) const;
 
         /// Solution factor as a function of p and z.
         virtual void R(const int n,
+                       const int* pvtTableIdx,
                        const double* p,
                        const double* z,
                        double* output_R) const;
 
         /// Solution factor and p-derivative as functions of p and z.
         virtual void dRdp(const int n,
+                          const int* pvtTableIdx,
                           const double* p,
                           const double* z,
                           double* output_R,
                           double* output_dRdp) const;
     private:
-        // PVT properties of dry gas or dead oil
-        UniformTableLinear<double> b_;
-        UniformTableLinear<double> viscosity_;
+        int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
+        {
+            if (!pvtTableIdx)
+                return 0;
+            return pvtTableIdx[cellIdx];
+        }
+
+        // PVT properties of dry gas or dead oil. We need to store one
+        // table per PVT region.
+        std::vector<UniformTableLinear<double> > b_;
+        std::vector<UniformTableLinear<double> > viscosity_;
     };
 
 }
 
-#endif // OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
+#endif // OPM_PVTDEADSPLINE_HEADER_INCLUDED
 
diff --git a/opm/core/props/pvt/PvtInterface.cpp b/opm/core/props/pvt/PvtInterface.cpp
new file mode 100644
index 00000000..105b0e91
--- /dev/null
+++ b/opm/core/props/pvt/PvtInterface.cpp
@@ -0,0 +1,79 @@
+/*
+  Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "config.h"
+#include <opm/core/props/pvt/PvtInterface.hpp>
+
+#include <cassert>
+
+namespace Opm
+{
+
+    PvtInterface::PvtInterface()
+        : num_phases_(MaxNumPhases)
+    {
+        for (int i = 0; i < MaxNumPhases; ++i) {
+            phase_pos_[i] = i;
+        }
+    }
+
+    PvtInterface::~PvtInterface()
+    {
+    }
+
+    void PvtInterface::setPhaseConfiguration(const int num_phases, const int* phase_pos)
+    {
+        num_phases_ = num_phases;
+        for (int i = 0; i < MaxNumPhases; ++i) {
+            phase_pos_[i] = phase_pos[i];
+        }
+    }
+
+    void extractPvtTableIndex(std::vector<int> &pvtTableIdx,
+                              Opm::DeckConstPtr deck,
+                              size_t numCompressed,
+                              const int *compressedToCartesianCellIdx)
+    {
+        if (deck->hasKeyword("PVTNUM")) {
+            // we have a "real" multi-pvt deck.
+
+            // get the PVT number from the deck
+            Opm::DeckKeywordConstPtr pvtnumKeyword = deck->getKeyword("PVTNUM");
+            const std::vector<int> &pvtnumData = pvtnumKeyword->getIntData();
+
+            // convert this into an array of compressed cells. since
+            // Eclipse uses Fortran-style indices which start at 1
+            // instead of 0, we subtract 1.
+            pvtTableIdx.resize(numCompressed);
+            for (size_t cellIdx = 0; cellIdx < numCompressed; ++ cellIdx) {
+                size_t cartesianCellIdx = compressedToCartesianCellIdx?compressedToCartesianCellIdx[cellIdx]:cellIdx;
+                assert(cartesianCellIdx < pvtnumData.size());
+
+                pvtTableIdx[cellIdx] = pvtnumData[cartesianCellIdx] - 1;
+            }
+        }
+        else {
+            // create the pvtIdxArray: all cells use the one and only
+            // PVT table...
+            pvtTableIdx.resize(numCompressed);
+            std::fill(pvtTableIdx.begin(), pvtTableIdx.end(), 0);
+        }
+    }
+
+} // namespace Opm
diff --git a/opm/core/props/pvt/SinglePvtInterface.hpp b/opm/core/props/pvt/PvtInterface.hpp
similarity index 69%
rename from opm/core/props/pvt/SinglePvtInterface.hpp
rename to opm/core/props/pvt/PvtInterface.hpp
index 83d3febb..85b5b74f 100644
--- a/opm/core/props/pvt/SinglePvtInterface.hpp
+++ b/opm/core/props/pvt/PvtInterface.hpp
@@ -17,22 +17,22 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#ifndef OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED
-#define OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED
-
+#ifndef OPM_PVTINTERFACE_HEADER_INCLUDED
+#define OPM_PVTINTERFACE_HEADER_INCLUDED
 
 #include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
 
 
 namespace Opm
 {
 
-    class SinglePvtInterface : public BlackoilPhases
+    class PvtInterface : public BlackoilPhases
     {
     public:
-        SinglePvtInterface();
+        PvtInterface();
 
-        virtual ~SinglePvtInterface();
+        virtual ~PvtInterface();
 
         /// \param[in]  num_phases   The number of active phases.
         /// \param[in]  phase_pos    Array of BlackpoilPhases::MaxNumPhases
@@ -44,13 +44,19 @@ namespace Opm
 
         /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
         /// or pressure (p) and gas resolution factor (r).
-        /// For all the virtual methods, the following apply: p, r and z
-        /// are expected to be of size n, size n and n*num_phases, respectively.
-        /// Output arrays shall be of size n, and must be valid before
-        /// calling the method.
+        /// For all the virtual methods, the following apply:
+        /// - pvtRegionIdx is an array of size n and represents the
+        ///   index of the PVT table which should be used to calculate
+        ///   the output. NULL can also be passed and is interpreted
+        ///   such that the first table should be used for the output
+        /// - p, r and z are expected to be of size n, size n and
+        ///   n*num_phases, respectively.
+        /// - Output arrays shall be of size n, and must be valid before
+        ///   calling the method.
 
         /// Viscosity as a function of p and z.
         virtual void mu(const int n,
+                        const int* pvtRegionIdx,
                         const double* p,
                         const double* z,
                         double* output_mu) const = 0;
@@ -58,6 +64,7 @@ namespace Opm
         /// Viscosity as a function of p and r.
         /// The fluid is considered saturated if r >= rsSat(p).
         virtual void mu(const int n,
+                        const int* pvtRegionIdx,
                               const double* p,
                               const double* r,
                               double* output_mu,
@@ -67,6 +74,7 @@ namespace Opm
         /// Viscosity as a function of p and r.
         /// State condition determined by 'cond'.
         virtual void mu(const int n,
+                        const int* pvtRegionIdx,
                               const double* p,
                               const double* r,
                               const PhasePresence* cond,
@@ -76,12 +84,14 @@ namespace Opm
 
         /// Formation volume factor as a function of p and z.
         virtual void B(const int n,
+                       const int* pvtRegionIdx,
                        const double* p,
                        const double* z,
                        double* output_B) const = 0;
 
         /// Formation volume factor and p-derivative as functions of p and z.
         virtual void dBdp(const int n,
+                          const int* pvtRegionIdx,
                           const double* p,
                           const double* z,
                           double* output_B,
@@ -90,6 +100,7 @@ namespace Opm
         /// The inverse of the volume factor b = 1 / B as a function of p and r.
         /// The fluid is considered saturated if r >= rsSat(p).
         virtual void b(const int n,
+                       const int* pvtRegionIdx,
                           const double* p,
                           const double* r,
                           double* output_b,
@@ -99,6 +110,7 @@ namespace Opm
         /// The inverse of the volume factor b = 1 / B as a function of p and r.
         /// State condition determined by 'cond'.
         virtual void b(const int n,
+                       const int* pvtRegionIdx,
                           const double* p,
                           const double* r,
                           const PhasePresence* cond,
@@ -108,12 +120,14 @@ namespace Opm
 
         /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
         virtual void rsSat(const int n,
+                           const int* pvtRegionIdx,
                           const double* p,
                           double* output_rsSat,
                           double* output_drsSatdp) const = 0;
 
         /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
         virtual void rvSat(const int n,
+                           const int* pvtRegionIdx,
                           const double* p,
                           double* output_rvSat,
                           double* output_drvSatdp) const = 0;
@@ -121,12 +135,14 @@ namespace Opm
 
         /// Solution factor as a function of p and z.
         virtual void R(const int n,
+                       const int* pvtRegionIdx,
                        const double* p,
                        const double* z,
                        double* output_R) const = 0;
 
         /// Solution factor and p-derivative as functions of p and z.
         virtual void dRdp(const int n,
+                          const int* pvtRegionIdx,
                           const double* p,
                           const double* z,
                           double* output_R,
@@ -138,7 +154,25 @@ namespace Opm
         int phase_pos_[MaxNumPhases];
     };
 
+    /*!
+     * \brief Helper function to create an array containing the (C-Style)
+     *        PVT table index for each compressed cell from an Eclipse deck.
+     *
+     * This function assumes that the degrees of freedom where PVT
+     * properties need to be calculated are grid cells. The main point
+     * of this function is to avoid code duplication because the
+     * Eclipse deck only contains Fortran-style PVT table indices
+     * which start at 1 instead of 0 and -- more relevantly -- it uses
+     * logically cartesian cell indices to specify the table index of
+     * a cell while the classes which use the PvtInterface
+     * implementations usually use compressed cells.
+     */
+    void extractPvtTableIndex(std::vector<int>& pvtTableIdx,
+                              Opm::DeckConstPtr deck,
+                              size_t numCompressed,
+                              const int* compressedToCartesianIdx);
+
 } // namespace Opm
 
-#endif // OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED
+#endif // OPM_PVTINTERFACE_HEADER_INCLUDED
 
diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp
new file mode 100644
index 00000000..7c59cc3a
--- /dev/null
+++ b/opm/core/props/pvt/PvtLiveGas.cpp
@@ -0,0 +1,526 @@
+//===========================================================================
+//
+// File: MiscibilityLiveGas.cpp
+//
+// Created: Wed Feb 10 09:21:53 2010
+//
+// Author: Bjørn Spjelkavik <bsp@sintef.no>
+//
+// Revision: $Id$
+//
+//===========================================================================
+/*
+  Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "config.h"
+#include <opm/core/props/pvt/PvtLiveGas.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/linearInterpolation.hpp>
+
+#include <algorithm>
+#include <iostream>
+
+namespace Opm
+{
+
+    using Opm::linearInterpolation;
+    using Opm::linearInterpolationDerivative;
+
+
+    //------------------------------------------------------------------------
+    // Member functions
+    //-------------------------------------------------------------------------
+    PvtLiveGas::PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword)
+    {
+        int numTables = Opm::PvtgTable::numTables(pvtgKeyword);
+        saturated_gas_table_.resize(numTables);
+        undersat_gas_tables_.resize(numTables);
+
+        for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) {
+            Opm::PvtgTable pvtgTable(pvtgKeyword, pvtTableIdx);
+
+            // GAS, PVTG
+            saturated_gas_table_[pvtTableIdx].resize(4);
+            saturated_gas_table_[pvtTableIdx][0] = pvtgTable.getOuterTable()->getPressureColumn();
+            saturated_gas_table_[pvtTableIdx][1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn();
+            saturated_gas_table_[pvtTableIdx][2] = pvtgTable.getOuterTable()->getGasViscosityColumn();
+            saturated_gas_table_[pvtTableIdx][3] = pvtgTable.getOuterTable()->getOilSolubilityColumn();
+
+            int sz = pvtgTable.getOuterTable()->numRows();
+            undersat_gas_tables_[pvtTableIdx].resize(sz);
+            for (int i=0; i<sz; ++i) {
+                const auto &undersatTable = *pvtgTable.getInnerTable(i);
+
+                undersat_gas_tables_[pvtTableIdx][i].resize(3);
+                undersat_gas_tables_[pvtTableIdx][i][0] = undersatTable.getOilSolubilityColumn();
+                undersat_gas_tables_[pvtTableIdx][i][1] = undersatTable.getGasFormationFactorColumn();
+                undersat_gas_tables_[pvtTableIdx][i][2] = undersatTable.getGasViscosityColumn();
+            }
+
+            // Bg -> 1/Bg
+            for (int i=0; i<sz; ++i) {
+                saturated_gas_table_[pvtTableIdx][1][i] = 1.0/saturated_gas_table_[pvtTableIdx][1][i];
+                for (size_t j=0; j<undersat_gas_tables_[pvtTableIdx][i][1].size(); ++j) {
+                    undersat_gas_tables_[pvtTableIdx][i][1][j] = 1.0/undersat_gas_tables_[pvtTableIdx][i][1][j];
+                }
+            }
+        }
+    }
+
+    // Destructor
+     PvtLiveGas::~PvtLiveGas()
+    {
+    }
+
+
+    void PvtLiveGas::mu(const int n,
+                        const int* pvtRegionIdx,
+                        const double* p,
+                        const double* z,
+                        double* output_mu) const
+    {
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            output_mu[i] = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 2, false);
+        }
+    }
+
+    /// Viscosity and its derivatives as a function of p and r.
+    void PvtLiveGas::mu(const int /*n*/,
+                        const int* /*pvtRegionIdx*/,
+                              const double* /*p*/,
+                              const double* /*r*/,
+                              double* /*output_mu*/,
+                              double* /*output_dmudp*/,
+                              double* /*output_dmudr*/) const
+    {
+        OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
+    }
+
+    /// Viscosity and its derivatives as a function of p and r.
+    void PvtLiveGas::mu(const int n,
+                        const int* pvtRegionIdx,
+                               const double* p,
+                               const double* r,
+                               const PhasePresence* cond,
+                               double* output_mu,
+                               double* output_dmudp,
+                               double* output_dmudr) const
+    {
+        for (int i = 0; i < n; ++i) {
+            const PhasePresence& cnd = cond[i];
+            int tableIdx = getTableIndex_(pvtRegionIdx, i);
+            output_mu[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 2, 0);
+            output_dmudp[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 2, 1);
+            output_dmudr[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 2, 2);
+        }
+
+    }
+
+
+    /// Formation volume factor as a function of p and z.
+    void PvtLiveGas::B(const int n,
+                       const int* pvtRegionIdx,
+                             const double* p,
+                             const double* z,
+                             double* output_B) const
+    {
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            output_B[i] = evalB(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i));
+        }
+
+    }
+
+
+    /// Formation volume factor and p-derivative as functions of p and z.
+    void PvtLiveGas::dBdp(const int n,
+                       const int* pvtRegionIdx,
+                                const double* p,
+                                const double* z,
+                                double* output_B,
+                                double* output_dBdp) const
+    {
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            evalBDeriv(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), output_B[i], output_dBdp[i]);
+        }
+    }
+
+    /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
+    void PvtLiveGas::b(const int /*n*/,
+                       const int* /*pvtRegionIdx*/,
+                             const double* /*p*/,
+                             const double* /*r*/,
+                             double* /*output_b*/,
+                             double* /*output_dbdp*/,
+                             double* /*output_dbdr*/) const
+
+    {
+        OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
+    }
+
+    /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
+    void PvtLiveGas::b(const int n,
+                       const int* pvtRegionIdx,
+                          const double* p,
+                          const double* r,
+                          const PhasePresence* cond,
+                          double* output_b,
+                          double* output_dbdp,
+                          double* output_dbdr) const
+
+    {
+        // #pragma omp parallel for
+                for (int i = 0; i < n; ++i) {
+                    const PhasePresence& cnd = cond[i];
+
+                    int tableIdx = getTableIndex_(pvtRegionIdx, i);
+                    output_b[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 0);
+                    output_dbdp[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 1);
+                    output_dbdr[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 2);
+
+                }
+    }
+
+    /// Gas resolution and its derivatives at bublepoint as a function of p.
+    void PvtLiveGas::rvSat(const int n,
+                           const int* pvtRegionIdx,
+                             const double* p,
+                             double* output_rvSat,
+                             double* output_drvSatdp) const
+    {
+        for (int i = 0; i < n; ++i) {
+            int pvtTableIdx = getTableIndex_(pvtRegionIdx, i);
+            output_rvSat[i] = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
+                    saturated_gas_table_[pvtTableIdx][3],p[i]);
+            output_drvSatdp[i] = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0],
+                    saturated_gas_table_[pvtTableIdx][3],p[i]);
+
+        }
+    }
+
+    void PvtLiveGas::rsSat(const int n,
+                           const int* pvtRegionIdx,
+                             const double* /*p*/,
+                             double* output_rsSat,
+                             double* output_drsSatdp) const
+    {
+        std::fill(output_rsSat, output_rsSat + n, 0.0);
+        std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
+    }
+
+    /// Solution factor as a function of p and z.
+    void PvtLiveGas::R(const int n,
+                       const int* pvtRegionIdx,
+                             const double* p,
+                             const double* z,
+                             double* output_R) const
+    {
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            output_R[i] = evalR(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i));
+        }
+
+    }
+
+
+    /// Solution factor and p-derivative as functions of p and z.
+    void PvtLiveGas::dRdp(const int n,
+                          const int* pvtRegionIdx,
+                                const double* p,
+                                const double* z,
+                                double* output_R,
+                                double* output_dRdp) const
+    {
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            evalRDeriv(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), output_R[i], output_dRdp[i]);
+        }
+    }
+
+
+    // ---- Private methods ----
+
+    double PvtLiveGas::evalB(const double press, const double* surfvol, int pvtTableIdx) const
+    {
+        if (surfvol[phase_pos_[Vapour]] == 0.0) {
+            // To handle no-gas case.
+            return 1.0;
+        }
+        return 1.0/miscible_gas(press, surfvol, pvtTableIdx, 1, false);
+    }
+
+    void PvtLiveGas::evalBDeriv(const double press, const double* surfvol, int pvtTableIdx,
+                                double& Bval, double& dBdpval) const
+    {
+        if (surfvol[phase_pos_[Vapour]] == 0.0) {
+            // To handle no-gas case.
+            Bval = 1.0;
+            dBdpval = 0.0;
+            return;
+        }
+        Bval = evalB(press, surfvol, pvtTableIdx);
+        dBdpval =  -Bval*Bval*miscible_gas(press, surfvol, pvtTableIdx, 1, true);
+    }
+
+    double PvtLiveGas::evalR(const double press, const double* surfvol, int pvtTableIdx) const
+    {
+        if (surfvol[phase_pos_[Liquid]] == 0.0) {
+            // To handle no-gas case.
+            return 0.0;
+        }
+        double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
+                                             saturated_gas_table_[pvtTableIdx][3], press);
+        double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
+        if (satR < maxR ) {
+            // Saturated case
+            return satR;
+        } else {
+            // Undersaturated case
+            return maxR;
+        }
+    }
+
+    void PvtLiveGas::evalRDeriv(const double press, const double* surfvol, int pvtTableIdx,
+                                double& Rval, double& dRdpval) const
+    {
+        if (surfvol[phase_pos_[Liquid]] == 0.0) {
+            // To handle no-gas case.
+            Rval = 0.0;
+            dRdpval = 0.0;
+            return;
+        }
+        double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
+                                             saturated_gas_table_[pvtTableIdx][3], press);
+        double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
+        if (satR < maxR ) {
+            // Saturated case
+            Rval = satR;
+            dRdpval = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0],
+                                            saturated_gas_table_[pvtTableIdx][3],
+                                            press);
+        } else {
+            // Undersaturated case
+            Rval = maxR;
+            dRdpval = 0.0;
+        }
+    }
+
+    double PvtLiveGas::miscible_gas(const double press,
+                                    const double* surfvol,
+                                    const int pvtTableIdx,
+                                    const int item,
+                                    const bool deriv) const
+    {
+        const std::vector<std::vector<double> > &saturatedGasTable =
+            saturated_gas_table_[pvtTableIdx];
+        const std::vector<std::vector<std::vector<double> > > &undersatGasTables =
+            undersat_gas_tables_[pvtTableIdx];
+
+        int section;
+        double Rval = linearInterpolation(saturatedGasTable[0],
+                                                saturatedGasTable[3], press,
+                                                section);
+        double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
+        if (deriv) {
+            if (Rval < maxR ) {  // Saturated case
+                return linearInterpolationDerivative(saturatedGasTable[0],
+                                                saturatedGasTable[item],
+                                                press);
+            } else {  // Undersaturated case
+                int is = section;
+                if (undersatGasTables[is][0].size() < 2) {
+                    double val = (saturatedGasTable[item][is+1]
+                                  - saturatedGasTable[item][is]) /
+                        (saturatedGasTable[0][is+1] -
+                         saturatedGasTable[0][is]);
+                    return val;
+                }
+                double val1 =
+                    linearInterpolation(undersatGasTables[is][0],
+                                              undersatGasTables[is][item],
+                                              maxR);
+                double val2 =
+                    linearInterpolation(undersatGasTables[is+1][0],
+                                              undersatGasTables[is+1][item],
+                                              maxR);
+                double val = (val2 - val1)/
+                    (saturatedGasTable[0][is+1] - saturatedGasTable[0][is]);
+                return val;
+            }
+        } else {
+            if (Rval < maxR ) {  // Saturated case
+                return linearInterpolation(saturatedGasTable[0],
+                                                 saturatedGasTable[item],
+                                                 press);
+            } else {  // Undersaturated case
+                int is = section;
+                // Extrapolate from first table section
+                if (is == 0 && press < saturatedGasTable[0][0]) {
+                    return linearInterpolation(undersatGasTables[0][0],
+                                                     undersatGasTables[0][item],
+                                                     maxR);
+                }
+
+                // Extrapolate from last table section
+                int ltp = saturatedGasTable[0].size() - 1;
+                if (is+1 == ltp && press > saturatedGasTable[0][ltp]) {
+                    return linearInterpolation(undersatGasTables[ltp][0],
+                                                     undersatGasTables[ltp][item],
+                                                     maxR);
+                }
+
+                // Interpolate between table sections
+                double w = (press - saturatedGasTable[0][is]) /
+                    (saturatedGasTable[0][is+1] -
+                     saturatedGasTable[0][is]);
+                if (undersatGasTables[is][0].size() < 2) {
+                    double val = saturatedGasTable[item][is] +
+                        w*(saturatedGasTable[item][is+1] -
+                           saturatedGasTable[item][is]);
+                    return val;
+                }
+                double val1 =
+                    linearInterpolation(undersatGasTables[is][0],
+                                              undersatGasTables[is][item],
+                                              maxR);
+                double val2 =
+                    linearInterpolation(undersatGasTables[is+1][0],
+                                              undersatGasTables[is+1][item],
+                                              maxR);
+                double val = val1 + w*(val2 - val1);
+                return val;
+            }
+        }
+    }
+
+    double PvtLiveGas::miscible_gas(const double press,
+                                    const double r,
+                                    const PhasePresence& cond,
+                                    const int pvtTableIdx,
+                                    const int item,
+                                    const int deriv) const
+    {
+        const std::vector<std::vector<double> > &saturatedGasTable =
+            saturated_gas_table_[pvtTableIdx];
+        const std::vector<std::vector<std::vector<double> > > &undersatGasTables =
+            undersat_gas_tables_[pvtTableIdx];
+
+        const bool isSat = cond.hasFreeOil();
+
+        // Derivative w.r.t p
+        if (deriv == 1) {
+            if (isSat) {  // Saturated case
+                return linearInterpolationDerivative(saturatedGasTable[0],
+                                                saturatedGasTable[item],
+                                                press);
+            } else {  // Undersaturated case
+                int is = tableIndex(saturatedGasTable[0], press);
+                if (undersatGasTables[is][0].size() < 2) {
+                    double val = (saturatedGasTable[item][is+1]
+                                  - saturatedGasTable[item][is]) /
+                        (saturatedGasTable[0][is+1] -
+                         saturatedGasTable[0][is]);
+                    return val;
+                }
+                double val1 =
+                    linearInterpolation(undersatGasTables[is][0],
+                                              undersatGasTables[is][item],
+                                              r);
+                double val2 =
+                    linearInterpolation(undersatGasTables[is+1][0],
+                                              undersatGasTables[is+1][item],
+                                              r);
+                double val = (val2 - val1)/
+                    (saturatedGasTable[0][is+1] - saturatedGasTable[0][is]);
+                return val;
+            }
+        } else if (deriv == 2){
+            if (isSat) {
+                return 0;
+            } else {
+                int is = tableIndex(saturatedGasTable[0], press);
+                double w = (press - saturatedGasTable[0][is]) /
+                    (saturatedGasTable[0][is+1] - saturatedGasTable[0][is]);
+                assert(undersatGasTables[is][0].size() >= 2);
+                assert(undersatGasTables[is+1][0].size() >= 2);
+                double val1 =
+                    linearInterpolationDerivative(undersatGasTables[is][0],
+                                              undersatGasTables[is][item],
+                                              r);
+                double val2 =
+                    linearInterpolationDerivative(undersatGasTables[is+1][0],
+                                              undersatGasTables[is+1][item],
+                                              r);
+
+                double val = val1 + w * (val2 - val1);
+                return val;
+
+            }
+        } else {
+            if (isSat) {  // Saturated case
+                return linearInterpolation(saturatedGasTable[0],
+                                                 saturatedGasTable[item],
+                                                 press);
+            } else {  // Undersaturated case
+                int is = tableIndex(saturatedGasTable[0], press);
+                // Extrapolate from first table section
+                if (is == 0 && press < saturatedGasTable[0][0]) {
+                    return linearInterpolation(undersatGasTables[0][0],
+                                                     undersatGasTables[0][item],
+                                                     r);
+                }
+
+                // Extrapolate from last table section
+                //int ltp = saturatedGasTable[0].size() - 1;
+                //if (is+1 == ltp && press > saturatedGasTable[0][ltp]) {
+                //    return linearInterpolation(undersatGasTables[ltp][0],
+                //                                    undersatGasTables[ltp][item],
+                //                                    r);
+                //}
+
+                // Interpolate between table sections
+                double w = (press - saturatedGasTable[0][is]) /
+                    (saturatedGasTable[0][is+1] -
+                     saturatedGasTable[0][is]);
+                if (undersatGasTables[is][0].size() < 2) {
+                    double val = saturatedGasTable[item][is] +
+                        w*(saturatedGasTable[item][is+1] -
+                           saturatedGasTable[item][is]);
+                    return val;
+                }
+                double val1 =
+                    linearInterpolation(undersatGasTables[is][0],
+                                              undersatGasTables[is][item],
+                                              r);
+                double val2 =
+                    linearInterpolation(undersatGasTables[is+1][0],
+                                              undersatGasTables[is+1][item],
+                                              r);
+                double val = val1 + w*(val2 - val1);
+                return val;
+            }
+        }
+    }
+
+
+
+
+} // namespace Opm
diff --git a/opm/core/props/pvt/SinglePvtLiveGas.hpp b/opm/core/props/pvt/PvtLiveGas.hpp
similarity index 74%
rename from opm/core/props/pvt/SinglePvtLiveGas.hpp
rename to opm/core/props/pvt/PvtLiveGas.hpp
index b61bd3e4..86ed258e 100644
--- a/opm/core/props/pvt/SinglePvtLiveGas.hpp
+++ b/opm/core/props/pvt/PvtLiveGas.hpp
@@ -17,10 +17,10 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#ifndef OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
-#define OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
+#ifndef OPM_PVTLIVEGAS_HEADER_INCLUDED
+#define OPM_PVTLIVEGAS_HEADER_INCLUDED
 
-#include <opm/core/props/pvt/SinglePvtInterface.hpp>
+#include <opm/core/props/pvt/PvtInterface.hpp>
 
 #include <opm/parser/eclipse/Utility/PvtgTable.hpp>
 
@@ -35,14 +35,15 @@ namespace Opm
     /// are expected to be of size n, size n and n*num_phases, respectively.
     /// Output arrays shall be of size n, and must be valid before
     /// calling the method.
-    class SinglePvtLiveGas : public SinglePvtInterface
+    class PvtLiveGas : public PvtInterface
     {
     public:
-        SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable);
-        virtual ~SinglePvtLiveGas();
+        PvtLiveGas(Opm::DeckKeywordConstPtr pvtgKeyword);
+        virtual ~PvtLiveGas();
 
         /// Viscosity as a function of p and z.
         virtual void mu(const int n,
+                        const int* pvtRegionIdx,
                         const double* p,
                         const double* z,
                         double* output_mu) const;
@@ -50,6 +51,7 @@ namespace Opm
         /// Viscosity and its derivatives as a function of p and r.
         /// The fluid is considered saturated if r >= rsSat(p).
         virtual void mu(const int n,
+                        const int* pvtRegionIdx,
                         const double* p,
                         const double* r,
                         double* output_mu,
@@ -59,6 +61,7 @@ namespace Opm
         /// Viscosity as a function of p and r.
         /// State condition determined by 'cond'.
         virtual void mu(const int n,
+                        const int* pvtRegionIdx,
                         const double* p,
                         const double* r,
                         const PhasePresence* cond,
@@ -68,12 +71,14 @@ namespace Opm
 
         /// Formation volume factor as a function of p and z.
         virtual void B(const int n,
+                       const int* pvtRegionIdx,
                        const double* p,
                        const double* z,
                        double* output_B) const;
 
         /// Formation volume factor and p-derivative as functions of p and z.
         virtual void dBdp(const int n,
+                          const int* pvtRegionIdx,
                           const double* p,
                           const double* z,
                           double* output_B,
@@ -82,6 +87,7 @@ namespace Opm
         /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
         /// The fluid is considered saturated if r >= rsSat(p).
         virtual void b(const int n,
+                       const int* pvtRegionIdx,
                        const double* p,
                        const double* r,
                        double* output_b,
@@ -91,6 +97,7 @@ namespace Opm
         /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
         /// State condition determined by 'cond'.
         virtual void b(const int n,
+                       const int* pvtRegionIdx,
                        const double* p,
                        const double* r,
                        const PhasePresence* cond,
@@ -102,52 +109,65 @@ namespace Opm
 
         /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
         virtual void rsSat(const int n,
+                           const int* pvtRegionIdx,
                           const double* p,
                           double* output_rsSat,
                           double* output_drsSatdp) const;
 
         /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
         virtual void rvSat(const int n,
+                           const int* pvtRegionIdx,
                           const double* p,
                           double* output_rvSat,
                           double* output_drvSatdp) const;
 
         /// Solution factor as a function of p and z.
         virtual void R(const int n,
+                       const int* pvtRegionIdx,
                        const double* p,
                        const double* z,
                        double* output_R) const;
 
         /// Solution factor and p-derivative as functions of p and z.
         virtual void dRdp(const int n,
+                          const int* pvtRegionIdx,
                           const double* p,
                           const double* z,
                           double* output_R,
                           double* output_dRdp) const;
 
     protected:
-        double evalB(double press, const double* surfvol) const;
-        void evalBDeriv(double press, const double* surfvol, double& B, double& dBdp) const;
-        double evalR(double press, const double* surfvol) const;
-        void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const;
+        int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
+        {
+            if (!pvtTableIdx)
+                return 0;
+            return pvtTableIdx[cellIdx];
+        }
+
+        double evalB(double press, const double* surfvol, int pvtTableIdx) const;
+        void evalBDeriv(double press, const double* surfvol, int pvtTableIdx, double& B, double& dBdp) const;
+        double evalR(double press, const double* surfvol, int pvtTableIdx) const;
+        void evalRDeriv(double press, const double* surfvol, int pvtTableIdx, double& R, double& dRdp) const;
 
         // item:  1=>B  2=>mu;
         double miscible_gas(const double press,
                             const double* surfvol,
+                            const int pvtTableIdx,
                             const int item,
                             const bool deriv = false) const;
         double miscible_gas(const double press,
                             const double r,
                             const PhasePresence& cond,
+                            const int pvtTableIdx,
                             const int item,
                             const int deriv = 0) const;
-        // PVT properties of wet gas (with vaporised oil)
-        std::vector<std::vector<double> > saturated_gas_table_;
-        std::vector<std::vector<std::vector<double> > > undersat_gas_tables_;
-
+        // PVT properties of wet gas (with vaporised oil). We need to
+        // store one table per PVT region.
+        std::vector< std::vector<std::vector<double> > > saturated_gas_table_;
+        std::vector< std::vector<std::vector<std::vector<double> > > > undersat_gas_tables_;
     };
 
 }
 
-#endif // OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
+#endif // OPM_PVTLIVEGAS_HEADER_INCLUDED
 
diff --git a/opm/core/props/pvt/PvtLiveOil.cpp b/opm/core/props/pvt/PvtLiveOil.cpp
new file mode 100644
index 00000000..578e6369
--- /dev/null
+++ b/opm/core/props/pvt/PvtLiveOil.cpp
@@ -0,0 +1,591 @@
+/*
+  Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "config.h"
+#include <opm/core/props/pvt/PvtLiveOil.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/linearInterpolation.hpp>
+
+#include <algorithm>
+
+namespace Opm
+{
+
+    using Opm::linearInterpolation;
+    using Opm::linearInterpolationDerivative;
+    using Opm::tableIndex;
+
+    //------------------------------------------------------------------------
+    // Member functions
+    //-------------------------------------------------------------------------
+    PvtLiveOil::PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword)
+    {
+        int numTables = Opm::PvtoTable::numTables(pvtoKeyword);
+        saturated_oil_table_.resize(numTables);
+        undersat_oil_tables_.resize(numTables);
+
+        for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) {
+            Opm::PvtoTable pvtoTable(pvtoKeyword, pvtTableIdx);
+
+            const auto saturatedPvto = pvtoTable.getOuterTable();
+
+            // OIL, PVTO
+            saturated_oil_table_[pvtTableIdx].resize(4);
+            const int sz = saturatedPvto->numRows();
+            for (int k=0; k<4; ++k) {
+                saturated_oil_table_[pvtTableIdx][k].resize(sz);
+            }
+            for (int i=0; i<sz; ++i) {
+                saturated_oil_table_[pvtTableIdx][0][i] = saturatedPvto->getPressureColumn()[i]; // p
+                saturated_oil_table_[pvtTableIdx][1][i] = 1.0/saturatedPvto->getOilFormationFactorColumn()[i]; // 1/Bo
+                saturated_oil_table_[pvtTableIdx][2][i] = saturatedPvto->getOilViscosityColumn()[i]; // mu_o
+                saturated_oil_table_[pvtTableIdx][3][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs
+            }
+
+            undersat_oil_tables_[pvtTableIdx].resize(sz);
+            for (int i=0; i<sz; ++i) {
+                const auto undersaturatedPvto = pvtoTable.getInnerTable(i);
+
+                undersat_oil_tables_[pvtTableIdx][i].resize(3);
+                int tsize = undersaturatedPvto->numRows();
+                undersat_oil_tables_[pvtTableIdx][i][0].resize(tsize);
+                undersat_oil_tables_[pvtTableIdx][i][1].resize(tsize);
+                undersat_oil_tables_[pvtTableIdx][i][2].resize(tsize);
+                for (int j=0; j<tsize; ++j) {
+                    undersat_oil_tables_[pvtTableIdx][i][0][j] = undersaturatedPvto->getPressureColumn()[j];  // p
+                    undersat_oil_tables_[pvtTableIdx][i][1][j] = 1.0/undersaturatedPvto->getOilFormationFactorColumn()[j];  // 1/Bo
+                    undersat_oil_tables_[pvtTableIdx][i][2][j] = undersaturatedPvto->getOilViscosityColumn()[j];  // mu_o
+                }
+            }
+
+            // Complete undersaturated tables by extrapolating from existing data
+            // as is done in Eclipse and Mrst
+            int iNext = -1;
+            for (int i=0; i<sz; ++i) {
+                // Skip records already containing undersaturated data
+                if (undersat_oil_tables_[pvtTableIdx][i][0].size() > 1) {
+                    continue;
+                }
+                // Look ahead for next record containing undersaturated data
+                if (iNext < i) {
+                    iNext = i+1;
+                    while (iNext<sz && undersat_oil_tables_[pvtTableIdx][iNext][0].size() < 2) {
+                        ++iNext;
+                    }
+                    if (iNext == sz) OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
+                }
+                // Add undersaturated data to current record while maintaining compressibility and viscosibility
+                typedef std::vector<std::vector<std::vector<double> > >::size_type sz_t;
+                for (sz_t j=1; j<undersat_oil_tables_[pvtTableIdx][iNext][0].size(); ++j) {
+                    double diffPressure = undersat_oil_tables_[pvtTableIdx][iNext][0][j]-undersat_oil_tables_[pvtTableIdx][iNext][0][j-1];
+                    double pressure = undersat_oil_tables_[pvtTableIdx][i][0].back()+diffPressure;
+                    undersat_oil_tables_[pvtTableIdx][i][0].push_back(pressure);
+                    double compr = (1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j]-1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j-1])
+                        / (0.5*(1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j]+1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j-1]));
+                    double B = (1.0/undersat_oil_tables_[pvtTableIdx][i][1].back())*(1.0+0.5*compr)/(1.0-0.5*compr);
+                    undersat_oil_tables_[pvtTableIdx][i][1].push_back(1.0/B);
+                    double visc = (undersat_oil_tables_[pvtTableIdx][iNext][2][j]-undersat_oil_tables_[pvtTableIdx][iNext][2][j-1])
+                        / (0.5*(undersat_oil_tables_[pvtTableIdx][iNext][2][j]+undersat_oil_tables_[pvtTableIdx][iNext][2][j-1]));
+                    double mu = (undersat_oil_tables_[pvtTableIdx][i][2].back())*(1.0+0.5*visc)/(1.0-0.5*visc);
+                    undersat_oil_tables_[pvtTableIdx][i][2].push_back(mu);
+                }
+            }
+        }
+    }
+
+    /// Destructor.
+    PvtLiveOil::~PvtLiveOil()
+    {
+    }
+
+
+    /// Viscosity as a function of p and z.
+    void PvtLiveOil::mu(const int n,
+                        const int* pvtTableIdx,
+                              const double* p,
+                              const double* z,
+                              double* output_mu) const
+    {
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            int tableIdx = getTableIndex_(pvtTableIdx, i);
+
+            output_mu[i] = miscible_oil(p[i], z + num_phases_*i, tableIdx, 2, false);
+        }
+    }
+
+    /// Viscosity and its derivatives as a function of p and r.
+    void PvtLiveOil::mu(const int n,
+                        const int* pvtTableIdx,
+                               const double* p,
+                               const double* r,
+                               double* output_mu,
+                               double* output_dmudp,
+                               double* output_dmudr) const
+    {
+        // #pragma omp parallel for
+                for (int i = 0; i < n; ++i) {
+                    int tableIdx = getTableIndex_(pvtTableIdx, i);
+
+                    output_mu[i] = miscible_oil(p[i], r[i], tableIdx, 2, 0);
+                    output_dmudp[i] = miscible_oil(p[i], r[i], tableIdx, 2, 1);
+                    output_dmudr[i] = miscible_oil(p[i], r[i], tableIdx, 2, 2);
+
+                }
+    }
+
+    /// Viscosity and its derivatives as a function of p and r.
+    void PvtLiveOil::mu(const int n,
+                        const int* pvtTableIdx,
+                               const double* p,
+                               const double* r,
+                               const PhasePresence* cond,
+                               double* output_mu,
+                               double* output_dmudp,
+                               double* output_dmudr) const
+    {
+        // #pragma omp parallel for
+                for (int i = 0; i < n; ++i) {
+                    int tableIdx = getTableIndex_(pvtTableIdx, i);
+                    const PhasePresence& cnd = cond[i];
+
+                    output_mu[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 2, 0);
+                    output_dmudp[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 2, 1);
+                    output_dmudr[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 2, 2);
+
+                }
+    }
+
+
+    /// Formation volume factor as a function of p and z.
+    void PvtLiveOil::B(const int n,
+                       const int* pvtTableIdx,
+                             const double* p,
+                             const double* z,
+                             double* output_B) const
+    {
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            int tableIdx = getTableIndex_(pvtTableIdx, i);
+
+            output_B[i] = evalB(tableIdx, p[i], z + num_phases_*i);
+        }
+
+    }
+
+
+    /// Formation volume factor and p-derivative as functions of p and z.
+    void PvtLiveOil::dBdp(const int n,
+                        const int* pvtTableIdx,
+                                const double* p,
+                                const double* z,
+                                double* output_B,
+                                double* output_dBdp) const
+    {
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            int tableIdx = getTableIndex_(pvtTableIdx, i);
+
+            evalBDeriv(tableIdx, p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
+        }
+    }
+
+    void PvtLiveOil::b(const int n,
+                        const int* pvtTableIdx,
+                          const double* p,
+                          const double* r,
+                          double* output_b,
+                          double* output_dbdp,
+                          double* output_dbdr) const
+
+    {
+        // #pragma omp parallel for
+                for (int i = 0; i < n; ++i) {
+                    int tableIdx = getTableIndex_(pvtTableIdx, i);
+
+                    output_b[i] = miscible_oil(tableIdx, p[i], r[i], 1, 0);
+                    output_dbdp[i] = miscible_oil(tableIdx, p[i], r[i], 1, 1);
+                    output_dbdr[i] = miscible_oil(tableIdx, p[i], r[i], 1, 2);
+
+                }
+    }
+
+    void PvtLiveOil::b(const int n,
+                        const int* pvtTableIdx,
+                          const double* p,
+                          const double* r,
+                          const PhasePresence* cond,
+                          double* output_b,
+                          double* output_dbdp,
+                          double* output_dbdr) const
+
+    {
+        // #pragma omp parallel for
+                for (int i = 0; i < n; ++i) {
+                    const PhasePresence& cnd = cond[i];
+                    int tableIdx = getTableIndex_(pvtTableIdx, i);
+
+                    output_b[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 0);
+                    output_dbdp[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 1);
+                    output_dbdr[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 2);
+
+                }
+    }
+
+    void PvtLiveOil::rsSat(const int n,
+                           const int* pvtTableIdx,
+                             const double* p,
+                             double* output_rsSat,
+                             double* output_drsSatdp) const
+    {
+
+        for (int i = 0; i < n; ++i) {
+            int tableIdx = getTableIndex_(pvtTableIdx, i);
+            output_rsSat[i] = linearInterpolation(saturated_oil_table_[tableIdx][0],
+                    saturated_oil_table_[tableIdx][3],p[i]);
+            output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[tableIdx][0],
+                    saturated_oil_table_[tableIdx][3],p[i]);
+
+        }
+    }
+
+    void PvtLiveOil::rvSat(const int n,
+                           const int* pvtTableIdx,
+                             const double* /*p*/,
+                             double* output_rvSat,
+                             double* output_drvSatdp) const
+    {
+        std::fill(output_rvSat, output_rvSat + n, 0.0);
+        std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
+    }
+
+    /// Solution factor as a function of p and z.
+    void PvtLiveOil::R(const int n,
+                       const int* pvtTableIdx,
+                             const double* p,
+                             const double* z,
+                             double* output_R) const
+    {
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            int tableIdx = getTableIndex_(pvtTableIdx, i);
+            output_R[i] = evalR(tableIdx, p[i], z + num_phases_*i);
+        }
+
+    }
+
+
+    /// Solution factor and p-derivative as functions of p and z.
+    void PvtLiveOil::dRdp(const int n,
+                          const int* pvtTableIdx,
+                                const double* p,
+                                const double* z,
+                                double* output_R,
+                                double* output_dRdp) const
+    {
+// #pragma omp parallel for
+        for (int i = 0; i < n; ++i) {
+            int tableIdx = getTableIndex_(pvtTableIdx, i);
+            evalRDeriv(tableIdx, p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
+        }
+    }
+
+
+
+
+    // ---- Private methods ----
+
+    double PvtLiveOil::evalB(size_t pvtTableIdx,
+                                   double press, const double* surfvol) const
+    {
+        // if (surfvol[phase_pos_[Liquid]] == 0.0) return 1.0; // To handle no-oil case.
+        return 1.0/miscible_oil(press, surfvol, pvtTableIdx, 1, false);
+    }
+
+
+    void PvtLiveOil::evalBDeriv(size_t pvtTableIdx,
+                                      const double press, const double* surfvol,
+                                      double& Bval, double& dBdpval) const
+    {
+        Bval = evalB(pvtTableIdx, press, surfvol);
+        dBdpval = -Bval*Bval*miscible_oil(press, surfvol, pvtTableIdx, 1, true);
+    }
+
+    double PvtLiveOil::evalR(size_t pvtTableIdx,
+                                   double press, const double* surfvol) const
+    {
+        if (surfvol[phase_pos_[Vapour]] == 0.0) {
+            return 0.0;
+        }
+        double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
+                                             saturated_oil_table_[pvtTableIdx][3], press);
+        double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
+        if (Rval < maxR ) {  // Saturated case
+            return Rval;
+        } else {
+            return maxR;  // Undersaturated case
+        }
+    }
+
+    void PvtLiveOil::evalRDeriv(size_t pvtTableIdx,
+                                      const double press, const double* surfvol,
+                                      double& Rval, double& dRdpval) const
+    {
+        if (surfvol[phase_pos_[Vapour]] == 0.0) {
+            Rval = 0.0;
+            dRdpval = 0.0;
+            return;
+        }
+        Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
+                                   saturated_oil_table_[pvtTableIdx][3], press);
+        double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
+        if (Rval < maxR ) {
+            // Saturated case
+            dRdpval = linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
+                                                    saturated_oil_table_[pvtTableIdx][3],
+                                                    press);
+        } else {
+            // Undersaturated case
+            Rval = maxR;
+            dRdpval = 0.0;
+        }
+    }
+
+
+    double PvtLiveOil::miscible_oil(const double press,
+                                    const double* surfvol,
+                                    const int pvtTableIdx,
+                                    const int item,
+                                    const bool deriv) const
+    {
+        int section;
+        double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
+                                          saturated_oil_table_[pvtTableIdx][3],
+                                          press, section);
+        double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
+        if (deriv) {
+            if (Rval < maxR ) {  // Saturated case
+                return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
+                                                     saturated_oil_table_[pvtTableIdx][item],
+                                                     press);
+            } else {  // Undersaturated case
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], maxR);
+                double w = (maxR - saturated_oil_table_[pvtTableIdx][3][is]) /
+                    (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
+                assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
+                assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
+                double val1 =
+                    linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0],
+                                                  undersat_oil_tables_[pvtTableIdx][is][item],
+                                                  press);
+                double val2 =
+                    linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0],
+                                                  undersat_oil_tables_[pvtTableIdx][is+1][item],
+                                                  press);
+                double val = val1 + w*(val2 - val1);
+                return val;
+            }
+        } else {
+            if (Rval < maxR ) {  // Saturated case
+                return linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
+                                           saturated_oil_table_[pvtTableIdx][item],
+                                           press);
+            } else {  // Undersaturated case
+                // Interpolate between table sections
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], maxR);
+                double w = (maxR - saturated_oil_table_[pvtTableIdx][3][is]) /
+                    (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
+                assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
+                assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
+                double val1 =
+                    linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
+                                        undersat_oil_tables_[pvtTableIdx][is][item],
+                                        press);
+                double val2 =
+                    linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
+                                        undersat_oil_tables_[pvtTableIdx][is+1][item],
+                                        press);
+                double val = val1 + w*(val2 - val1);
+                return val;
+            }
+        }
+    }
+
+    double PvtLiveOil::miscible_oil(const double press,
+                                    const double r,
+                                    const int pvtTableIdx,
+                                    const int item,
+                                    const int deriv) const
+    {
+        int section;
+        double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
+                                          saturated_oil_table_[pvtTableIdx][3],
+                                          press, section);
+        // derivative with respect to frist component (pressure)
+        if (deriv == 1) {
+            if (Rval <= r ) {  // Saturated case
+                return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
+                                                     saturated_oil_table_[pvtTableIdx][item],
+                                                     press);
+            } else {  // Undersaturated case
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
+                double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) /
+                    (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
+                assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
+                assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
+                double val1 =
+                    linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0],
+                                             undersat_oil_tables_[pvtTableIdx][is][item],
+                                             press);
+                double val2 =
+                    linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0],
+                                             undersat_oil_tables_[pvtTableIdx][is+1][item],
+                                             press);
+                double val = val1 + w*(val2 - val1);
+                return val;
+            }
+            // derivative with respect to second component (r)
+        } else if (deriv == 2)  {
+            if (Rval <= r ) {  // Saturated case
+                return 0;
+            } else {  // Undersaturated case
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
+                assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
+                assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
+                double val1 =
+                    linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
+                                              undersat_oil_tables_[pvtTableIdx][is][item],
+                                              press);
+                double val2 =
+                    linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
+                                              undersat_oil_tables_[pvtTableIdx][is+1][item],
+                                              press);
+
+                double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][3][is+1]-saturated_oil_table_[pvtTableIdx][3][is]);
+                return val;
+            }
+
+
+            }    else {
+            if (Rval <= r ) {  // Saturated case
+                return linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
+                                           saturated_oil_table_[pvtTableIdx][item],
+                                           press);
+            } else {  // Undersaturated case
+                // Interpolate between table sections
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
+                double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) /
+                    (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
+                assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
+                assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
+                double val1 =
+                    linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
+                                              undersat_oil_tables_[pvtTableIdx][is][item],
+                                              press);
+                double val2 =
+                    linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
+                                              undersat_oil_tables_[pvtTableIdx][is+1][item],
+                                              press);
+                double val = val1 + w*(val2 - val1);
+                return val;
+            }
+        }
+    }
+
+    double PvtLiveOil::miscible_oil(const double press,
+                                    const double r,
+                                    const PhasePresence& cond,
+                                    const int pvtTableIdx,
+                                    const int item,
+                                    const int deriv) const
+    {
+        const bool isSat = cond.hasFreeGas();
+
+        // derivative with respect to frist component (pressure)
+        if (deriv == 1) {
+            if (isSat) {  // Saturated case
+                return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
+                                                     saturated_oil_table_[pvtTableIdx][item],
+                                                     press);
+            } else {  // Undersaturated case
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
+                double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) /
+                    (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
+                assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
+                assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
+                double val1 =
+                    linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0],
+                                             undersat_oil_tables_[pvtTableIdx][is][item],
+                                             press);
+                double val2 =
+                    linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0],
+                                             undersat_oil_tables_[pvtTableIdx][is+1][item],
+                                             press);
+                double val = val1 + w*(val2 - val1);
+                return val;
+            }
+            // derivative with respect to second component (r)
+        } else if (deriv == 2)  {
+            if (isSat) {  // Saturated case
+                return 0;
+            } else {  // Undersaturated case
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
+                assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
+                assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
+                double val1 =
+                    linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
+                                              undersat_oil_tables_[pvtTableIdx][is][item],
+                                              press);
+                double val2 =
+                    linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
+                                              undersat_oil_tables_[pvtTableIdx][is+1][item],
+                                              press);
+
+                double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][3][is+1]-saturated_oil_table_[pvtTableIdx][3][is]);
+                return val;
+            }
+
+
+            }    else {
+            if (isSat) {  // Saturated case
+                return linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
+                                           saturated_oil_table_[pvtTableIdx][item],
+                                           press);
+            } else {  // Undersaturated case
+                // Interpolate between table sections
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
+                double w = (r - saturated_oil_table_[pvtTableIdx][3][is]) /
+                    (saturated_oil_table_[pvtTableIdx][3][is+1] - saturated_oil_table_[pvtTableIdx][3][is]);
+                assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
+                assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
+                double val1 =
+                    linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
+                                              undersat_oil_tables_[pvtTableIdx][is][item],
+                                              press);
+                double val2 =
+                    linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
+                                              undersat_oil_tables_[pvtTableIdx][is+1][item],
+                                              press);
+                double val = val1 + w*(val2 - val1);
+                return val;
+            }
+        }
+    }
+
+} // namespace Opm
diff --git a/opm/core/props/pvt/SinglePvtLiveOil.hpp b/opm/core/props/pvt/PvtLiveOil.hpp
similarity index 74%
rename from opm/core/props/pvt/SinglePvtLiveOil.hpp
rename to opm/core/props/pvt/PvtLiveOil.hpp
index 8f5328f1..b6f61644 100644
--- a/opm/core/props/pvt/SinglePvtLiveOil.hpp
+++ b/opm/core/props/pvt/PvtLiveOil.hpp
@@ -17,11 +17,12 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#ifndef OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
-#define OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
+#ifndef OPM_PVTLIVEOIL_HEADER_INCLUDED
+#define OPM_PVTLIVEOIL_HEADER_INCLUDED
 
-#include <opm/core/props/pvt/SinglePvtInterface.hpp>
+#include <opm/core/props/pvt/PvtInterface.hpp>
 
+#include <opm/parser/eclipse/Deck/Deck.hpp>
 #include <opm/parser/eclipse/Utility/PvtoTable.hpp>
 
 #include <vector>
@@ -35,14 +36,15 @@ namespace Opm
     /// are expected to be of size n, size n and n*num_phases, respectively.
     /// Output arrays shall be of size n, and must be valid before
     /// calling the method.
-    class SinglePvtLiveOil : public SinglePvtInterface
+    class PvtLiveOil : public PvtInterface
     {
     public:
-        SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable);
-        virtual ~SinglePvtLiveOil();
+        PvtLiveOil(Opm::DeckKeywordConstPtr pvtoKeyword);
+        virtual ~PvtLiveOil();
 
         /// Viscosity as a function of p and z.
         virtual void mu(const int n,
+                        const int* pvtTableIdx,
                         const double* p,
                         const double* z,
                         double* output_mu) const;
@@ -50,6 +52,7 @@ namespace Opm
         /// Viscosity and its derivatives as a function of p and r.
         /// The fluid is considered saturated if r >= rsSat(p).
         virtual void mu(const int n,
+                        const int* pvtTableIdx,
                         const double* p,
                         const double* r,
                         double* output_mu,
@@ -59,6 +62,7 @@ namespace Opm
         /// Viscosity as a function of p and r.
         /// State condition determined by 'cond'.
         virtual void mu(const int n,
+                        const int* pvtTableIdx,
                         const double* p,
                         const double* r,
                         const PhasePresence* cond,
@@ -68,12 +72,14 @@ namespace Opm
 
         /// Formation volume factor as a function of p and z.
         virtual void B(const int n,
+                       const int* pvtTableIdx,
                        const double* p,
                        const double* z,
                        double* output_B) const;
 
         /// Formation volume factor and p-derivative as functions of p and z.
         virtual void dBdp(const int n,
+                          const int* pvtTableIdx,
                           const double* p,
                           const double* z,
                           double* output_B,
@@ -82,6 +88,7 @@ namespace Opm
         /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
         /// The fluid is considered saturated if r >= rsSat(p).
         virtual void b(const int n,
+                       const int* pvtTableIdx,
                        const double* p,
                        const double* r,
                        double* output_b,
@@ -91,6 +98,7 @@ namespace Opm
         /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
         /// State condition determined by 'cond'.
         virtual void b(const int n,
+                       const int* pvtTableIdx,
                        const double* p,
                        const double* r,
                        const PhasePresence* cond,
@@ -100,58 +108,73 @@ namespace Opm
 
         /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
         virtual void rsSat(const int n,
+                           const int* pvtTableIdx,
                           const double* p,
                           double* output_rsSat,
                           double* output_drsSatdp) const;
 
         /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
         virtual void rvSat(const int n,
+                           const int* pvtTableIdx,
                           const double* p,
                           double* output_rvSat,
                           double* output_drvSatdp) const;
 
         /// Solution factor as a function of p and z.
         virtual void R(const int n,
+                       const int* pvtTableIdx,
                        const double* p,
                        const double* z,
                        double* output_R) const;
 
         /// Solution factor and p-derivative as functions of p and z.
         virtual void dRdp(const int n,
+                          const int* pvtTableIdx,
                           const double* p,
                           const double* z,
                           double* output_R,
                           double* output_dRdp) const;
 
     private:
-        double evalB(double press, const double* surfvol) const;
-        void evalBDeriv(double press, const double* surfvol, double& B, double& dBdp) const;
-        double evalR(double press, const double* surfvol) const;
-        void evalRDeriv(double press, const double* surfvol, double& R, double& dRdp) const;
+        int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
+        {
+            if (!pvtTableIdx)
+                return 0;
+            return pvtTableIdx[cellIdx];
+        }
+
+        double evalB(size_t pvtTableIdx, double press, const double* surfvol) const;
+        void evalBDeriv(size_t pvtTableIdx, double press, const double* surfvol, double& B, double& dBdp) const;
+        double evalR(size_t pvtTableIdx, double press, const double* surfvol) const;
+        void evalRDeriv(size_t pvtTableIdx, double press, const double* surfvol, double& R, double& dRdp) const;
 
         // item:  1=>1/B  2=>mu;
         double miscible_oil(const double press,
                             const double* surfvol,
+                            const int pvtTableIdx,
                             const int item,
                             const bool deriv = false) const;
 
         double miscible_oil(const double press,
                             const double r,
+                            const int pvtTableIdx,
                             const int item,
                             const int deriv = 0) const;
 
         double miscible_oil(const double press,
                             const double r,
                             const PhasePresence& cond,
+                            const int pvtTableIdx,
                             const int item,
                             const int deriv = 0) const;
 
-        // PVT properties of live oil (with dissolved gas)
-        std::vector<std::vector<double> > saturated_oil_table_;
-        std::vector<std::vector<std::vector<double> > > undersat_oil_tables_;
+        // PVT properties of live oil (with dissolved gas). We need to
+        // store one table per PVT region.
+        std::vector<std::vector<std::vector<double> > > saturated_oil_table_;
+        std::vector<std::vector<std::vector<std::vector<double> > > > undersat_oil_tables_;
     };
 
 }
 
-#endif // OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
+#endif // OPM_PVTLIVEOIL_HEADER_INCLUDED
 
diff --git a/opm/core/props/pvt/PvtPropertiesBasic.hpp b/opm/core/props/pvt/PvtPropertiesBasic.hpp
index e9198ff9..46add622 100644
--- a/opm/core/props/pvt/PvtPropertiesBasic.hpp
+++ b/opm/core/props/pvt/PvtPropertiesBasic.hpp
@@ -95,6 +95,8 @@ namespace Opm
                   double* output_dRdp) const;
 
     private:
+        // The PVT properties. We need to store one value per PVT
+        // region.
         std::vector<double> density_;
         std::vector<double> viscosity_;
         std::vector<double> formation_volume_factor_;
diff --git a/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp b/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
index c4c764ac..27176351 100644
--- a/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
+++ b/opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
@@ -33,9 +33,9 @@ namespace Opm
     {
     }
 
-    void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck )
+    void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck)
     {
-        // If we need multiple regions, this class and the SinglePvt* classes must change.
+        // So far, this class only supports a single PVT region. TODO?
         int region_number = 0;
 
         PhaseUsage phase_usage = phaseUsageFromDeck(deck);
diff --git a/opm/core/props/pvt/SinglePvtConstCompr.hpp b/opm/core/props/pvt/SinglePvtConstCompr.hpp
deleted file mode 100644
index 180611ca..00000000
--- a/opm/core/props/pvt/SinglePvtConstCompr.hpp
+++ /dev/null
@@ -1,291 +0,0 @@
-/*
-  Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
-
-  This file is part of the Open Porous Media project (OPM).
-
-  OPM is free software: you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation, either version 3 of the License, or
-  (at your option) any later version.
-
-  OPM is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-
-  You should have received a copy of the GNU General Public License
-  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#ifndef OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
-#define OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
-
-#include <opm/core/props/pvt/SinglePvtInterface.hpp>
-#include <opm/core/utility/ErrorMacros.hpp>
-
-#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
-#include <opm/parser/eclipse/Utility/PvcdoTable.hpp>
-
-#include <vector>
-#include <algorithm>
-
-
-namespace Opm
-{
-
-    /// Class for constant compressible phases (PVTW or PVCDO).
-    /// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
-    /// or pressure (p) and gas resolution factor (r).
-    /// For all the virtual methods, the following apply: p, r and z
-    /// are expected to be of size n, size n and n*num_phases, respectively.
-    /// Output arrays shall be of size n, and must be valid before
-    /// calling the method.
-    class SinglePvtConstCompr : public SinglePvtInterface
-    {
-    public:
-        SinglePvtConstCompr(const Opm::PvtwTable &pvtwTable)
-        {
-            if (pvtwTable.numRows() != 1)
-                OPM_THROW(std::runtime_error,
-                          "The table specified by the PVTW keyword is required"
-                          "to exhibit exactly one row.");
-            ref_press_ = pvtwTable.getPressureColumn()[0];
-            ref_B_     = pvtwTable.getFormationFactorColumn()[0];
-            comp_      = pvtwTable.getCompressibilityColumn()[0];
-            viscosity_ = pvtwTable.getViscosityColumn()[0];
-            visc_comp_ = pvtwTable.getViscosibilityColumn()[0];
-        }
-
-        SinglePvtConstCompr(const Opm::PvcdoTable &pvcdoTable)
-        {
-            if (pvcdoTable.numRows() != 1)
-                OPM_THROW(std::runtime_error,
-                          "The table specified by the PVCDO keyword is required"
-                          "to exhibit exactly one row.");
-            ref_press_ = pvcdoTable.getPressureColumn()[0];
-            ref_B_     = pvcdoTable.getFormationFactorColumn()[0];
-            comp_      = pvcdoTable.getCompressibilityColumn()[0];
-            viscosity_ = pvcdoTable.getViscosityColumn()[0];
-            visc_comp_ = pvcdoTable.getViscosibilityColumn()[0];
-        }
-
-        SinglePvtConstCompr(double visc)
-            : ref_press_(0.0),
-              ref_B_(1.0),
-              comp_(0.0),
-              viscosity_(visc),
-              visc_comp_(0.0)
-        {
-        }
-
-        virtual ~SinglePvtConstCompr()
-        {
-        }
-
-        virtual void mu(const int n,
-                        const double* p,
-                        const double* /*z*/,
-                        double* output_mu) const
-        {
-            if (visc_comp_) {
-// #pragma omp parallel for
-                for (int i = 0; i < n; ++i) {
-                    // Computing a polynomial approximation to the exponential.
-                    double x = -visc_comp_*(p[i] - ref_press_);
-                    output_mu[i] = viscosity_/(1.0 + x + 0.5*x*x);
-                }
-            } else {
-                std::fill(output_mu, output_mu + n, viscosity_);
-            }
-        }
-
-        virtual void mu(const int n,
-                        const double* p,
-                        const double* /*r*/,
-                        double* output_mu,
-                        double* output_dmudp,
-                        double* output_dmudr) const
-        {
-            if (visc_comp_) {
-                // #pragma omp parallel for
-                for (int i = 0; i < n; ++i) {
-                    // Computing a polynomial approximation to the exponential.
-                    double x = -visc_comp_*(p[i] - ref_press_);
-                    double d = (1.0 + x + 0.5*x*x);
-                    output_mu[i] = viscosity_/d;
-                    output_dmudp[i] = (viscosity_/(d*d))*(1+x) * visc_comp_;
-                }
-            } else {
-                std::fill(output_mu, output_mu + n, viscosity_);
-                std::fill(output_dmudp, output_dmudp + n, 0.0);
-            }
-            std::fill(output_dmudr, output_dmudr + n, 0.0);
-        }
-
-        virtual void mu(const int n,
-                        const double* p,
-                        const double* /*r*/,
-                        const PhasePresence* /*cond*/,
-                        double* output_mu,
-                        double* output_dmudp,
-                        double* output_dmudr) const
-        {
-            if (visc_comp_) {
-                // #pragma omp parallel for
-                for (int i = 0; i < n; ++i) {
-                    // Computing a polynomial approximation to the exponential.
-                    double x = -visc_comp_*(p[i] - ref_press_);
-                    double d = (1.0 + x + 0.5*x*x);
-                    output_mu[i] = viscosity_/d;
-                    output_dmudp[i] = (viscosity_/(d*d))*(1+x) * visc_comp_;
-                }
-            } else {
-                std::fill(output_mu, output_mu + n, viscosity_);
-                std::fill(output_dmudp, output_dmudp + n, 0.0);
-            }
-            std::fill(output_dmudr, output_dmudr + n, 0.0);
-        }
-
-        virtual void B(const int n,
-                       const double* p,
-                       const double* /*z*/,
-                       double* output_B) const
-        {
-            if (comp_) {
-// #pragma omp parallel for
-                for (int i = 0; i < n; ++i) {
-                    // Computing a polynomial approximation to the exponential.
-                    double x = comp_*(p[i] - ref_press_);
-                    output_B[i] = ref_B_/(1.0 + x + 0.5*x*x);
-                }
-            } else {
-                std::fill(output_B, output_B + n, ref_B_);
-            }
-        }
-
-        virtual void dBdp(const int n,
-                          const double* p,
-                          const double* /*z*/,
-                          double* output_B,
-                          double* output_dBdp) const
-        {
-            if (comp_) {
-// #pragma omp parallel for
-                for (int i = 0; i < n; ++i) {
-                    double x = comp_*(p[i] - ref_press_);
-                    double d = (1.0 + x + 0.5*x*x);
-                    output_B[i] = ref_B_/d;
-                    output_dBdp[i] = (-ref_B_/(d*d))*(1 + x) * comp_;
-                }
-            } else {
-                std::fill(output_B, output_B + n, ref_B_);
-                std::fill(output_dBdp, output_dBdp + n, 0.0);
-            }
-        }
-
-        virtual void b(const int n,
-                               const double* p,
-                               const double* /*r*/,
-                               double* output_b,
-                               double* output_dbdp,
-                               double* output_dbdr) const
-                {
-                    if (comp_) {
-        // #pragma omp parallel for
-                        for (int i = 0; i < n; ++i) {
-                            // Computing a polynomial approximation to the exponential.
-                            double x = comp_*(p[i] - ref_press_);
-                            double d = (1.0 + x + 0.5*x*x);
-
-                            // b = 1/B = d/ref_B_B;
-                            output_b[i] = d/ref_B_;
-                            output_dbdp[i] = (1 + x) * comp_/ref_B_;
-                        }
-                    } else {
-
-                        std::fill(output_b, output_b + n, 1/ref_B_);
-                        std::fill(output_dbdp, output_dbdp + n, 0.0);
-                    }
-                    std::fill(output_dbdr, output_dbdr + n, 0.0);
-
-                }
-
-        virtual void b(const int n,
-                               const double* p,
-                               const double* /*r*/,
-                               const PhasePresence* /*cond*/,
-                               double* output_b,
-                               double* output_dbdp,
-                               double* output_dbdr) const
-                {
-                    if (comp_) {
-        // #pragma omp parallel for
-                        for (int i = 0; i < n; ++i) {
-                            // Computing a polynomial approximation to the exponential.
-                            double x = comp_*(p[i] - ref_press_);
-                            double d = (1.0 + x + 0.5*x*x);
-
-                            // b = 1/B = d/ref_B_B;
-                            output_b[i] = d/ref_B_;
-                            output_dbdp[i] = (1 + x) * comp_/ref_B_;
-                        }
-                    } else {
-
-                        std::fill(output_b, output_b + n, 1/ref_B_);
-                        std::fill(output_dbdp, output_dbdp + n, 0.0);
-                    }
-                    std::fill(output_dbdr, output_dbdr + n, 0.0);
-
-                }
-
-        virtual void rsSat(const int n,
-                             const double* /*p*/,
-                             double* output_rsSat,
-                             double* output_drsSatdp) const
-
-              {
-                  std::fill(output_rsSat, output_rsSat + n, 0.0);
-                  std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
-              }
-
-        virtual void rvSat(const int n,
-                             const double* /*p*/,
-                             double* output_rvSat,
-                             double* output_drvSatdp) const
-
-              {
-                  std::fill(output_rvSat, output_rvSat + n, 0.0);
-                  std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
-              }
-
-        virtual void R(const int n,
-                       const double* /*p*/,
-                       const double* /*z*/,
-                       double* output_R) const
-        {
-            std::fill(output_R, output_R + n, 0.0);
-        }
-
-        virtual void dRdp(const int n,
-                          const double* /*p*/,
-                          const double* /*z*/,
-                          double* output_R,
-                          double* output_dRdp) const
-        {
-            std::fill(output_R, output_R + n, 0.0);
-            std::fill(output_dRdp, output_dRdp + n, 0.0);
-        }
-
-    private:
-        double ref_press_;
-        double ref_B_;
-        double comp_;
-        double viscosity_;
-        double visc_comp_;
-    };
-
-}
-
-
-#endif // OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
-
diff --git a/opm/core/props/pvt/SinglePvtDeadSpline.cpp b/opm/core/props/pvt/SinglePvtDeadSpline.cpp
deleted file mode 100644
index 9f64197d..00000000
--- a/opm/core/props/pvt/SinglePvtDeadSpline.cpp
+++ /dev/null
@@ -1,226 +0,0 @@
-/*
-  Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
-
-  This file is part of the Open Porous Media project (OPM).
-
-  OPM is free software: you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation, either version 3 of the License, or
-  (at your option) any later version.
-
-  OPM is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-
-  You should have received a copy of the GNU General Public License
-  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#include "config.h"
-#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
-#include <opm/core/utility/buildUniformMonotoneTable.hpp>
-
-#include <opm/parser/eclipse/Utility/SingleRecordTable.hpp>
-
-#include <algorithm>
-
-// Extra includes for debug dumping of tables.
-// #include <boost/lexical_cast.hpp>
-// #include <string>
-// #include <fstream>
-
-namespace Opm
-{
-    //------------------------------------------------------------------------
-    // Member functions
-    //-------------------------------------------------------------------------
-
-    SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples)
-    {
-        int numRows = pvdoTable.numRows();
-
-        // Copy data
-        const std::vector<double> &press = pvdoTable.getPressureColumn();
-        const std::vector<double> &B = pvdoTable.getFormationFactorColumn();
-        const std::vector<double> &visc = pvdoTable.getViscosityColumn();
-
-        std::vector<double> B_inv(numRows);
-        for (int i = 0; i < numRows; ++i) {
-            B_inv[i] = 1.0 / B[i];
-        }
-
-        buildUniformMonotoneTable(press, B_inv, samples, b_);
-        buildUniformMonotoneTable(press, visc, samples, viscosity_);
-    }
-
-    SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples)
-    {
-        int numRows = pvdgTable.numRows();
-
-        // Copy data
-        const std::vector<double> &press = pvdgTable.getPressureColumn();
-        const std::vector<double> &B = pvdgTable.getFormationFactorColumn();
-        const std::vector<double> &visc = pvdgTable.getViscosityColumn();
-
-        std::vector<double> B_inv(numRows);
-        for (int i = 0; i < numRows; ++i) {
-            B_inv[i] = 1.0 / B[i];
-        }
-
-        buildUniformMonotoneTable(press, B_inv, samples, b_);
-        buildUniformMonotoneTable(press, visc, samples, viscosity_);
-    }
-
-    // Destructor
-    SinglePvtDeadSpline::~SinglePvtDeadSpline()
-    {
-    }
-
-
-
-    void SinglePvtDeadSpline::mu(const int n,
-                                 const double* p,
-                                 const double* /*z*/,
-                                 double* output_mu) const
-    {
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            output_mu[i] = viscosity_(p[i]);
-        }
-    }
-
-    void SinglePvtDeadSpline::mu(const int n,
-                               const double* p,
-                               const double* /*r*/,
-                               double* output_mu,
-                               double* output_dmudp,
-                               double* output_dmudr) const
-        {
-    // #pragma omp parallel for
-
-            for (int i = 0; i < n; ++i) {
-                output_mu[i] = viscosity_(p[i]);
-                output_dmudp[i] = viscosity_.derivative(p[i]);
-            }
-            std::fill(output_dmudr, output_dmudr + n, 0.0);
-
-        }
-
-    void SinglePvtDeadSpline::mu(const int n,
-                               const double* p,
-                               const double* /*r*/,
-                               const PhasePresence* /*cond*/,
-                               double* output_mu,
-                               double* output_dmudp,
-                               double* output_dmudr) const
-        {
-    // #pragma omp parallel for
-
-            for (int i = 0; i < n; ++i) {
-                output_mu[i] = viscosity_(p[i]);
-                output_dmudp[i] = viscosity_.derivative(p[i]);
-            }
-            std::fill(output_dmudr, output_dmudr + n, 0.0);
-
-        }
-
-    void SinglePvtDeadSpline::B(const int n,
-                                const double* p,
-                                const double* /*z*/,
-                                double* output_B) const
-    {
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            output_B[i] = 1.0/b_(p[i]);
-        }
-    }
-
-    void SinglePvtDeadSpline::dBdp(const int n,
-                                   const double* p,
-                                   const double* /*z*/,
-                                   double* output_B,
-                                   double* output_dBdp) const
-    {
-        B(n, p, 0, output_B);
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            double Bg = output_B[i];
-            output_dBdp[i] = -Bg*Bg*b_.derivative(p[i]);
-        }
-    }
-
-    void SinglePvtDeadSpline::b(const int n,
-                              const double* p,
-                              const double* /*r*/,
-                              double* output_b,
-                              double* output_dbdp,
-                              double* output_dbdr) const
-
-        {
-    // #pragma omp parallel for
-            for (int i = 0; i < n; ++i) {
-                output_b[i] = b_(p[i]);
-                output_dbdp[i] = b_.derivative(p[i]);
-
-            }
-            std::fill(output_dbdr, output_dbdr + n, 0.0);
-
-        }
-
-    void SinglePvtDeadSpline::b(const int n,
-                              const double* p,
-                              const double* /*r*/,
-                              const PhasePresence* /*cond*/,
-                              double* output_b,
-                              double* output_dbdp,
-                              double* output_dbdr) const
-
-        {
-    // #pragma omp parallel for
-            for (int i = 0; i < n; ++i) {
-                output_b[i] = b_(p[i]);
-                output_dbdp[i] = b_.derivative(p[i]);
-
-            }
-            std::fill(output_dbdr, output_dbdr + n, 0.0);
-
-        }
-
-    void SinglePvtDeadSpline::rsSat(const int n,
-                             const double* /*p*/,
-                             double* output_rsSat,
-                             double* output_drsSatdp) const
-    {
-        std::fill(output_rsSat, output_rsSat + n, 0.0);
-        std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
-    }
-
-    void SinglePvtDeadSpline::rvSat(const int n,
-                             const double* /*p*/,
-                             double* output_rvSat,
-                             double* output_drvSatdp) const
-    {
-        std::fill(output_rvSat, output_rvSat + n, 0.0);
-        std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
-    }
-
-    void SinglePvtDeadSpline::R(const int n,
-                                const double* /*p*/,
-                                const double* /*z*/,
-                                double* output_R) const
-    {
-        std::fill(output_R, output_R + n, 0.0);
-    }
-
-    void SinglePvtDeadSpline::dRdp(const int n,
-                                   const double* /*p*/,
-                                   const double* /*z*/,
-                                   double* output_R,
-                                   double* output_dRdp) const
-    {
-        std::fill(output_R, output_R + n, 0.0);
-        std::fill(output_dRdp, output_dRdp + n, 0.0);
-    }
-
-}
diff --git a/opm/core/props/pvt/SinglePvtInterface.cpp b/opm/core/props/pvt/SinglePvtInterface.cpp
deleted file mode 100644
index bbd905e6..00000000
--- a/opm/core/props/pvt/SinglePvtInterface.cpp
+++ /dev/null
@@ -1,48 +0,0 @@
-/*
-  Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
-
-  This file is part of the Open Porous Media project (OPM).
-
-  OPM is free software: you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation, either version 3 of the License, or
-  (at your option) any later version.
-
-  OPM is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-
-  You should have received a copy of the GNU General Public License
-  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#include "config.h"
-#include <opm/core/props/pvt/SinglePvtInterface.hpp>
-
-
-namespace Opm
-{
-
-    SinglePvtInterface::SinglePvtInterface()
-        : num_phases_(MaxNumPhases)
-    {
-        for (int i = 0; i < MaxNumPhases; ++i) {
-            phase_pos_[i] = i;
-        }
-    }
-
-    SinglePvtInterface::~SinglePvtInterface()
-    {
-    }
-
-    void SinglePvtInterface::setPhaseConfiguration(const int num_phases, const int* phase_pos)
-    {
-        num_phases_ = num_phases;
-        for (int i = 0; i < MaxNumPhases; ++i) {
-            phase_pos_[i] = phase_pos[i];
-        }
-    }
-
-
-} // namespace Opm
diff --git a/opm/core/props/pvt/SinglePvtLiveGas.cpp b/opm/core/props/pvt/SinglePvtLiveGas.cpp
deleted file mode 100644
index abc9a641..00000000
--- a/opm/core/props/pvt/SinglePvtLiveGas.cpp
+++ /dev/null
@@ -1,492 +0,0 @@
-//===========================================================================
-//
-// File: MiscibilityLiveGas.cpp
-//
-// Created: Wed Feb 10 09:21:53 2010
-//
-// Author: Bjørn Spjelkavik <bsp@sintef.no>
-//
-// Revision: $Id$
-//
-//===========================================================================
-/*
-  Copyright 2010 SINTEF ICT, Applied Mathematics.
-
-  This file is part of the Open Porous Media project (OPM).
-
-  OPM is free software: you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation, either version 3 of the License, or
-  (at your option) any later version.
-
-  OPM is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-
-  You should have received a copy of the GNU General Public License
-  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#include "config.h"
-#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
-#include <opm/core/utility/ErrorMacros.hpp>
-#include <opm/core/utility/linearInterpolation.hpp>
-
-#include <algorithm>
-#include <iostream>
-
-namespace Opm
-{
-
-    using Opm::linearInterpolation;
-    using Opm::linearInterpolationDerivative;
-
-
-    //------------------------------------------------------------------------
-    // Member functions
-    //-------------------------------------------------------------------------
-    SinglePvtLiveGas::SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable)
-    {
-        // GAS, PVTG
-        saturated_gas_table_.resize(4);
-        saturated_gas_table_[0] = pvtgTable.getOuterTable()->getPressureColumn();
-        saturated_gas_table_[1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn();
-        saturated_gas_table_[2] = pvtgTable.getOuterTable()->getGasViscosityColumn();
-        saturated_gas_table_[3] = pvtgTable.getOuterTable()->getOilSolubilityColumn();
-
-        int sz = pvtgTable.getOuterTable()->numRows();
-        undersat_gas_tables_.resize(sz);
-        for (int i=0; i<sz; ++i) {
-            const auto &undersatTable = *pvtgTable.getInnerTable(i);
-
-            undersat_gas_tables_[i].resize(3);
-            undersat_gas_tables_[i][0] = undersatTable.getOilSolubilityColumn();
-            undersat_gas_tables_[i][1] = undersatTable.getGasFormationFactorColumn();
-            undersat_gas_tables_[i][2] = undersatTable.getGasViscosityColumn();
-        }
-        
-        // Bg -> 1/Bg      
-        for (int i=0; i<sz; ++i) {
-            saturated_gas_table_[1][i] = 1.0/saturated_gas_table_[1][i];
-            for (size_t j=0; j<undersat_gas_tables_[i][1].size(); ++j) {
-                undersat_gas_tables_[i][1][j] = 1.0/undersat_gas_tables_[i][1][j];
-            }
-        }
-    }
-
-    // Destructor
-     SinglePvtLiveGas::~SinglePvtLiveGas()
-    {
-    }
-
-
-    void SinglePvtLiveGas::mu(const int n,
-                              const double* p,
-                              const double* z,
-                              double* output_mu) const
-    {
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            output_mu[i] = miscible_gas(p[i], z + num_phases_*i, 2, false);
-        }
-    }
-
-    /// Viscosity and its derivatives as a function of p and r.
-    void SinglePvtLiveGas::mu(const int /*n*/,
-                              const double* /*p*/,
-                              const double* /*r*/,
-                              double* /*output_mu*/,
-                              double* /*output_dmudp*/,
-                              double* /*output_dmudr*/) const
-    {
-        OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
-    }
-
-    /// Viscosity and its derivatives as a function of p and r.
-    void SinglePvtLiveGas::mu(const int n,
-                               const double* p,
-                               const double* r,
-                               const PhasePresence* cond,
-                               double* output_mu,
-                               double* output_dmudp,
-                               double* output_dmudr) const
-    {
-        for (int i = 0; i < n; ++i) {
-            const PhasePresence& cnd = cond[i];
-            output_mu[i] = miscible_gas(p[i], r[i], cnd,2, 0);
-            output_dmudp[i] = miscible_gas(p[i], r[i], cnd, 2, 1);
-            output_dmudr[i] = miscible_gas(p[i], r[i], cnd, 2, 2);
-        }
-
-    }
-
-
-    /// Formation volume factor as a function of p and z.
-    void SinglePvtLiveGas::B(const int n,
-                             const double* p,
-                             const double* z,
-                             double* output_B) const
-    {
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            output_B[i] = evalB(p[i], z + num_phases_*i);
-        }
-
-    }
-
-
-    /// Formation volume factor and p-derivative as functions of p and z.
-    void SinglePvtLiveGas::dBdp(const int n,
-                                const double* p,
-                                const double* z,
-                                double* output_B,
-                                double* output_dBdp) const
-    {
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            evalBDeriv(p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
-        }
-    }
-
-    /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
-    void SinglePvtLiveGas::b(const int /*n*/,
-                             const double* /*p*/,
-                             const double* /*r*/,
-                             double* /*output_b*/,
-                             double* /*output_dbdp*/,
-                             double* /*output_dbdr*/) const
-
-    {
-        OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
-    }
-
-    /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
-    void SinglePvtLiveGas::b(const int n,
-                          const double* p,
-                          const double* r,
-                          const PhasePresence* cond,
-                          double* output_b,
-                          double* output_dbdp,
-                          double* output_dbdr) const
-
-    {
-        // #pragma omp parallel for
-                for (int i = 0; i < n; ++i) {
-                    const PhasePresence& cnd = cond[i];
-
-                    output_b[i] = miscible_gas(p[i], r[i], cnd, 1, 0);
-                    output_dbdp[i] = miscible_gas(p[i], r[i], cnd, 1, 1);
-                    output_dbdr[i] = miscible_gas(p[i], r[i], cnd, 1, 2);
-
-                }
-    }
-
-    /// Gas resolution and its derivatives at bublepoint as a function of p.
-    void SinglePvtLiveGas::rvSat(const int n,
-                             const double* p,
-                             double* output_rvSat,
-                             double* output_drvSatdp) const
-    {
-        for (int i = 0; i < n; ++i) {
-            output_rvSat[i] = linearInterpolation(saturated_gas_table_[0],
-                    saturated_gas_table_[3],p[i]);
-            output_drvSatdp[i] = linearInterpolationDerivative(saturated_gas_table_[0],
-                    saturated_gas_table_[3],p[i]);
-
-        }
-    }
-
-    void SinglePvtLiveGas::rsSat(const int n,
-                             const double* /*p*/,
-                             double* output_rsSat,
-                             double* output_drsSatdp) const
-    {
-        std::fill(output_rsSat, output_rsSat + n, 0.0);
-        std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
-    }
-
-    /// Solution factor as a function of p and z.
-    void SinglePvtLiveGas::R(const int n,
-                             const double* p,
-                             const double* z,
-                             double* output_R) const
-    {
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            output_R[i] = evalR(p[i], z + num_phases_*i);
-        }
-
-    }
-
-
-    /// Solution factor and p-derivative as functions of p and z.
-    void SinglePvtLiveGas::dRdp(const int n,
-                                const double* p,
-                                const double* z,
-                                double* output_R,
-                                double* output_dRdp) const
-    {
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            evalRDeriv(p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
-        }
-    }
-
-
-    // ---- Private methods ----
-
-    double SinglePvtLiveGas::evalB(const double press, const double* surfvol) const
-    {
-        if (surfvol[phase_pos_[Vapour]] == 0.0) {
-            // To handle no-gas case.
-            return 1.0;
-        }
-        return 1.0/miscible_gas(press, surfvol, 1, false);
-    }
-
-    void SinglePvtLiveGas::evalBDeriv(const double press, const double* surfvol,
-                                      double& Bval, double& dBdpval) const
-    {
-        if (surfvol[phase_pos_[Vapour]] == 0.0) {
-            // To handle no-gas case.
-            Bval = 1.0;
-            dBdpval = 0.0;
-            return;
-        }
-        Bval = evalB(press, surfvol);
-        dBdpval =  -Bval*Bval*miscible_gas(press, surfvol, 1, true);
-    }
-
-    double SinglePvtLiveGas::evalR(const double press, const double* surfvol) const
-    {
-        if (surfvol[phase_pos_[Liquid]] == 0.0) {
-            // To handle no-gas case.
-            return 0.0;
-        }
-        double satR = linearInterpolation(saturated_gas_table_[0],
-                                             saturated_gas_table_[3], press);
-        double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
-        if (satR < maxR ) {
-            // Saturated case
-            return satR;
-        } else {
-            // Undersaturated case
-            return maxR;
-        }
-    }
-
-    void SinglePvtLiveGas::evalRDeriv(const double press, const double* surfvol,
-                                      double& Rval, double& dRdpval) const
-    {
-        if (surfvol[phase_pos_[Liquid]] == 0.0) {
-            // To handle no-gas case.
-            Rval = 0.0;
-            dRdpval = 0.0;
-            return;
-        }
-        double satR = linearInterpolation(saturated_gas_table_[0],
-                                             saturated_gas_table_[3], press);
-        double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
-        if (satR < maxR ) {
-            // Saturated case
-            Rval = satR;
-            dRdpval = linearInterpolationDerivative(saturated_gas_table_[0],
-                                            saturated_gas_table_[3],
-                                            press);
-        } else {
-            // Undersaturated case
-            Rval = maxR;
-            dRdpval = 0.0;
-        }
-    }
-
-    double SinglePvtLiveGas::miscible_gas(const double press,
-                                          const double* surfvol,
-                                          const int item,
-                                          const bool deriv) const
-    {
-        int section;
-        double Rval = linearInterpolation(saturated_gas_table_[0],
-                                                saturated_gas_table_[3], press,
-                                                section);
-        double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
-        if (deriv) {
-            if (Rval < maxR ) {  // Saturated case
-                return linearInterpolationDerivative(saturated_gas_table_[0],
-                                                saturated_gas_table_[item],
-                                                press);
-            } else {  // Undersaturated case
-                int is = section;
-                if (undersat_gas_tables_[is][0].size() < 2) {
-                    double val = (saturated_gas_table_[item][is+1]
-                                  - saturated_gas_table_[item][is]) /
-                        (saturated_gas_table_[0][is+1] -
-                         saturated_gas_table_[0][is]);
-                    return val;
-                }
-                double val1 =
-                    linearInterpolation(undersat_gas_tables_[is][0],
-                                              undersat_gas_tables_[is][item],
-                                              maxR);
-                double val2 =
-                    linearInterpolation(undersat_gas_tables_[is+1][0],
-                                              undersat_gas_tables_[is+1][item],
-                                              maxR);
-                double val = (val2 - val1)/
-                    (saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
-                return val;
-            }
-        } else {
-            if (Rval < maxR ) {  // Saturated case
-                return linearInterpolation(saturated_gas_table_[0],
-                                                 saturated_gas_table_[item],
-                                                 press);
-            } else {  // Undersaturated case
-                int is = section;
-                // Extrapolate from first table section
-                if (is == 0 && press < saturated_gas_table_[0][0]) {
-                    return linearInterpolation(undersat_gas_tables_[0][0],
-                                                     undersat_gas_tables_[0][item],
-                                                     maxR);
-                }
-
-                // Extrapolate from last table section
-                int ltp = saturated_gas_table_[0].size() - 1;
-                if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
-                    return linearInterpolation(undersat_gas_tables_[ltp][0],
-                                                     undersat_gas_tables_[ltp][item],
-                                                     maxR);
-                }
-
-                // Interpolate between table sections
-                double w = (press - saturated_gas_table_[0][is]) /
-                    (saturated_gas_table_[0][is+1] -
-                     saturated_gas_table_[0][is]);
-                if (undersat_gas_tables_[is][0].size() < 2) {
-                    double val = saturated_gas_table_[item][is] +
-                        w*(saturated_gas_table_[item][is+1] -
-                           saturated_gas_table_[item][is]);
-                    return val;
-                }
-                double val1 =
-                    linearInterpolation(undersat_gas_tables_[is][0],
-                                              undersat_gas_tables_[is][item],
-                                              maxR);
-                double val2 =
-                    linearInterpolation(undersat_gas_tables_[is+1][0],
-                                              undersat_gas_tables_[is+1][item],
-                                              maxR);
-                double val = val1 + w*(val2 - val1);
-                return val;
-            }
-        }
-    }
-
-    double SinglePvtLiveGas::miscible_gas(const double press,
-                                          const double r,
-                                          const PhasePresence& cond,
-                                          const int item,
-                                          const int deriv) const
-{
-    const bool isSat = cond.hasFreeOil();
-
-        // Derivative w.r.t p
-        if (deriv == 1) {
-            if (isSat) {  // Saturated case
-                return linearInterpolationDerivative(saturated_gas_table_[0],
-                                                saturated_gas_table_[item],
-                                                press);
-            } else {  // Undersaturated case
-                int is = tableIndex(saturated_gas_table_[0], press);
-                if (undersat_gas_tables_[is][0].size() < 2) {
-                    double val = (saturated_gas_table_[item][is+1]
-                                  - saturated_gas_table_[item][is]) /
-                        (saturated_gas_table_[0][is+1] -
-                         saturated_gas_table_[0][is]);
-                    return val;
-                }
-                double val1 =
-                    linearInterpolation(undersat_gas_tables_[is][0],
-                                              undersat_gas_tables_[is][item],
-                                              r);
-                double val2 =
-                    linearInterpolation(undersat_gas_tables_[is+1][0],
-                                              undersat_gas_tables_[is+1][item],
-                                              r);
-                double val = (val2 - val1)/
-                    (saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
-                return val;
-            }
-        } else if (deriv == 2){
-            if (isSat) {
-                return 0;
-            } else {
-                int is = tableIndex(saturated_gas_table_[0], press);
-                double w = (press - saturated_gas_table_[0][is]) /
-                    (saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
-                assert(undersat_gas_tables_[is][0].size() >= 2);
-                assert(undersat_gas_tables_[is+1][0].size() >= 2);
-                double val1 =
-                    linearInterpolationDerivative(undersat_gas_tables_[is][0],
-                                              undersat_gas_tables_[is][item],
-                                              r);
-                double val2 =
-                    linearInterpolationDerivative(undersat_gas_tables_[is+1][0],
-                                              undersat_gas_tables_[is+1][item],
-                                              r);
-
-                double val = val1 + w * (val2 - val1);
-                return val;
-
-            }
-        } else {
-            if (isSat) {  // Saturated case
-                return linearInterpolation(saturated_gas_table_[0],
-                                                 saturated_gas_table_[item],
-                                                 press);
-            } else {  // Undersaturated case
-                int is = tableIndex(saturated_gas_table_[0], press);
-                // Extrapolate from first table section
-                if (is == 0 && press < saturated_gas_table_[0][0]) {
-                    return linearInterpolation(undersat_gas_tables_[0][0],
-                                                     undersat_gas_tables_[0][item],
-                                                     r);
-                }
-
-                // Extrapolate from last table section
-                //int ltp = saturated_gas_table_[0].size() - 1;
-                //if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
-                //    return linearInterpolation(undersat_gas_tables_[ltp][0],
-                //                                    undersat_gas_tables_[ltp][item],
-                //                                    r);
-                //}
-
-                // Interpolate between table sections
-                double w = (press - saturated_gas_table_[0][is]) /
-                    (saturated_gas_table_[0][is+1] -
-                     saturated_gas_table_[0][is]);
-                if (undersat_gas_tables_[is][0].size() < 2) {
-                    double val = saturated_gas_table_[item][is] +
-                        w*(saturated_gas_table_[item][is+1] -
-                           saturated_gas_table_[item][is]);
-                    return val;
-                }
-                double val1 =
-                    linearInterpolation(undersat_gas_tables_[is][0],
-                                              undersat_gas_tables_[is][item],
-                                              r);
-                double val2 =
-                    linearInterpolation(undersat_gas_tables_[is+1][0],
-                                              undersat_gas_tables_[is+1][item],
-                                              r);
-                double val = val1 + w*(val2 - val1);
-                return val;
-            }
-        }
-    }
-
-
-
-
-} // namespace Opm
diff --git a/opm/core/props/pvt/SinglePvtLiveOil.cpp b/opm/core/props/pvt/SinglePvtLiveOil.cpp
deleted file mode 100644
index b7e4b113..00000000
--- a/opm/core/props/pvt/SinglePvtLiveOil.cpp
+++ /dev/null
@@ -1,551 +0,0 @@
-/*
-  Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
-
-  This file is part of the Open Porous Media project (OPM).
-
-  OPM is free software: you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation, either version 3 of the License, or
-  (at your option) any later version.
-
-  OPM is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-
-  You should have received a copy of the GNU General Public License
-  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#include "config.h"
-#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
-#include <opm/core/utility/ErrorMacros.hpp>
-#include <opm/core/utility/linearInterpolation.hpp>
-
-#include <algorithm>
-
-namespace Opm
-{
-
-    using Opm::linearInterpolation;
-    using Opm::linearInterpolationDerivative;
-    using Opm::tableIndex;
-
-    //------------------------------------------------------------------------
-    // Member functions
-    //-------------------------------------------------------------------------
-
-    SinglePvtLiveOil::SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable)
-    {
-        const auto saturatedPvto = pvtoTable.getOuterTable();
-
-        // OIL, PVTO
-        saturated_oil_table_.resize(4);
-        const int sz =  saturatedPvto->numRows();
-        for (int k=0; k<4; ++k) {
-            saturated_oil_table_[k].resize(sz);
-        }
-        for (int i=0; i<sz; ++i) {
-            saturated_oil_table_[0][i] = saturatedPvto->getPressureColumn()[i]; // p
-            saturated_oil_table_[1][i] = 1.0/saturatedPvto->getOilFormationFactorColumn()[i]; // 1/Bo
-            saturated_oil_table_[2][i] = saturatedPvto->getOilViscosityColumn()[i]; // mu_o
-            saturated_oil_table_[3][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs
-        }
-
-        undersat_oil_tables_.resize(sz);
-        for (int i=0; i<sz; ++i) {
-            const auto undersaturatedPvto = pvtoTable.getInnerTable(i);
-
-            undersat_oil_tables_[i].resize(3);
-            int tsize = undersaturatedPvto->numRows();
-            undersat_oil_tables_[i][0].resize(tsize);
-            undersat_oil_tables_[i][1].resize(tsize);
-            undersat_oil_tables_[i][2].resize(tsize);
-            for (int j=0; j<tsize; ++j) {
-                undersat_oil_tables_[i][0][j] = undersaturatedPvto->getPressureColumn()[j];  // p
-                undersat_oil_tables_[i][1][j] = 1.0/undersaturatedPvto->getOilFormationFactorColumn()[j];  // 1/Bo
-                undersat_oil_tables_[i][2][j] = undersaturatedPvto->getOilViscosityColumn()[j];  // mu_o
-            }
-        }
-
-        // Complete undersaturated tables by extrapolating from existing data
-        // as is done in Eclipse and Mrst
-        int iNext = -1;
-        for (int i=0; i<sz; ++i) {
-            // Skip records already containing undersaturated data
-            if (undersat_oil_tables_[i][0].size() > 1) {
-                continue;
-            }
-            // Look ahead for next record containing undersaturated data
-            if (iNext < i) {
-                iNext = i+1;
-                while (iNext<sz && undersat_oil_tables_[iNext][0].size() < 2) {
-                    ++iNext;
-                }
-                if (iNext == sz) OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
-            }
-            // Add undersaturated data to current record while maintaining compressibility and viscosibility
-            typedef std::vector<std::vector<std::vector<double> > >::size_type sz_t;
-            for (sz_t j=1; j<undersat_oil_tables_[iNext][0].size(); ++j) {
-                double diffPressure = undersat_oil_tables_[iNext][0][j]-undersat_oil_tables_[iNext][0][j-1];
-                double pressure = undersat_oil_tables_[i][0].back()+diffPressure;
-                undersat_oil_tables_[i][0].push_back(pressure);
-                double compr = (1.0/undersat_oil_tables_[iNext][1][j]-1.0/undersat_oil_tables_[iNext][1][j-1])
-                        / (0.5*(1.0/undersat_oil_tables_[iNext][1][j]+1.0/undersat_oil_tables_[iNext][1][j-1]));
-                double B = (1.0/undersat_oil_tables_[i][1].back())*(1.0+0.5*compr)/(1.0-0.5*compr);
-                undersat_oil_tables_[i][1].push_back(1.0/B);
-                double visc = (undersat_oil_tables_[iNext][2][j]-undersat_oil_tables_[iNext][2][j-1])
-                        / (0.5*(undersat_oil_tables_[iNext][2][j]+undersat_oil_tables_[iNext][2][j-1]));
-                double mu = (undersat_oil_tables_[i][2].back())*(1.0+0.5*visc)/(1.0-0.5*visc);
-                undersat_oil_tables_[i][2].push_back(mu);
-            }
-        }
-    }
-
-    /// Destructor.
-    SinglePvtLiveOil::~SinglePvtLiveOil()
-    {
-    }
-
-
-    /// Viscosity as a function of p and z.
-    void SinglePvtLiveOil::mu(const int n,
-                              const double* p,
-                              const double* z,
-                              double* output_mu) const
-    {
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            output_mu[i] = miscible_oil(p[i], z + num_phases_*i, 2, false);
-        }
-    }
-
-    /// Viscosity and its derivatives as a function of p and r.
-    void SinglePvtLiveOil::mu(const int n,
-                               const double* p,
-                               const double* r,
-                               double* output_mu,
-                               double* output_dmudp,
-                               double* output_dmudr) const
-    {
-        // #pragma omp parallel for
-                for (int i = 0; i < n; ++i) {
-                    output_mu[i] = miscible_oil(p[i], r[i], 2, 0);
-                    output_dmudp[i] = miscible_oil(p[i], r[i], 2, 1);
-                    output_dmudr[i] = miscible_oil(p[i], r[i], 2, 2);
-
-                }
-    }
-
-    /// Viscosity and its derivatives as a function of p and r.
-    void SinglePvtLiveOil::mu(const int n,
-                               const double* p,
-                               const double* r,
-                               const PhasePresence* cond,
-                               double* output_mu,
-                               double* output_dmudp,
-                               double* output_dmudr) const
-    {
-        // #pragma omp parallel for
-                for (int i = 0; i < n; ++i) {
-                    const PhasePresence& cnd = cond[i];
-
-                    output_mu[i] = miscible_oil(p[i], r[i], cnd, 2, 0);
-                    output_dmudp[i] = miscible_oil(p[i], r[i], cnd, 2, 1);
-                    output_dmudr[i] = miscible_oil(p[i], r[i], cnd, 2, 2);
-
-                }
-    }
-
-
-    /// Formation volume factor as a function of p and z.
-    void SinglePvtLiveOil::B(const int n,
-                             const double* p,
-                             const double* z,
-                             double* output_B) const
-    {
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            output_B[i] = evalB(p[i], z + num_phases_*i);
-        }
-
-    }
-
-
-    /// Formation volume factor and p-derivative as functions of p and z.
-    void SinglePvtLiveOil::dBdp(const int n,
-                                const double* p,
-                                const double* z,
-                                double* output_B,
-                                double* output_dBdp) const
-    {
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            evalBDeriv(p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
-        }
-    }
-
-    void SinglePvtLiveOil::b(const int n,
-                          const double* p,
-                          const double* r,
-                          double* output_b,
-                          double* output_dbdp,
-                          double* output_dbdr) const
-
-    {
-        // #pragma omp parallel for
-                for (int i = 0; i < n; ++i) {
-                    output_b[i] = miscible_oil(p[i], r[i], 1, 0);
-                    output_dbdp[i] = miscible_oil(p[i], r[i], 1, 1);
-                    output_dbdr[i] = miscible_oil(p[i], r[i], 1, 2);
-
-                }
-    }
-
-    void SinglePvtLiveOil::b(const int n,
-                          const double* p,
-                          const double* r,
-                          const PhasePresence* cond,
-                          double* output_b,
-                          double* output_dbdp,
-                          double* output_dbdr) const
-
-    {
-        // #pragma omp parallel for
-                for (int i = 0; i < n; ++i) {
-                    const PhasePresence& cnd = cond[i];
-
-                    output_b[i] = miscible_oil(p[i], r[i], cnd, 1, 0);
-                    output_dbdp[i] = miscible_oil(p[i], r[i], cnd, 1, 1);
-                    output_dbdr[i] = miscible_oil(p[i], r[i], cnd, 1, 2);
-
-                }
-    }
-
-    void SinglePvtLiveOil::rsSat(const int n,
-                             const double* p,
-                             double* output_rsSat,
-                             double* output_drsSatdp) const
-    {
-
-        for (int i = 0; i < n; ++i) {
-            output_rsSat[i] = linearInterpolation(saturated_oil_table_[0],
-                    saturated_oil_table_[3],p[i]);
-            output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[0],
-                    saturated_oil_table_[3],p[i]);
-
-        }
-    }
-
-    void SinglePvtLiveOil::rvSat(const int n,
-                             const double* /*p*/,
-                             double* output_rvSat,
-                             double* output_drvSatdp) const
-    {
-        std::fill(output_rvSat, output_rvSat + n, 0.0);
-        std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
-    }
-
-    /// Solution factor as a function of p and z.
-    void SinglePvtLiveOil::R(const int n,
-                             const double* p,
-                             const double* z,
-                             double* output_R) const
-    {
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            output_R[i] = evalR(p[i], z + num_phases_*i);
-        }
-
-    }
-
-
-    /// Solution factor and p-derivative as functions of p and z.
-    void SinglePvtLiveOil::dRdp(const int n,
-                                const double* p,
-                                const double* z,
-                                double* output_R,
-                                double* output_dRdp) const
-    {
-// #pragma omp parallel for
-        for (int i = 0; i < n; ++i) {
-            evalRDeriv(p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
-        }
-    }
-
-
-
-
-    // ---- Private methods ----
-
-    double SinglePvtLiveOil::evalB(double press, const double* surfvol) const
-    {
-        // if (surfvol[phase_pos_[Liquid]] == 0.0) return 1.0; // To handle no-oil case.
-        return 1.0/miscible_oil(press, surfvol, 1, false);
-    }
-
-
-    void SinglePvtLiveOil::evalBDeriv(const double press, const double* surfvol,
-                                      double& Bval, double& dBdpval) const
-    {
-        Bval = evalB(press, surfvol);
-        dBdpval = -Bval*Bval*miscible_oil(press, surfvol, 1, true);
-    }
-
-    double SinglePvtLiveOil::evalR(double press, const double* surfvol) const
-    {
-        if (surfvol[phase_pos_[Vapour]] == 0.0) {
-            return 0.0;
-        }
-        double Rval = linearInterpolation(saturated_oil_table_[0],
-                                             saturated_oil_table_[3], press);
-        double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
-        if (Rval < maxR ) {  // Saturated case
-            return Rval;
-        } else {
-            return maxR;  // Undersaturated case
-        }
-    }
-
-    void SinglePvtLiveOil::evalRDeriv(const double press, const double* surfvol,
-                                      double& Rval, double& dRdpval) const
-    {
-        if (surfvol[phase_pos_[Vapour]] == 0.0) {
-            Rval = 0.0;
-            dRdpval = 0.0;
-            return;
-        }
-        Rval = linearInterpolation(saturated_oil_table_[0],
-                                      saturated_oil_table_[3], press);
-        double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
-        if (Rval < maxR ) {
-            // Saturated case
-            dRdpval = linearInterpolationDerivative(saturated_oil_table_[0],
-                                            saturated_oil_table_[3],
-                                            press);
-        } else {
-            // Undersaturated case
-            Rval = maxR;
-            dRdpval = 0.0;
-        }
-    }
-
-
-    double SinglePvtLiveOil::miscible_oil(const double press,
-                                          const double* surfvol,
-                                          const int item,
-                                          const bool deriv) const
-    {
-        int section;
-        double Rval = linearInterpolation(saturated_oil_table_[0],
-                                                saturated_oil_table_[3],
-                                                press, section);
-        double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
-        if (deriv) {
-            if (Rval < maxR ) {  // Saturated case
-                return linearInterpolationDerivative(saturated_oil_table_[0],
-                                                saturated_oil_table_[item],
-                                                press);
-            } else {  // Undersaturated case
-                int is = tableIndex(saturated_oil_table_[3], maxR);
-                double w = (maxR - saturated_oil_table_[3][is]) /
-                    (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
-                assert(undersat_oil_tables_[is][0].size() >= 2);
-                assert(undersat_oil_tables_[is+1][0].size() >= 2);
-                double val1 =
-                    linearInterpolationDerivative(undersat_oil_tables_[is][0],
-                                             undersat_oil_tables_[is][item],
-                                             press);
-                double val2 =
-                    linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
-                                             undersat_oil_tables_[is+1][item],
-                                             press);
-                double val = val1 + w*(val2 - val1);
-                return val;
-            }
-        } else {
-            if (Rval < maxR ) {  // Saturated case
-                return linearInterpolation(saturated_oil_table_[0],
-                                                 saturated_oil_table_[item],
-                                                 press);
-            } else {  // Undersaturated case
-                // Interpolate between table sections
-                int is = tableIndex(saturated_oil_table_[3], maxR);
-                double w = (maxR - saturated_oil_table_[3][is]) /
-                    (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
-                assert(undersat_oil_tables_[is][0].size() >= 2);
-                assert(undersat_oil_tables_[is+1][0].size() >= 2);
-                double val1 =
-                    linearInterpolation(undersat_oil_tables_[is][0],
-                                              undersat_oil_tables_[is][item],
-                                              press);
-                double val2 =
-                    linearInterpolation(undersat_oil_tables_[is+1][0],
-                                              undersat_oil_tables_[is+1][item],
-                                              press);
-                double val = val1 + w*(val2 - val1);
-                return val;
-            }
-        }
-    }
-
-    double SinglePvtLiveOil::miscible_oil(const double press,
-                                          const double r,
-                                          const int item,
-                                          const int deriv) const
-    {
-        int section;
-        double Rval = linearInterpolation(saturated_oil_table_[0],
-                                                saturated_oil_table_[3],
-                                                press, section);
-        // derivative with respect to frist component (pressure)
-        if (deriv == 1) {
-            if (Rval <= r ) {  // Saturated case
-                return linearInterpolationDerivative(saturated_oil_table_[0],
-                                                saturated_oil_table_[item],
-                                                press);
-            } else {  // Undersaturated case
-                int is = tableIndex(saturated_oil_table_[3], r);
-                double w = (r - saturated_oil_table_[3][is]) /
-                    (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
-                assert(undersat_oil_tables_[is][0].size() >= 2);
-                assert(undersat_oil_tables_[is+1][0].size() >= 2);
-                double val1 =
-                    linearInterpolationDerivative(undersat_oil_tables_[is][0],
-                                             undersat_oil_tables_[is][item],
-                                             press);
-                double val2 =
-                    linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
-                                             undersat_oil_tables_[is+1][item],
-                                             press);
-                double val = val1 + w*(val2 - val1);
-                return val;
-            }
-            // derivative with respect to second component (r)
-        } else if (deriv == 2)  {
-            if (Rval <= r ) {  // Saturated case
-                return 0;
-            } else {  // Undersaturated case
-                int is = tableIndex(saturated_oil_table_[3], r);
-                assert(undersat_oil_tables_[is][0].size() >= 2);
-                assert(undersat_oil_tables_[is+1][0].size() >= 2);
-                double val1 =
-                    linearInterpolation(undersat_oil_tables_[is][0],
-                                              undersat_oil_tables_[is][item],
-                                              press);
-                double val2 =
-                    linearInterpolation(undersat_oil_tables_[is+1][0],
-                                              undersat_oil_tables_[is+1][item],
-                                              press);
-
-                double val = (val2 - val1)/(saturated_oil_table_[3][is+1]-saturated_oil_table_[3][is]);
-                return val;
-            }
-
-
-            }    else {
-            if (Rval <= r ) {  // Saturated case
-                return linearInterpolation(saturated_oil_table_[0],
-                                                 saturated_oil_table_[item],
-                                                 press);
-            } else {  // Undersaturated case
-                // Interpolate between table sections
-                int is = tableIndex(saturated_oil_table_[3], r);
-                double w = (r - saturated_oil_table_[3][is]) /
-                    (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
-                assert(undersat_oil_tables_[is][0].size() >= 2);
-                assert(undersat_oil_tables_[is+1][0].size() >= 2);
-                double val1 =
-                    linearInterpolation(undersat_oil_tables_[is][0],
-                                              undersat_oil_tables_[is][item],
-                                              press);
-                double val2 =
-                    linearInterpolation(undersat_oil_tables_[is+1][0],
-                                              undersat_oil_tables_[is+1][item],
-                                              press);
-                double val = val1 + w*(val2 - val1);
-                return val;
-            }
-        }
-    }
-
-    double SinglePvtLiveOil::miscible_oil(const double press,
-                                          const double r,
-                                          const PhasePresence& cond,
-                                          const int item,
-                                          const int deriv) const
-    {
-        const bool isSat = cond.hasFreeGas();
-
-        // derivative with respect to frist component (pressure)
-        if (deriv == 1) {
-            if (isSat) {  // Saturated case
-                return linearInterpolationDerivative(saturated_oil_table_[0],
-                                                saturated_oil_table_[item],
-                                                press);
-            } else {  // Undersaturated case
-                int is = tableIndex(saturated_oil_table_[3], r);
-                double w = (r - saturated_oil_table_[3][is]) /
-                    (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
-                assert(undersat_oil_tables_[is][0].size() >= 2);
-                assert(undersat_oil_tables_[is+1][0].size() >= 2);
-                double val1 =
-                    linearInterpolationDerivative(undersat_oil_tables_[is][0],
-                                             undersat_oil_tables_[is][item],
-                                             press);
-                double val2 =
-                    linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
-                                             undersat_oil_tables_[is+1][item],
-                                             press);
-                double val = val1 + w*(val2 - val1);
-                return val;
-            }
-            // derivative with respect to second component (r)
-        } else if (deriv == 2)  {
-            if (isSat) {  // Saturated case
-                return 0;
-            } else {  // Undersaturated case
-                int is = tableIndex(saturated_oil_table_[3], r);
-                assert(undersat_oil_tables_[is][0].size() >= 2);
-                assert(undersat_oil_tables_[is+1][0].size() >= 2);
-                double val1 =
-                    linearInterpolation(undersat_oil_tables_[is][0],
-                                              undersat_oil_tables_[is][item],
-                                              press);
-                double val2 =
-                    linearInterpolation(undersat_oil_tables_[is+1][0],
-                                              undersat_oil_tables_[is+1][item],
-                                              press);
-
-                double val = (val2 - val1)/(saturated_oil_table_[3][is+1]-saturated_oil_table_[3][is]);
-                return val;
-            }
-
-
-            }    else {
-            if (isSat) {  // Saturated case
-                return linearInterpolation(saturated_oil_table_[0],
-                                                 saturated_oil_table_[item],
-                                                 press);
-            } else {  // Undersaturated case
-                // Interpolate between table sections
-                int is = tableIndex(saturated_oil_table_[3], r);
-                double w = (r - saturated_oil_table_[3][is]) /
-                    (saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
-                assert(undersat_oil_tables_[is][0].size() >= 2);
-                assert(undersat_oil_tables_[is+1][0].size() >= 2);
-                double val1 =
-                    linearInterpolation(undersat_oil_tables_[is][0],
-                                              undersat_oil_tables_[is][item],
-                                              press);
-                double val2 =
-                    linearInterpolation(undersat_oil_tables_[is+1][0],
-                                              undersat_oil_tables_[is+1][item],
-                                              press);
-                double val = val1 + w*(val2 - val1);
-                return val;
-            }
-        }
-    }
-
-} // namespace Opm
diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp
index dd90e368..e97d5077 100644
--- a/opm/core/simulator/EquilibrationHelpers.hpp
+++ b/opm/core/simulator/EquilibrationHelpers.hpp
@@ -139,7 +139,7 @@ namespace Opm
                 props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp);
 
                 std::vector<double> rho(np, 0.0);
-                props_.density(1, &A[0], &rho[0]);
+                props_.density(1, &A[0], &c_[0], &rho[0]);
 
                 return rho;
             }
diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp
index 1cf88ee4..a4ec1ee4 100644
--- a/opm/core/simulator/initStateEquil.hpp
+++ b/opm/core/simulator/initStateEquil.hpp
@@ -279,7 +279,7 @@ namespace Opm
                             if (rec[i].live_oil_table_index > 0) {
                                 if (deck->hasKeyword("RSVD") &&
                                     size_t(rec[i].live_oil_table_index) <= deck->getKeyword("RSVD")->size()) { 
-                                    Opm::SingleRecordTable rsvd(deck->getKeyword("RSVD"),std::vector<std::string>{"vd", "rs"},rec[i].live_oil_table_index-1);                                
+                                    Opm::SingleRecordTable rsvd(deck->getKeyword("RSVD"),std::vector<std::string>{"vd", "rs"},rec[i].live_oil_table_index-1);
                                     std::vector<double> vd(rsvd.getColumn("vd"));
                                     std::vector<double> rs(rsvd.getColumn("rs"));
                                     rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, vd, rs));
@@ -369,8 +369,6 @@ namespace Opm
                                  const UnstructuredGrid&                 G    ,
                                  const double grav)
                 {
-                    typedef Miscibility::NoMixing NoMix;
-
                     for (typename RMap::RegionId
                              r = 0, nr = reg.numRegions();
                          r < nr; ++r)
diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp
index 6b1885ea..70a6248d 100644
--- a/opm/core/simulator/initState_impl.hpp
+++ b/opm/core/simulator/initState_impl.hpp
@@ -184,7 +184,7 @@ namespace Opm
                 double A[4] = { 0.0 };
                 props_.matrix(1, &pressure, surfvol[phase], cells, A, 0);
                 double rho[2] = { 0.0 };
-                props_.density(1, A, rho);
+                props_.density(1, A, cells, rho);
                 return rho[phase];
             }
         };
diff --git a/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp
index 4370a8ef..82f8ea80 100644
--- a/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp
+++ b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp
@@ -418,7 +418,7 @@ namespace Opm
         assert(np == 2);
         const int dim = grid_.dimensions;
         density_.resize(nc*np);
-        props_.density(grid_.number_of_cells, &A_[0], &density_[0]);
+        props_.density(grid_.number_of_cells, &A_[0], /*cellIndices=*/NULL, &density_[0]);
         std::fill(gravflux_.begin(), gravflux_.end(), 0.0);
         for (int f = 0; f < nf; ++f) {
             const int* c = &grid_.face_cells[2*f];
diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp
index 2454eb73..8a0ebac8 100644
--- a/opm/core/utility/miscUtilitiesBlackoil.cpp
+++ b/opm/core/utility/miscUtilitiesBlackoil.cpp
@@ -285,7 +285,8 @@ namespace Opm
 
     /// @brief Computes saturation from surface volume
     void computeSaturation(const BlackoilPropertiesInterface& props,
-                          BlackoilState& state){
+                           BlackoilState& state)
+    {
 
         const int np = props.numPhases();
         const int nc = props.numCells();
diff --git a/tests/capillary.DATA b/tests/capillary.DATA
index 89c93899..64429cf0 100644
--- a/tests/capillary.DATA
+++ b/tests/capillary.DATA
@@ -19,6 +19,10 @@ PVDG
 200 0.005 0.2
 /
 
+PVTW
+1.0 1.0 4.0E-5 0.96 0.0
+/
+
 SWOF
 0.2 0 1 0.4
 1   1 0 0.1
diff --git a/tests/deadfluids.DATA b/tests/deadfluids.DATA
index 89232054..bc88ccb5 100644
--- a/tests/deadfluids.DATA
+++ b/tests/deadfluids.DATA
@@ -19,6 +19,10 @@ PVDG
 200 0.02 0.2
 /
 
+PVTW
+1.0 1.0 4.0E-5 0.96 0.0
+/
+
 SWOF
 0 0 1 0
 1 1 0 0
diff --git a/tests/liveoil.DATA b/tests/liveoil.DATA
index 48fcafe6..b3c3e271 100644
--- a/tests/liveoil.DATA
+++ b/tests/liveoil.DATA
@@ -9,6 +9,19 @@ WATER
 
 FIELD
 
+-- tests for the PVT functions need a grid because the OPM-API for the
+-- PVT functions assumes the presence of compressed grid cells,
+-- i.e. we need to be able to map from compressed to logical cartesian
+-- cell indices to find out the PVT region of a cell
+DIMENS
+3 3 3 /
+DXV
+1 2 3 /
+DYV
+1 2 3 /
+DZV
+1 2 3 /
+
 PVTO
 --     Rs       Pbub       Bo        Vo
       .0         14.7     1.0000     1.20  /
diff --git a/tests/test_blackoilfluid.cpp b/tests/test_blackoilfluid.cpp
index 5e2107f4..6f5cc3ca 100644
--- a/tests/test_blackoilfluid.cpp
+++ b/tests/test_blackoilfluid.cpp
@@ -1,15 +1,17 @@
 #include <config.h>
 
 #include <opm/core/grid/GridManager.hpp>
-#include <opm/core/props/pvt/SinglePvtConstCompr.hpp>
-#include <opm/core/props/pvt/SinglePvtDead.hpp>
-#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
-#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
-#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
+#include <opm/core/props/pvt/PvtConstCompr.hpp>
+#include <opm/core/props/pvt/PvtDead.hpp>
+#include <opm/core/props/pvt/PvtDeadSpline.hpp>
+#include <opm/core/props/pvt/PvtLiveOil.hpp>
+#include <opm/core/props/pvt/PvtLiveGas.hpp>
 #include <opm/core/props/phaseUsageFromDeck.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/core/props/pvt/BlackoilPvtProperties.hpp>
 #include <opm/core/utility/Units.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/grid.h>
 
 #include <opm/parser/eclipse/Utility/PvtoTable.hpp>
 #include <opm/parser/eclipse/Utility/PvtgTable.hpp>
@@ -36,63 +38,68 @@ using namespace Opm;
 using namespace std;
 
 
-std::vector<std::shared_ptr<SinglePvtInterface> > getProps(Opm::DeckConstPtr deck,PhaseUsage phase_usage_){
-
+std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, PhaseUsage phase_usage_){
+    Opm::GridManager grid(deck);
 
     enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
     int samples = 0;
 
-
-    std::vector<std::shared_ptr<SinglePvtInterface> > props_;
+    std::vector<std::shared_ptr<PvtInterface> > props_;
     // Set the properties.
     props_.resize(phase_usage_.num_phases);
 
     // Water PVT
     if (phase_usage_.phase_used[Aqua]) {
         if (deck->hasKeyword("PVTW")) {
-            Opm::PvtwTable pvtwTable(deck->getKeyword("PVTW"));
-            props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable));
+            std::shared_ptr<PvtConstCompr> pvtw(new PvtConstCompr);
+            pvtw->initFromWater(deck->getKeyword("PVTW"));
+
+            props_[phase_usage_.phase_pos[Aqua]] = pvtw;
         } else {
             // Eclipse 100 default.
-            props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
+            props_[phase_usage_.phase_pos[Aqua]].reset(new PvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
         }
     }
 
     // Oil PVT
     if (phase_usage_.phase_used[Liquid]) {
         if (deck->hasKeyword("PVDO")) {
-            Opm::PvdoTable pvdoTable(deck->getKeyword("PVDO"));
+            Opm::DeckKeywordConstPtr pvdoKeyword(deck->getKeyword("PVDO"));
             if (samples > 0) {
-                props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples));
+                std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
+                splinePvt->initFromOil(pvdoKeyword, samples);
+                props_[phase_usage_.phase_pos[Liquid]] = splinePvt;
             } else {
-                props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(pvdoTable));
+                std::shared_ptr<PvtDead> deadPvt(new PvtDead);
+                deadPvt->initFromOil(pvdoKeyword);
+                props_[phase_usage_.phase_pos[Liquid]] = deadPvt;
             }
         } else if (deck->hasKeyword("PVTO")) {
-            Opm::PvtoTable pvtoTable(deck->getKeyword("PVTO"), /*tableIdx=*/0);
-
-            props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable));
+            props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(deck->getKeyword("PVTO")));
         } else if (deck->hasKeyword("PVCDO")) {
-            Opm::PvcdoTable pvcdoTable(deck->getKeyword("PVCDO"), /*tableIdx=*/0);
+            std::shared_ptr<PvtConstCompr> pvcdo(new PvtConstCompr);
+            pvcdo->initFromOil(deck->getKeyword("PVCDO"));
 
-            props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable));
+            props_[phase_usage_.phase_pos[Liquid]] = pvcdo;
         } else {
-            OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n");
+            OPM_THROW(std::runtime_error, "Input is missing PVDO, PVCDO or PVTO\n");
         }
     }
     // Gas PVT
     if (phase_usage_.phase_used[Vapour]) {
         if (deck->hasKeyword("PVDG")) {
-            Opm::PvdgTable pvdgTable(deck->getKeyword("PVDG"));
-
+            Opm::DeckKeywordConstPtr pvdgKeyword(deck->getKeyword("PVDG"));
             if (samples > 0) {
-                props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples));
+                std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
+                splinePvt->initFromGas(pvdgKeyword, samples);
+                props_[phase_usage_.phase_pos[Vapour]] = splinePvt;
             } else {
-                props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(pvdgTable));
+                std::shared_ptr<PvtDead> deadPvt(new PvtDead);
+                deadPvt->initFromGas(pvdgKeyword);
+                props_[phase_usage_.phase_pos[Vapour]] = deadPvt;
             }
         } else if (deck->hasKeyword("PVTG")) {
-            Opm::PvtgTable pvtgTable(deck->getKeyword("PVTG"), /*tableIdx=*/0);
-
-            props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable));
+            props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(deck->getKeyword("PVTG")));
         } else {
             OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
         }
@@ -101,8 +108,8 @@ std::vector<std::shared_ptr<SinglePvtInterface> > getProps(Opm::DeckConstPtr dec
     return props_;
 }
 
-void testmu(const double reltol, int n, int np, std::vector<double> p, std::vector<double> r,std::vector<double> z,
-            std::vector<std::shared_ptr<SinglePvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
+void testmu(const double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<double> r,std::vector<double> z,
+            std::vector<std::shared_ptr<PvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
 {
     std::vector<double> mu(n);
     std::vector<double> dmudp(n);
@@ -115,8 +122,8 @@ void testmu(const double reltol, int n, int np, std::vector<double> p, std::vect
 
     // test mu
     for (int phase = 0; phase < np; ++phase) {
-        props_[phase]->mu(n, &p[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]);
-        props_[phase]->mu(n, &p[0], &z[0], &mu[0]);
+        props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]);
+        props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &z[0], &mu[0]);
         dmudp_diff = (mu_new[1]-mu_new[0])/(p[1]-p[0]);
         dmudr_diff = (mu_new[2]-mu_new[0])/(r[2]-r[0]);
         dmudp_diff_u = (mu_new[4]-mu_new[3])/(p[4]-p[3]);
@@ -137,8 +144,8 @@ void testmu(const double reltol, int n, int np, std::vector<double> p, std::vect
     }
 }
 
-void testb(const double reltol, int n, int np, std::vector<double> p, std::vector<double> r,std::vector<double> z,
-            std::vector<std::shared_ptr<SinglePvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
+void testb(const double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<double> r,std::vector<double> z,
+            std::vector<std::shared_ptr<PvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
 {
     // test b
     std::vector<double> b(n);
@@ -154,8 +161,8 @@ void testb(const double reltol, int n, int np, std::vector<double> p, std::vecto
     double dbdr_diff_u;
 
     for (int phase = 0; phase < np; ++phase) {
-        props_[phase]->b(n, &p[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]);
-        props_[phase]->dBdp(n, &p[0], &z[0], &B[0], &dBdp[0]);
+        props_[phase]->b(n, &pvtTableIdx[0], &p[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]);
+        props_[phase]->dBdp(n, &pvtTableIdx[0], &p[0], &z[0], &B[0], &dBdp[0]);
         dbdp_diff = (b[1]-b[0])/(p[1]-p[0]);
         dbdr_diff = (b[2]-b[0])/(r[2]-r[0]);
         dbdp_diff_u = (b[4]-b[3])/(p[4]-p[3]);
@@ -180,7 +187,7 @@ void testb(const double reltol, int n, int np, std::vector<double> p, std::vecto
     }
 }
 
-void testrsSat(double reltol, int n, int np, std::vector<double> p, std::vector<std::shared_ptr<SinglePvtInterface> > props_){
+void testrsSat(double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<std::shared_ptr<PvtInterface> > props_){
     // test bublepoint pressure
     std::vector<double> rs(n);
     std::vector<double> drsdp(n);
@@ -188,7 +195,7 @@ void testrsSat(double reltol, int n, int np, std::vector<double> p, std::vector<
     double drsdp_diff_u;
 
     for (int phase = 0; phase < np; ++phase) {
-        props_[phase] ->rsSat(n, &p[0], &rs[0], &drsdp[0]);
+        props_[phase] ->rsSat(n, &pvtTableIdx[0], &p[0], &rs[0], &drsdp[0]);
 
         drsdp_diff = (rs[1]-rs[0])/(p[1]-p[0]);
         drsdp_diff_u = (rs[4]-rs[3])/(p[4]-p[3]);
@@ -202,7 +209,7 @@ void testrsSat(double reltol, int n, int np, std::vector<double> p, std::vector<
     }
 }
 
-void testrvSat(double reltol, int n, int np, std::vector<double> p, std::vector<std::shared_ptr<SinglePvtInterface> > props_){
+void testrvSat(double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<std::shared_ptr<PvtInterface> > props_){
     // test rv saturated
     std::vector<double> rv(n);
     std::vector<double> drvdp(n);
@@ -210,7 +217,7 @@ void testrvSat(double reltol, int n, int np, std::vector<double> p, std::vector<
     double drvdp_diff_u;
 
     for (int phase = 0; phase < np; ++phase) {
-        props_[phase] ->rvSat(n, &p[0], &rv[0], &drvdp[0]);
+        props_[phase] ->rvSat(n, &pvtTableIdx[0], &p[0], &rv[0], &drvdp[0]);
 
         drvdp_diff = (rv[1]-rv[0])/(p[1]-p[0]);
         drvdp_diff_u = (rv[4]-rv[3])/(p[4]-p[3]);
@@ -234,15 +241,18 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
     Opm::ParserPtr parser(new Opm::Parser());
     Opm::DeckConstPtr deck(parser->parseFile(filename));
 
+    
+
     // setup pvt interface
     PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
-    std::vector<std::shared_ptr<SinglePvtInterface> > props_ = getProps(deck,phase_usage_);
+    std::vector<std::shared_ptr<PvtInterface> > props_ = getProps(deck,phase_usage_);
 
 
     // setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
     // approximation of the derivatives.
     const int n = 6;
     const int np = phase_usage_.num_phases;
+    std::vector<int> pvtRegionIdx(n, 0);
 
     // the tolerance for acceptable difference in values
     const double reltol = 1e-9;
@@ -287,13 +297,13 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
 
     }
 
-    testmu(reltol, n, np, p, r,z, props_, condition);
+    testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition);
 
-    testb(reltol,n,np,p,r,z,props_,condition);
+    testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition);
 
-    testrsSat(reltol,n,np,p,props_);
+    testrsSat(reltol,n,np,pvtRegionIdx,p,props_);
 
-    testrvSat(reltol,n,np,p,props_);
+    testrvSat(reltol,n,np,pvtRegionIdx,p,props_);
 
 
 }
@@ -311,13 +321,14 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
 
     // setup pvt interface
     PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
-    std::vector<std::shared_ptr<SinglePvtInterface> > props_ = getProps(deck,phase_usage_);
+    std::vector<std::shared_ptr<PvtInterface> > props_ = getProps(deck,phase_usage_);
 
 
     // setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
     // approximation of the derivatives.
     const int n = 6;
     const int np = phase_usage_.num_phases;
+    std::vector<int> pvtRegionIdx(n, 0);
 
     // the tolerance for acceptable difference in values
     const double reltol = 1e-9;
@@ -362,12 +373,12 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
 
     }
 
-    testmu(reltol, n, np, p, r,z, props_, condition);
+    testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition);
 
-    testb(reltol,n,np,p,r,z,props_,condition);
+    testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition);
 
-    testrsSat(reltol,n,np,p,props_);
+    testrsSat(reltol,n,np,pvtRegionIdx,p,props_);
 
-    testrvSat(reltol,n,np,p,props_);
+    testrvSat(reltol,n,np,pvtRegionIdx,p,props_);
 
 }
diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp
index e94f5d26..ffe3f736 100644
--- a/tests/test_equil.cpp
+++ b/tests/test_equil.cpp
@@ -428,13 +428,13 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary)
     // the true answer or something else.
     const double reltol = 1.0e-6;
     BOOST_CHECK_CLOSE(pressures[0][first] , 1.469769063e7   , reltol);
-    BOOST_CHECK_CLOSE(pressures[0][last ] , 1.545e7   , reltol);
-    BOOST_CHECK_CLOSE(pressures[1][last] , 1.546e7   , reltol);
+    BOOST_CHECK_CLOSE(pressures[0][last ] , 15452880.328284413   , reltol);
+    BOOST_CHECK_CLOSE(pressures[1][last] , 15462880.328284413   , reltol);
 
     const auto& sats = comp.saturation();
     const std::vector<double> s[3]{
-        { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.425893333333, 0.774026666666, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
-        { 0, 0, 0, 0.00736, 0.792746666666, 0.8, 0.8, 0.8, 0.8, 0.574106666666, 0.225973333333, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+        { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.42192000000000002, 0.77802666666666664, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
+        { 0, 0, 0, 0.00736, 0.792746666666, 0.8, 0.8, 0.8, 0.8, 0.57807999999999993, 0.22197333333333336, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
         { 0.8, 0.8, 0.8, 0.79264, 0.007253333333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
     };
     for (int phase = 0; phase < 3; ++phase) {
diff --git a/tests/wetgas.DATA b/tests/wetgas.DATA
index 3c094dc1..1a98a975 100644
--- a/tests/wetgas.DATA
+++ b/tests/wetgas.DATA
@@ -9,6 +9,19 @@ GAS
 
 FIELD
 
+-- tests for the PVT functions need a grid because the OPM-API for the
+-- PVT functions assumes the presence of compressed grid cells,
+-- i.e. we need to be able to map from compressed to logical cartesian
+-- cell indices to find out the PVT region of a cell
+DIMENS
+3 3 3 /
+DXV
+1 2 3 /
+DYV
+1 2 3 /
+DZV
+1 2 3 /
+
 -- PVT PROPERTIES OF DRY GAS (NO VAPOURISED OIL)
 -- FROM SPE3 Blackoil Kleppe
 --