diff --git a/opm/core/props/pvt/PvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp
index d4118aaf..db9cdb39 100644
--- a/opm/core/props/pvt/PvtDead.cpp
+++ b/opm/core/props/pvt/PvtDead.cpp
@@ -41,6 +41,7 @@ namespace Opm
         // resize the attributes of the object
         b_.resize(numRegions);
         viscosity_.resize(numRegions);
+        inverseBmu_.resize(numRegions);
 
         for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
             const Opm::PvdoTable& pvdoTable = pvdoTables[regionIdx];
@@ -51,12 +52,19 @@ namespace Opm
             const std::vector<double>& visc = pvdoTable.getViscosityColumn();
 
             const int sz = b.size();
-            std::vector<double> bInv(sz);
+            std::vector<double> inverseB(sz);
             for (int i = 0; i < sz; ++i) {
-                bInv[i] = 1.0 / b[i];
+                inverseB[i] = 1.0 / b[i];
             }
-            b_[regionIdx] = NonuniformTableLinear<double>(press, bInv);
+
+            std::vector<double> inverseBmu(sz);
+            for (int i = 0; i < sz; ++i) {
+                inverseBmu[i] = 1.0 / (b[i] * visc[i]);
+            }
+
+            b_[regionIdx] = NonuniformTableLinear<double>(press, inverseB);
             viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
+            inverseBmu_[regionIdx] = NonuniformTableLinear<double>(press, inverseBmu);
         }
     }
 
@@ -68,6 +76,7 @@ namespace Opm
         // resize the attributes of the object
         b_.resize(numRegions);
         viscosity_.resize(numRegions);
+        inverseBmu_.resize(numRegions);
 
         for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
             const Opm::PvdgTable& pvdgTable = pvdgTables[regionIdx];
@@ -78,12 +87,19 @@ namespace Opm
             const std::vector<double>& visc = pvdgTable.getViscosityColumn();
 
             const int sz = b.size();
-            std::vector<double> bInv(sz);
+            std::vector<double> inverseB(sz);
             for (int i = 0; i < sz; ++i) {
-                bInv[i] = 1.0 / b[i];
+                inverseB[i] = 1.0 / b[i];
             }
-            b_[regionIdx] = NonuniformTableLinear<double>(press, bInv);
+
+            std::vector<double> inverseBmu(sz);
+            for (int i = 0; i < sz; ++i) {
+                inverseBmu[i] = 1.0 / (b[i] * visc[i]);
+            }
+
+            b_[regionIdx] = NonuniformTableLinear<double>(press, inverseB);
             viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
+            inverseBmu_[regionIdx] = NonuniformTableLinear<double>(press, inverseBmu);
         }
     }
 
@@ -103,7 +119,9 @@ namespace Opm
 // #pragma omp parallel for
         for (int i = 0; i < n; ++i) {
             int regionIdx = getTableIndex_(pvtTableIdx, i);
-            output_mu[i] = viscosity_[regionIdx](p[i]);
+            double tempInvB = b_[regionIdx](p[i]);
+            double tempInvBmu = inverseBmu_[regionIdx](p[i]);
+            output_mu[i] = tempInvB / tempInvBmu;
         }
     }
 
@@ -118,8 +136,11 @@ namespace Opm
     // #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]);
+                double tempInvB = b_[regionIdx](p[i]);
+                double tempInvBmu = inverseBmu_[regionIdx](p[i]);
+                output_mu[i] = tempInvB / tempInvBmu;
+                output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i])
+                                 - tempInvB * inverseBmu_[regionIdx].derivative(p[i])) / (tempInvBmu * tempInvBmu);
             }
             std::fill(output_dmudr, output_dmudr + n, 0.0);
 
@@ -137,8 +158,12 @@ namespace Opm
     // #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]);
+                double tempInvB = b_[regionIdx](p[i]);
+                double tempInvBmu = inverseBmu_[regionIdx](p[i]);
+                output_mu[i] = tempInvB / tempInvBmu;
+                output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i])
+                                 - tempInvB * inverseBmu_[regionIdx].derivative(p[i]))
+                                 / (tempInvBmu * tempInvBmu);
             }
             std::fill(output_dmudr, output_dmudr + n, 0.0);
 
diff --git a/opm/core/props/pvt/PvtDead.hpp b/opm/core/props/pvt/PvtDead.hpp
index d83d4f69..ccc9acac 100644
--- a/opm/core/props/pvt/PvtDead.hpp
+++ b/opm/core/props/pvt/PvtDead.hpp
@@ -155,6 +155,7 @@ namespace Opm
         // table per PVT region.
         std::vector<NonuniformTableLinear<double> > b_;
         std::vector<NonuniformTableLinear<double> > viscosity_;
+        std::vector<NonuniformTableLinear<double> > inverseBmu_;
     };
 }
 
diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp
index cb1cb829..63446b3c 100644
--- a/opm/core/props/pvt/PvtLiveGas.cpp
+++ b/opm/core/props/pvt/PvtLiveGas.cpp
@@ -56,28 +56,40 @@ namespace Opm
             const Opm::PvtgTable& pvtgTable = pvtgTables[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();
+            // saturated_gas_table_[pvtTableIdx].resize(4);
+            // Adding one extra line to PVTG to store 1./(Bg*mu_g)
+            saturated_gas_table_[pvtTableIdx].resize(5);
+            saturated_gas_table_[pvtTableIdx][0] = pvtgTable.getOuterTable()->getPressureColumn(); // Pressure
+            saturated_gas_table_[pvtTableIdx][1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); // Bg
+            saturated_gas_table_[pvtTableIdx][2] = pvtgTable.getOuterTable()->getGasViscosityColumn(); // mu_g
+            // The number of the columns
+            int nRows = saturated_gas_table_[pvtTableIdx][2].size();
+            saturated_gas_table_[pvtTableIdx][3].resize(nRows); // allocate memory for 1/(Bg*mu_g)
+            saturated_gas_table_[pvtTableIdx][4] = pvtgTable.getOuterTable()->getOilSolubilityColumn(); // Rv
 
             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();
+                // undersat_gas_tables_[pvtTableIdx][i].resize(3);
+                undersat_gas_tables_[pvtTableIdx][i].resize(4);
+                undersat_gas_tables_[pvtTableIdx][i][0] = undersatTable.getOilSolubilityColumn(); // Rv
+                undersat_gas_tables_[pvtTableIdx][i][1] = undersatTable.getGasFormationFactorColumn(); // Bg
+                undersat_gas_tables_[pvtTableIdx][i][2] = undersatTable.getGasViscosityColumn();  // mu_g
+                int nRows = undersat_gas_tables_[pvtTableIdx][i][2].size();
+                undersat_gas_tables_[pvtTableIdx][i][3].resize(nRows); // allocate memory for 1/(Bg*mu_g)
             }
 
             // Bg -> 1/Bg
             for (int i=0; i<sz; ++i) {
-                saturated_gas_table_[pvtTableIdx][1][i] = 1.0/saturated_gas_table_[pvtTableIdx][1][i];
+                saturated_gas_table_[pvtTableIdx][3][i] = 1.0 / (saturated_gas_table_[pvtTableIdx][1][i]
+                                                               * saturated_gas_table_[pvtTableIdx][2][i]); // 1/(Bg*mu_g)
+                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];
+                    undersat_gas_tables_[pvtTableIdx][i][3][j] = 1.0 / (undersat_gas_tables_[pvtTableIdx][i][1][j]
+                                                                      * undersat_gas_tables_[pvtTableIdx][i][2][j]); // 1/(Bg*mu_g)
+                    undersat_gas_tables_[pvtTableIdx][i][1][j] = 1.0 / undersat_gas_tables_[pvtTableIdx][i][1][j];
                 }
             }
         }
@@ -97,7 +109,10 @@ namespace Opm
     {
 // #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);
+            double inverseB = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 1, false);
+            double inverseBMu = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 3, false);
+
+            output_mu[i] = inverseB / inverseBMu;
         }
     }
 
@@ -126,9 +141,24 @@ namespace Opm
         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);
+
+            double inverseBMu = miscible_gas(p[i], r[i], cnd, tableIdx, 3, 0);
+            double inverseB = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 0);
+
+            output_mu[i] = inverseB / inverseBMu;
+
+            double dinverseBmudp = miscible_gas(p[i], r[i], cnd, tableIdx, 3, 1);
+            double dinverseBdp = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 1);
+
+            output_dmudp[i] = (inverseBMu * dinverseBdp - inverseB * dinverseBmudp)
+                              / (inverseBMu * inverseBMu);
+
+            double dinverseBmudr = miscible_gas(p[i], r[i], cnd, tableIdx, 3, 2);
+            double dinverseBdr = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 2);
+
+            output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr)
+                              / (inverseBMu * inverseBMu);
+
         }
 
     }
@@ -209,9 +239,9 @@ namespace Opm
         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]);
+                    saturated_gas_table_[pvtTableIdx][4],p[i]);
             output_drvSatdp[i] = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0],
-                    saturated_gas_table_[pvtTableIdx][3],p[i]);
+                    saturated_gas_table_[pvtTableIdx][4],p[i]);
 
         }
     }
@@ -287,7 +317,7 @@ namespace Opm
             return 0.0;
         }
         double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
-                                             saturated_gas_table_[pvtTableIdx][3], press);
+                                             saturated_gas_table_[pvtTableIdx][4], press);
         double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
         if (satR < maxR ) {
             // Saturated case
@@ -308,13 +338,13 @@ namespace Opm
             return;
         }
         double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
-                                             saturated_gas_table_[pvtTableIdx][3], press);
+                                             saturated_gas_table_[pvtTableIdx][4], 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],
+                                            saturated_gas_table_[pvtTableIdx][4],
                                             press);
         } else {
             // Undersaturated case
@@ -336,7 +366,7 @@ namespace Opm
 
         int section;
         double Rval = linearInterpolation(saturatedGasTable[0],
-                                                saturatedGasTable[3], press,
+                                                saturatedGasTable[4], press,
                                                 section);
         double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
         if (deriv) {
diff --git a/opm/core/props/pvt/PvtLiveOil.cpp b/opm/core/props/pvt/PvtLiveOil.cpp
index a4ac13d5..012a1fdb 100644
--- a/opm/core/props/pvt/PvtLiveOil.cpp
+++ b/opm/core/props/pvt/PvtLiveOil.cpp
@@ -46,36 +46,47 @@ namespace Opm
             const auto saturatedPvto = pvtoTable.getOuterTable();
 
             // OIL, PVTO
-            saturated_oil_table_[pvtTableIdx].resize(4);
+            // saturated_oil_table_[pvtTableIdx].resize(4);
+            // adding a extra colummn to the PVTO to store 1/(B*mu)
+            saturated_oil_table_[pvtTableIdx].resize(5);
             const int sz = saturatedPvto->numRows();
-            for (int k=0; k<4; ++k) {
+            // for (int k=0; k<4; ++k) {
+            for (int k=0; k<5; ++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
+                saturated_oil_table_[pvtTableIdx][3][i] = 1.0 / (saturatedPvto->getOilFormationFactorColumn()[i]
+                                                               * saturatedPvto->getOilViscosityColumn()[i]); // 1/(Bo*mu_o)
+                saturated_oil_table_[pvtTableIdx][4][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);
+                // undersat_oil_tables_[pvtTableIdx][i].resize(3);
+                // adding a extra colummn to the PVTO to store 1/(B*mu)
+                undersat_oil_tables_[pvtTableIdx][i].resize(4);
                 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);
+                undersat_oil_tables_[pvtTableIdx][i][3].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
+                    undersat_oil_tables_[pvtTableIdx][i][3][j] = 1.0 / (undersaturatedPvto->getOilFormationFactorColumn()[j] *
+                                                                        undersaturatedPvto->getOilViscosityColumn()[j]);  // 1/(Bo*mu_o)
                 }
             }
 
             // Complete undersaturated tables by extrapolating from existing data
             // as is done in Eclipse and Mrst
+            // TODO: check if the following formulations applying to 1/(Bo*mu_o)
             int iNext = -1;
             for (int i=0; i<sz; ++i) {
                 // Skip records already containing undersaturated data
@@ -91,6 +102,7 @@ namespace Opm
                     if (iNext == sz) OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
                 }
                 // Add undersaturated data to current record while maintaining compressibility and viscosibility
+                // TODO: How to add 1/(B*mu) in this way?
                 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];
@@ -104,6 +116,11 @@ namespace Opm
                         / (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);
+
+                    // A try to expolate the 1/BMu with the expolated mu and B
+                    double inverseBMu = 1.0 / (B*mu);
+                    undersat_oil_tables_[pvtTableIdx][i][3].push_back(inverseBMu);
+
                 }
             }
         }
@@ -126,7 +143,10 @@ namespace Opm
         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);
+            double inverseB = miscible_oil(p[i], z + num_phases_*i, tableIdx, 1, false);
+            double inverseBMu = miscible_oil(p[i], z + num_phases_*i, tableIdx, 3, false);
+
+            output_mu[i] = inverseB / inverseBMu; 
         }
     }
 
@@ -143,9 +163,23 @@ namespace Opm
                 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);
+                    double inverseBMu = miscible_oil(p[i], r[i], tableIdx, 3, 0);
+                    double inverseB = miscible_oil(p[i], r[i], tableIdx, 1, 0);
+
+                    output_mu[i] = inverseB / inverseBMu;
+
+                    double dinverseBmudp = miscible_oil(p[i], r[i], tableIdx, 3, 1);
+                    double dinverseBdp = miscible_oil(p[i], r[i], tableIdx, 1, 1);
+
+                    output_dmudp[i] = (inverseBMu * dinverseBdp - inverseB * dinverseBmudp)
+                                      / (inverseBMu * inverseBMu);
+
+
+                    double dinverseBmudr = miscible_oil(p[i], r[i], tableIdx, 3, 2);
+                    double dinverseBdr = miscible_oil(p[i], r[i], tableIdx, 1, 2);
+
+                    output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr)
+                                      / (inverseBMu * inverseBMu);
 
                 }
     }
@@ -165,9 +199,23 @@ namespace Opm
                     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);
+                    double inverseBMu = miscible_oil(p[i], r[i], cnd, tableIdx, 3, 0);
+                    double inverseB = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 0);
+
+                    output_mu[i] = inverseB / inverseBMu;
+
+                    double dinverseBmudp = miscible_oil(p[i], r[i], cnd, tableIdx, 3, 1);
+                    double dinverseBdp = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 1);
+
+                    output_dmudp[i] = (inverseBMu * dinverseBdp - inverseB * dinverseBmudp)
+                                      / (inverseBMu * inverseBMu);
+
+
+                    double dinverseBmudr = miscible_oil(p[i], r[i], cnd, tableIdx, 3, 2);
+                    double dinverseBdr = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 2);
+
+                    output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr)
+                                      / (inverseBMu * inverseBMu);
 
                 }
     }
@@ -258,9 +306,9 @@ namespace Opm
         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]);
+                    saturated_oil_table_[tableIdx][4],p[i]);
             output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[tableIdx][0],
-                    saturated_oil_table_[tableIdx][3],p[i]);
+                    saturated_oil_table_[tableIdx][4],p[i]);
 
         }
     }
@@ -334,7 +382,7 @@ namespace Opm
             return 0.0;
         }
         double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
-                                             saturated_oil_table_[pvtTableIdx][3], press);
+                                             saturated_oil_table_[pvtTableIdx][4], press);
         double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
         if (Rval < maxR ) {  // Saturated case
             return Rval;
@@ -353,12 +401,12 @@ namespace Opm
             return;
         }
         Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
-                                   saturated_oil_table_[pvtTableIdx][3], press);
+                                   saturated_oil_table_[pvtTableIdx][4], 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],
+                                                    saturated_oil_table_[pvtTableIdx][4],
                                                     press);
         } else {
             // Undersaturated case
@@ -368,6 +416,7 @@ namespace Opm
     }
 
 
+    // TODO: Check if this function need to be adapted to the 1/(B*mu) interpolation.
     double PvtLiveOil::miscible_oil(const double press,
                                     const double* surfvol,
                                     const int pvtTableIdx,
@@ -376,7 +425,7 @@ namespace Opm
     {
         int section;
         double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
-                                          saturated_oil_table_[pvtTableIdx][3],
+                                          saturated_oil_table_[pvtTableIdx][4],
                                           press, section);
         double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
         if (deriv) {
@@ -385,9 +434,9 @@ namespace Opm
                                                      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]);
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], maxR);
+                double w = (maxR - saturated_oil_table_[pvtTableIdx][4][is]) /
+                    (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
                 assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
                 assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
                 double val1 =
@@ -408,9 +457,9 @@ namespace Opm
                                            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]);
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], maxR);
+                double w = (maxR - saturated_oil_table_[pvtTableIdx][4][is]) /
+                    (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
                 assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
                 assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
                 double val1 =
@@ -435,7 +484,7 @@ namespace Opm
     {
         int section;
         double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
-                                          saturated_oil_table_[pvtTableIdx][3],
+                                          saturated_oil_table_[pvtTableIdx][4],
                                           press, section);
         // derivative with respect to frist component (pressure)
         if (deriv == 1) {
@@ -444,9 +493,9 @@ namespace Opm
                                                      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]);
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
+                double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) /
+                    (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
                 assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
                 assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
                 double val1 =
@@ -465,7 +514,7 @@ namespace Opm
             if (Rval <= r ) {  // Saturated case
                 return 0;
             } else {  // Undersaturated case
-                int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
                 assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
                 assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
                 double val1 =
@@ -477,7 +526,7 @@ namespace Opm
                                               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]);
+                double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][4][is+1]-saturated_oil_table_[pvtTableIdx][4][is]);
                 return val;
             }
 
@@ -489,9 +538,9 @@ namespace Opm
                                            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]);
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
+                double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) /
+                    (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
                 assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
                 assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
                 double val1 =
@@ -524,9 +573,9 @@ namespace Opm
                                                      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]);
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
+                double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) /
+                    (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
                 assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
                 assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
                 double val1 =
@@ -545,7 +594,7 @@ namespace Opm
             if (isSat) {  // Saturated case
                 return 0;
             } else {  // Undersaturated case
-                int is = tableIndex(saturated_oil_table_[pvtTableIdx][3], r);
+                int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
                 assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
                 assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
                 double val1 =
@@ -557,35 +606,37 @@ namespace Opm
                                               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]);
+                double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][4][is+1]-saturated_oil_table_[pvtTableIdx][4][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;
+            }
+            // no derivative
+            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][4], r);
+                    double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) /
+                        (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][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/tests/test_blackoilfluid.cpp b/tests/test_blackoilfluid.cpp
index 699ebb9a..add7d495 100644
--- a/tests/test_blackoilfluid.cpp
+++ b/tests/test_blackoilfluid.cpp
@@ -247,8 +247,10 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
     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;
+    // the relative tolerance in percentage for acceptable difference in values
+    const double reltolper = 1e-9;
+    // the relative tolerance in percentage for acceptable difference in values for viscosity
+    const double reltolpermu = 1e-1;
 
     std::vector<double> p(n);
     std::vector<double> r(n);
@@ -290,13 +292,13 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
 
     }
 
-    testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition);
+    testmu(reltolpermu, n, np, pvtRegionIdx, p, r,z, props_, condition);
 
-    testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition);
+    testb(reltolper,n,np,pvtRegionIdx,p,r,z,props_,condition);
 
-    testrsSat(reltol,n,np,pvtRegionIdx,p,props_);
+    testrsSat(reltolper,n,np,pvtRegionIdx,p,props_);
 
-    testrvSat(reltol,n,np,pvtRegionIdx,p,props_);
+    testrvSat(reltolper,n,np,pvtRegionIdx,p,props_);
 
 
 }
@@ -324,8 +326,10 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
     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;
+    // the relative tolerance in percentage for acceptable difference in values
+    const double reltolper = 1e-9;
+    // the relative tolerance in percentage for acceptable difference in values for viscosity
+    const double reltolpermu = 1e-1;
 
     std::vector<double> p(n);
     std::vector<double> r(n);
@@ -367,12 +371,12 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
 
     }
 
-    testmu(reltol, n, np, pvtRegionIdx, p, r,z, props_, condition);
+    testmu(reltolpermu, n, np, pvtRegionIdx, p, r,z, props_, condition);
 
-    testb(reltol,n,np,pvtRegionIdx,p,r,z,props_,condition);
+    testb(reltolper,n,np,pvtRegionIdx,p,r,z,props_,condition);
 
-    testrsSat(reltol,n,np,pvtRegionIdx,p,props_);
+    testrsSat(reltolper,n,np,pvtRegionIdx,p,props_);
 
-    testrvSat(reltol,n,np,pvtRegionIdx,p,props_);
+    testrvSat(reltolper,n,np,pvtRegionIdx,p,props_);
 
 }