diff --git a/examples/compute_initial_state.cpp b/examples/compute_initial_state.cpp
index a88638f02..2eb27d1f4 100644
--- a/examples/compute_initial_state.cpp
+++ b/examples/compute_initial_state.cpp
@@ -36,6 +36,8 @@
 
 #include <boost/filesystem.hpp>
 
+#include <fstream>
+
 namespace
 {
     void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
@@ -88,9 +90,9 @@ try
     Opm::ParseContext parseContext;
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::DeckConstPtr deck = parser->parseFile(deck_filename , parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck, parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
     const double grav = param.getDefault("gravity", unit::gravity);
-    GridManager gm(deck);
+    GridManager gm(eclipseState->getInputGrid());
     const UnstructuredGrid& grid = *gm.c_grid();
     BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param);
     warnIfUnusedParams(param);
diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp
index fd2069206..aedeebefe 100644
--- a/examples/compute_tof_from_files.cpp
+++ b/examples/compute_tof_from_files.cpp
@@ -54,7 +54,7 @@
 #include <vector>
 #include <numeric>
 #include <iterator>
-
+#include <fstream>
 
 namespace
 {
diff --git a/examples/diagnose_relperm.cpp b/examples/diagnose_relperm.cpp
index 862fb96d7..99158109d 100644
--- a/examples/diagnose_relperm.cpp
+++ b/examples/diagnose_relperm.cpp
@@ -76,9 +76,9 @@ try
                               { ParseContext::INTERNAL_ERROR_UNINITIALIZED_THPRES, InputError::IGNORE}
                              });
     Opm::DeckConstPtr deck(parser->parseFile(eclipseFilename, parseContext));
-    eclState.reset(new EclipseState(deck, parseContext));
+    eclState.reset(new EclipseState(*deck, parseContext));
 
-    GridManager gm(deck);
+    GridManager gm(eclState->getInputGrid());
     const UnstructuredGrid& grid = *gm.c_grid();
     using boost::filesystem::path; 
     path fpath(eclipseFilename);
@@ -92,13 +92,13 @@ try
     std::string logFile = baseName + ".SATFUNCLOG";
     std::shared_ptr<EclipsePRTLog> prtLog = std::make_shared<EclipsePRTLog>(logFile, Log::DefaultMessageTypes);
     OpmLog::addBackend( "ECLIPSEPRTLOG" , prtLog );
-    Opm::time::StopWatch timer;
-    timer.start();
+    prtLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(true, false));
+    std::shared_ptr<StreamLog> streamLog = std::make_shared<EclipsePRTLog>(std::cout, Log::DefaultMessageTypes);
+    OpmLog::addBackend( "STREAMLOG" , streamLog );
+    streamLog->setMessageLimiter(std::make_shared<MessageLimiter>(10));
+    streamLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(true, true));
     RelpermDiagnostics diagnostic;
     diagnostic.diagnosis(eclState, deck, grid);
-    timer.stop();
-    double tt = timer.secsSinceStart();
-    std::cout << "relperm diagnostics: " << tt << " seconds." << std::endl;
 }
 catch (const std::exception &e) {
     std::cerr << "Program threw an exception: " << e.what() << "\n";
diff --git a/opm/core/linalg/LinearSolverFactory.cpp b/opm/core/linalg/LinearSolverFactory.cpp
index d2bd31204..f8eca83d1 100644
--- a/opm/core/linalg/LinearSolverFactory.cpp
+++ b/opm/core/linalg/LinearSolverFactory.cpp
@@ -61,8 +61,19 @@ namespace Opm
 
     LinearSolverFactory::LinearSolverFactory(const parameter::ParameterGroup& param)
     {
+#if HAVE_SUITESPARSE_UMFPACK_H
+        std::string default_solver = "umfpack";
+#elif HAVE_DUNE_ISTL
+        std::string default_solver = "istl";
+#elif HAVE_PETSC
+        std::string default_solver = "petsc";
+#else
+        std::string default_solver = "no_solver_available";
+        OPM_THROW(std::runtime_error, "No linear solver available, you must have UMFPACK , dune-istl or Petsc installed to use LinearSolverFactory.");
+#endif
+
         const std::string ls =
-            param.getDefault<std::string>("linsolver", "umfpack");
+            param.getDefault("linsolver", default_solver);
 
         if (ls == "umfpack") {
 #if HAVE_SUITESPARSE_UMFPACK_H
diff --git a/opm/core/linalg/ParallelIstlInformation.hpp b/opm/core/linalg/ParallelIstlInformation.hpp
index 97c191b39..e96f88e17 100644
--- a/opm/core/linalg/ParallelIstlInformation.hpp
+++ b/opm/core/linalg/ParallelIstlInformation.hpp
@@ -179,12 +179,24 @@ public:
         }
         return ownerMask_;
     }
+
+    /// \brief Get the owner Mask.
+    ///
+    /// \return A vector with entries 0, and 1. 0 marks an index that we cannot
+    ///         compute correct results for. 1 marks an index that this process
+    ///         is responsible for and computes correct results in parallel.
+    const std::vector<double>& getOwnerMask() const
+    {
+        return ownerMask_;
+    }
+
     /// \brief Compute one or more global reductions.
     ///
     /// This function can either be used with a container, an operator, and an initial value
     /// to compute a reduction. Or with tuples of them to compute multiple reductions with only
     /// one global communication.
     /// The possible functors needed can be constructed with Opm::Reduction::makeGlobalMaxFunctor(),
+    /// Opm::Reduction::makeLInfinityNormFunctor(),
     /// Opm::Reduction::makeGlobalMinFunctor(), and 
     /// Opm::Reduction::makeGlobalSumFunctor().
     /// \tparam type of the container or the tuple of  containers.
@@ -574,6 +586,45 @@ private:
             (std::pointer_to_binary_function<const T&,const T&,const T&>
              ((const T&(*)(const T&, const T&))std::max<T>));
     }
+
+    namespace detail
+    {
+        /// \brief Computes the maximum of the absolute values of two values.
+        template<typename T, typename Enable = void>
+        struct MaxAbsFunctor
+        {
+            using result_type = T;
+            result_type operator()(const T& t1,
+                                   const T& t2)
+            {
+                return std::max(std::abs(t1), std::abs(t2));
+            }
+        };
+
+        // Specialization for unsigned integers. They need their own
+        // version since abs(x) is ambiguous (as well as somewhat
+        // meaningless).
+        template<typename T>
+        struct MaxAbsFunctor<T, typename std::enable_if<std::is_unsigned<T>::value>::type>
+        {
+            using result_type = T;
+            result_type operator()(const T& t1,
+                                   const T& t2)
+            {
+                return std::max(t1, t2);
+            }
+        };
+    }
+
+    /// \brief Create a functor for computing a global L infinity norm
+    ///
+    /// To be used with ParallelISTLInformation::computeReduction.
+    template<class T>
+    MaskIDOperator<detail::MaxAbsFunctor<T> >
+    makeLInfinityNormFunctor()
+    {
+        return MaskIDOperator<detail::MaxAbsFunctor<T> >();
+    }
     /// \brief Create a functor for computing a global minimum.
     ///
     /// To be used with ParallelISTLInformation::computeReduction.
diff --git a/opm/core/props/BlackoilPhases.hpp b/opm/core/props/BlackoilPhases.hpp
index ebfd1f72a..981588ff7 100644
--- a/opm/core/props/BlackoilPhases.hpp
+++ b/opm/core/props/BlackoilPhases.hpp
@@ -60,6 +60,9 @@ namespace Opm
         void setFreeOil  () { insert(BlackoilPhases::Liquid); }
         void setFreeGas  () { insert(BlackoilPhases::Vapour); }
 
+        bool operator==(const PhasePresence& other) const { return present_ == other.present_; }
+        bool operator!=(const PhasePresence& other) const { return !this->operator==(other); }
+
     private:
         unsigned char present_;
 
diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp
index 97d0a6e65..b2b6b3780 100644
--- a/opm/core/props/BlackoilPropertiesFromDeck.cpp
+++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp
@@ -244,15 +244,15 @@ namespace Opm
         const auto& pu = phaseUsage();
         const int np = numPhases();
 
-        typedef Opm::LocalAd::Evaluation<double, /*size=*/1> LadEval;
+        typedef Opm::DenseAd::Evaluation<double, /*size=*/1> Eval;
 
-        LadEval pLad = 0.0;
-        LadEval TLad = 0.0;
-        LadEval RsLad = 0.0;
-        LadEval RvLad = 0.0;
-        LadEval muLad = 0.0;
+        Eval pEval = 0.0;
+        Eval TEval = 0.0;
+        Eval RsEval = 0.0;
+        Eval RvEval = 0.0;
+        Eval muEval = 0.0;
 
-        pLad.derivatives[0] = 1.0;
+        pEval.derivatives[0] = 1.0;
 
         R_.resize(n*np);
         this->compute_R_(n, p, T, z, cells, &R_[0]);
@@ -260,30 +260,30 @@ namespace Opm
         for (int i = 0; i < n; ++ i) {
             int cellIdx = cells[i];
             int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
-            pLad.value = p[i];
-            TLad.value = T[i];
+            pEval.value = p[i];
+            TEval.value = T[i];
 
             if (pu.phase_used[BlackoilPhases::Aqua]) {
-                muLad = waterPvt_.viscosity(pvtRegionIdx, TLad, pLad);
+                muEval = waterPvt_.viscosity(pvtRegionIdx, TEval, pEval);
                 int offset = pu.num_phases*cellIdx + pu.phase_pos[BlackoilPhases::Aqua];
-                mu[offset] = muLad.value;
-                dmudp[offset] = muLad.derivatives[0];
+                mu[offset] = muEval.value;
+                dmudp[offset] = muEval.derivatives[0];
             }
 
             if (pu.phase_used[BlackoilPhases::Liquid]) {
-                RsLad.value = R_[i*np + pu.phase_pos[BlackoilPhases::Liquid]];
-                muLad = oilPvt_.viscosity(pvtRegionIdx, TLad, pLad, RsLad);
+                RsEval.value = R_[i*np + pu.phase_pos[BlackoilPhases::Liquid]];
+                muEval = oilPvt_.viscosity(pvtRegionIdx, TEval, pEval, RsEval);
                 int offset = pu.num_phases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid];
-                mu[offset] = muLad.value;
-                dmudp[offset] = muLad.derivatives[0];
+                mu[offset] = muEval.value;
+                dmudp[offset] = muEval.derivatives[0];
             }
 
             if (pu.phase_used[BlackoilPhases::Vapour]) {
-                RvLad.value = R_[i*np + pu.phase_pos[BlackoilPhases::Vapour]];
-                muLad = gasPvt_.viscosity(pvtRegionIdx, TLad, pLad, RvLad);
+                RvEval.value = R_[i*np + pu.phase_pos[BlackoilPhases::Vapour]];
+                muEval = gasPvt_.viscosity(pvtRegionIdx, TEval, pEval, RvEval);
                 int offset = pu.num_phases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour];
-                mu[offset] = muLad.value;
-                dmudp[offset] = muLad.derivatives[0];
+                mu[offset] = muEval.value;
+                dmudp[offset] = muEval.derivatives[0];
             }
         }
     }
@@ -396,27 +396,27 @@ namespace Opm
     {
         const auto& pu = phaseUsage();
 
-        typedef double LadEval;
+        typedef double Eval;
 
-        LadEval pLad = 0.0;
-        LadEval TLad = 0.0;
-        LadEval RsLad = 0.0;
-        LadEval RvLad = 0.0;
+        Eval pEval = 0.0;
+        Eval TEval = 0.0;
+        Eval RsEval = 0.0;
+        Eval RvEval = 0.0;
 
         for (int i = 0; i < n; ++ i) {
             int cellIdx = cells[i];
             int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
-            pLad = p[i];
-            TLad = T[i];
+            pEval = p[i];
+            TEval = T[i];
 
             int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid];
             int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour];
             int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua];
 
             if (pu.phase_used[BlackoilPhases::Aqua]) {
-                LadEval BLad = 1.0/waterPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
+                Eval BEval = 1.0/waterPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval);
 
-                B[waterOffset] = BLad;
+                B[waterOffset] = BEval;
             }
 
             if (pu.phase_used[BlackoilPhases::Liquid]) {
@@ -424,18 +424,18 @@ namespace Opm
                 double maxRs = 0.0;
                 if (pu.phase_used[BlackoilPhases::Vapour]) {
                     currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset];
-                    maxRs = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad, pLad);
+                    maxRs = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TEval, pEval);
                 }
-                LadEval BLad;
+                Eval BEval;
                 if (currentRs >= maxRs) {
-                    BLad = 1.0/oilPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
+                    BEval = 1.0/oilPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval);
                 }
                 else {
-                    RsLad = currentRs;
-                    BLad = 1.0/oilPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RsLad);
+                    RsEval = currentRs;
+                    BEval = 1.0/oilPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval, RsEval);
                 }
 
-                B[oilOffset] = BLad;
+                B[oilOffset] = BEval;
             }
 
             if (pu.phase_used[BlackoilPhases::Vapour]) {
@@ -443,18 +443,18 @@ namespace Opm
                 double maxRv = 0.0;
                 if (pu.phase_used[BlackoilPhases::Liquid]) {
                     currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset];
-                    maxRv = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad, pLad);
+                    maxRv = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TEval, pEval);
                 }
-                LadEval BLad;
+                Eval BEval;
                 if (currentRv >= maxRv) {
-                    BLad = 1.0/gasPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
+                    BEval = 1.0/gasPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval);
                 }
                 else {
-                    RvLad = currentRv;
-                    BLad = 1.0/gasPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RvLad);
+                    RvEval = currentRv;
+                    BEval = 1.0/gasPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval, RvEval);
                 }
 
-                B[gasOffset] = BLad;
+                B[gasOffset] = BEval;
             }
         }
     }
@@ -469,30 +469,30 @@ namespace Opm
     {
         const auto& pu = phaseUsage();
 
-        typedef Opm::LocalAd::Evaluation<double, /*size=*/1> LadEval;
+        typedef Opm::DenseAd::Evaluation<double, /*size=*/1> Eval;
 
-        LadEval pLad = 0.0;
-        LadEval TLad = 0.0;
-        LadEval RsLad = 0.0;
-        LadEval RvLad = 0.0;
+        Eval pEval = 0.0;
+        Eval TEval = 0.0;
+        Eval RsEval = 0.0;
+        Eval RvEval = 0.0;
 
-        pLad.derivatives[0] = 1.0;
+        pEval.derivatives[0] = 1.0;
 
         for (int i = 0; i < n; ++ i) {
             int cellIdx = cells[i];
             int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
-            pLad.value = p[i];
-            TLad.value = T[i];
+            pEval.value = p[i];
+            TEval.value = T[i];
 
             int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid];
             int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour];
             int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua];
 
             if (pu.phase_used[BlackoilPhases::Aqua]) {
-                LadEval BLad = 1.0/waterPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
+                Eval BEval = 1.0/waterPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval);
 
-                B[waterOffset] = BLad.value;
-                dBdp[waterOffset] = BLad.derivatives[0];
+                B[waterOffset] = BEval.value;
+                dBdp[waterOffset] = BEval.derivatives[0];
             }
 
             if (pu.phase_used[BlackoilPhases::Liquid]) {
@@ -500,19 +500,19 @@ namespace Opm
                 double maxRs = 0.0;
                 if (pu.phase_used[BlackoilPhases::Vapour]) {
                     currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset];
-                    maxRs = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad.value, pLad.value);
+                    maxRs = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TEval.value, pEval.value);
                 }
-                LadEval BLad;
+                Eval BEval;
                 if (currentRs >= maxRs) {
-                    BLad = 1.0/oilPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
+                    BEval = 1.0/oilPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval);
                 }
                 else {
-                    RsLad.value = currentRs;
-                    BLad = 1.0/oilPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RsLad);
+                    RsEval.value = currentRs;
+                    BEval = 1.0/oilPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval, RsEval);
                 }
 
-                B[oilOffset] = BLad.value;
-                dBdp[oilOffset] = BLad.derivatives[0];
+                B[oilOffset] = BEval.value;
+                dBdp[oilOffset] = BEval.derivatives[0];
             }
 
             if (pu.phase_used[BlackoilPhases::Vapour]) {
@@ -520,19 +520,19 @@ namespace Opm
                 double maxRv = 0.0;
                 if (pu.phase_used[BlackoilPhases::Liquid]) {
                     currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset];
-                    maxRv = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad.value, pLad.value);
+                    maxRv = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TEval.value, pEval.value);
                 }
-                LadEval BLad;
+                Eval BEval;
                 if (currentRv >= maxRv) {
-                    BLad = 1.0/gasPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
+                    BEval = 1.0/gasPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval);
                 }
                 else {
-                    RvLad.value = currentRv;
-                    BLad = 1.0/gasPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RvLad);
+                    RvEval.value = currentRv;
+                    BEval = 1.0/gasPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval, RvEval);
                 }
 
-                B[gasOffset] = BLad.value;
-                dBdp[gasOffset] = BLad.derivatives[0];
+                B[gasOffset] = BEval.value;
+                dBdp[gasOffset] = BEval.derivatives[0];
             }
         }
     }
@@ -546,16 +546,16 @@ namespace Opm
     {
         const auto& pu = phaseUsage();
 
-        typedef double LadEval;
+        typedef double Eval;
 
-        LadEval pLad = 0.0;
-        LadEval TLad = 0.0;
+        Eval pEval = 0.0;
+        Eval TEval = 0.0;
 
         for (int i = 0; i < n; ++ i) {
             int cellIdx = cells[i];
             int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
-            pLad = p[i];
-            TLad = T[i];
+            pEval = p[i];
+            TEval = T[i];
 
             int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid];
             int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour];
@@ -566,29 +566,29 @@ namespace Opm
             }
 
             if (pu.phase_used[BlackoilPhases::Liquid]) {
-                LadEval RsSatLad = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad, pLad);
+                Eval RsSatEval = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TEval, pEval);
 
                 double currentRs = 0.0;
                 if (pu.phase_used[BlackoilPhases::Vapour]) {
                     currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset];
                 }
 
-                RsSatLad = std::min(RsSatLad, currentRs);
+                RsSatEval = std::min(RsSatEval, currentRs);
 
-                R[oilOffset] = RsSatLad;
+                R[oilOffset] = RsSatEval;
             }
 
             if (pu.phase_used[BlackoilPhases::Vapour]) {
-                LadEval RvSatLad = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad, pLad);
+                Eval RvSatEval = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TEval, pEval);
 
                 double currentRv = 0.0;
                 if (pu.phase_used[BlackoilPhases::Liquid]) {
                     currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset];
                 }
 
-                RvSatLad = std::min(RvSatLad, currentRv);
+                RvSatEval = std::min(RvSatEval, currentRv);
 
-                R[gasOffset] = RvSatLad;
+                R[gasOffset] = RvSatEval;
             }
         }
     }
@@ -603,19 +603,19 @@ namespace Opm
     {
         const auto& pu = phaseUsage();
 
-        typedef Opm::LocalAd::Evaluation<double, /*size=*/1> LadEval;
-        typedef Opm::MathToolbox<LadEval> Toolbox;
+        typedef Opm::DenseAd::Evaluation<double, /*size=*/1> Eval;
+        typedef Opm::MathToolbox<Eval> Toolbox;
 
-        LadEval pLad = 0.0;
-        LadEval TLad = 0.0;
+        Eval pEval = 0.0;
+        Eval TEval = 0.0;
 
-        pLad.derivatives[0] = 1.0;
+        pEval.derivatives[0] = 1.0;
 
         for (int i = 0; i < n; ++ i) {
             int cellIdx = cells[i];
             int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
-            pLad.value = p[i];
-            TLad.value = T[i];
+            pEval.value = p[i];
+            TEval.value = T[i];
 
             int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid];
             int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour];
@@ -626,31 +626,31 @@ namespace Opm
             }
 
             if (pu.phase_used[BlackoilPhases::Liquid]) {
-                LadEval RsSatLad = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad, pLad);
+                Eval RsSatEval = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TEval, pEval);
 
-                LadEval currentRs = 0.0;
+                Eval currentRs = 0.0;
                 if (pu.phase_used[BlackoilPhases::Vapour]) {
                     currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset];
                 }
 
-                RsSatLad = Toolbox::min(RsSatLad, currentRs);
+                RsSatEval = Toolbox::min(RsSatEval, currentRs);
 
-                R[oilOffset] = RsSatLad.value;
-                dRdp[oilOffset] = RsSatLad.derivatives[0];
+                R[oilOffset] = RsSatEval.value;
+                dRdp[oilOffset] = RsSatEval.derivatives[0];
             }
 
             if (pu.phase_used[BlackoilPhases::Vapour]) {
-                LadEval RvSatLad = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad, pLad);
+                Eval RvSatEval = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TEval, pEval);
 
-                LadEval currentRv = 0.0;
+                Eval currentRv = 0.0;
                 if (pu.phase_used[BlackoilPhases::Liquid]) {
                     currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset];
                 }
 
-                RvSatLad = Toolbox::min(RvSatLad, currentRv);
+                RvSatEval = Toolbox::min(RvSatEval, currentRv);
 
-                R[gasOffset] = RvSatLad.value;
-                dRdp[gasOffset] = RvSatLad.derivatives[0];
+                R[gasOffset] = RvSatEval.value;
+                dRdp[gasOffset] = RvSatEval.derivatives[0];
             }
         }
     }
diff --git a/opm/core/props/phaseUsageFromDeck.hpp b/opm/core/props/phaseUsageFromDeck.hpp
index 8feefd1f0..88571e452 100644
--- a/opm/core/props/phaseUsageFromDeck.hpp
+++ b/opm/core/props/phaseUsageFromDeck.hpp
@@ -88,8 +88,14 @@ namespace Opm
         }
         pu.num_phases = 0;
         for (int i = 0; i < BlackoilPhases::MaxNumPhases; ++i) {
-            pu.phase_pos[i] = pu.num_phases;
-            pu.num_phases += pu.phase_used[i];
+            if (pu.phase_used[i]) {
+                pu.phase_pos[i] = pu.num_phases;
+                pu.num_phases += 1;
+            }
+            else {
+                //Set to ridiculous value on purpose: should never be used
+                pu.phase_pos[i] = 2000000000;
+            }
         }
 
         // Only 2 or 3 phase systems handled.
diff --git a/opm/core/props/rock/RockCompressibility.cpp b/opm/core/props/rock/RockCompressibility.cpp
index 03f384f44..b4f48ad5c 100644
--- a/opm/core/props/rock/RockCompressibility.cpp
+++ b/opm/core/props/rock/RockCompressibility.cpp
@@ -22,6 +22,7 @@
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/utility/Units.hpp>
 #include <opm/common/ErrorMacros.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
 #include <opm/core/utility/linearInterpolation.hpp>
 
 #include <opm/parser/eclipse/EclipseState/Tables/RocktabTable.hpp>
@@ -64,14 +65,18 @@ namespace Opm
             if (rockKeyword.size() != 1) {
                 // here it would be better not to use std::cout directly but to add the
                 // warning to some "warning list"...
-                std::cout << "Can only handle a single region in ROCK ("<<rockKeyword.size()<<" regions specified)."
-                          << " Ignoring all except for the first.\n";
+                OpmLog::warning("Can only handle a single region in ROCK ("
+                                + std::to_string(rockKeyword.size())
+                                + " regions specified)."
+                                + " Ignoring all except for the first.\n" 
+                                + "In file " + rockKeyword.getFileName()
+                                + ", line " + std::to_string(rockKeyword.getLineNumber()) + "\n");
             }
 
             pref_ = rockKeyword.getRecord(0).getItem("PREF").getSIDouble(0);
             rock_comp_ = rockKeyword.getRecord(0).getItem("COMPRESSIBILITY").getSIDouble(0);
         } else {
-            std::cout << "**** warning: no rock compressibility data found in deck (ROCK or ROCKTAB)." << std::endl;
+            OpmLog::warning("No rock compressibility data found in deck (ROCK or ROCKTAB).");
         }
     }
 
@@ -96,7 +101,8 @@ namespace Opm
         if (p_.empty()) {
             // Approximating poro multiplier with a quadratic curve,
             // we must use its derivative.
-            return rock_comp_ + 2 * rock_comp_ * rock_comp_ * (pressure - pref_);
+            const double cpnorm = rock_comp_*(pressure - pref_);
+            return rock_comp_ + cpnorm*rock_comp_;
         } else {
             return Opm::linearInterpolationDerivative(p_, poromult_, pressure);
         }
diff --git a/opm/core/props/satfunc/RelpermDiagnostics.cpp b/opm/core/props/satfunc/RelpermDiagnostics.cpp
index 5cf71b59a..1a6a53072 100644
--- a/opm/core/props/satfunc/RelpermDiagnostics.cpp
+++ b/opm/core/props/satfunc/RelpermDiagnostics.cpp
@@ -24,29 +24,6 @@
 
 namespace Opm{
 
-
-    RelpermDiagnostics::Counter::Counter()
-        :error(0)
-        ,warning(0)
-        ,problem(0)
-        ,bug(0)
-    {
-    }
-
-
-
-
-
-
-    std::vector<std::string>
-    RelpermDiagnostics::getMessages() const
-    {
-        return messages_;
-    }
-
-
-
-
     void RelpermDiagnostics::phaseCheck_(DeckConstPtr deck)
     {
         bool hasWater = deck->hasKeyword("WATER");
@@ -56,31 +33,26 @@ namespace Opm{
             
         if (hasWater && hasGas && !hasOil && !hasSolvent) {
             const std::string msg = "System:  Water-Gas system.";
-            std::cout << msg << std::endl;
             OpmLog::info(msg);
             fluidSystem_ = FluidSystem::WaterGas;
         }
         if (hasWater && hasOil && !hasGas && !hasSolvent) { 
             const std::string msg = "System:  Oil-Water system.";
-            std::cout << msg << std::endl;
             OpmLog::info(msg);
             fluidSystem_ = FluidSystem::OilWater; 
         }
         if (hasOil && hasGas && !hasWater && !hasSolvent) { 
             const std::string msg = "System:  Oil-Gas system.";
-            std::cout << msg << std::endl;
             OpmLog::info(msg);
             fluidSystem_ = FluidSystem::OilGas; 
         }
         if (hasOil && hasWater && hasGas && !hasSolvent) {
             const std::string msg = "System:  Black-oil system.";
-            std::cout << msg << std::endl;
             OpmLog::info(msg);
             fluidSystem_ = FluidSystem::BlackOil;
         }
         if (hasSolvent) {
             const std::string msg = "System:  Solvent model.";
-            std::cout << msg << std::endl;
             OpmLog::info(msg);
             fluidSystem_ = FluidSystem::Solvent;
         }
@@ -107,31 +79,25 @@ namespace Opm{
         bool family2 = ((!swfnTables.empty() && !sgfnTables.empty()) || !sgwfnTables.empty()) && (!sof3Tables.empty() || !sof2Tables.empty());
 
         if (family1 && family2) {
-            const std::string msg = "-- Error:   Saturation families should not be mixed.\n Use either SGOF and SWOF or SGFN, SWFN and SOF3.";
-            messages_.push_back(msg);
+            const std::string msg = "Saturation families should not be mixed.\n Use either SGOF and SWOF or SGFN, SWFN and SOF3.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         if (!family1 && !family2) {
-            const std::string msg = "-- Error,   Saturations function must be specified using either \n \
+            const std::string msg = "Saturations function must be specified using either \n \
                              family 1 or family 2 keywords \n \
                              Use either SGOF and SWOF or SGFN, SWFN and SOF3.";
-            messages_.push_back(msg);
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         if (family1 && !family2) {
             satFamily_ = SaturationFunctionFamily::FamilyI;
             const std::string msg = "Relative permeability input format: Saturation Family I.";
-            std::cout << msg << std::endl;
             OpmLog::info(msg);
         } 
         if (!family1 && family2) {
             satFamily_ = SaturationFunctionFamily::FamilyII;
             const std::string msg = "Relative permeability input format: Saturation Family II.";
-            std::cout << msg << std::endl;
             OpmLog::info(msg);
         }
     }
@@ -144,7 +110,6 @@ namespace Opm{
     {
         const int numSatRegions = deck->getKeyword("TABDIMS").getRecord(0).getItem("NTSFUN").get< int >(0);
         const std::string msg = "Number of saturation regions: " + std::to_string(numSatRegions) + "\n";
-        std::cout << msg << std::endl;
         OpmLog::info(msg);
         const auto& tableManager = eclState->getTableManager();
         const TableContainer& swofTables    = tableManager.getSwofTables();
@@ -217,32 +182,24 @@ namespace Opm{
         const std::string regionIdx = std::to_string(satnumIdx);
         //Check sw column.
         if (sw.front() < 0.0 || sw.back() > 1.0) {
-            const std::string msg = "-- Error:   In SWOF table SATNUM = "+ regionIdx + ", saturation should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SWOF table SATNUM = "+ regionIdx + ", saturation should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         //TODO check endpoint sw.back() == 1. - Sor.
         //Check krw column.
         if (krw.front() != 0.0) {
-            const std::string msg = "-- Error:   In SWOF table SATNUM = " + regionIdx + ", first value of krw should be 0.";
-            messages_.push_back(msg);
+            const std::string msg = "In SWOF table SATNUM = " + regionIdx + ", first value of krw should be 0.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         if (krw.front() < 0.0 || krw.back() > 1.0) {
-            const std::string msg = "-- Error:   In SWOF table SATNUM = " + regionIdx + ", krw should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SWOF table SATNUM = " + regionIdx + ", krw should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         ///Check krow column.
         if (krow.front() > 1.0 || krow.back() < 0.0) {
-            const std::string msg = "-- Error:   In SWOF table SATNUM = "+ regionIdx + ", krow should be in range [0, 1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SWOF table SATNUM = "+ regionIdx + ", krow should be in range [0, 1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         ///TODO check if run with gas.
     }
@@ -260,38 +217,28 @@ namespace Opm{
         const std::string regionIdx = std::to_string(satnumIdx);
         //Check sw column.
         if (sg.front() < 0.0 || sg.back() > 1.0) {
-            const std::string msg = "-- Error:   In SGOF table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SGOF table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         if (sg.front() != 0.0) {
-            const std::string msg = "-- Error:   In SGOF table SATNUM = " + regionIdx + ", first value of sg should be 0.";
-            messages_.push_back(msg);
+            const std::string msg = "In SGOF table SATNUM = " + regionIdx + ", first value of sg should be 0.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         //TODO check endpoint sw.back() == 1. - Sor.
         //Check krw column.
         if (krg.front() != 0.0) {
-            const std::string msg = "-- Error:   In SGOF table SATNUM = " + regionIdx + ", first value of krg should be 0.";
-            messages_.push_back(msg);
+            const std::string msg = "In SGOF table SATNUM = " + regionIdx + ", first value of krg should be 0.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         if (krg.front() < 0.0 || krg.back() > 1.0) {
-            const std::string msg = "-- Error:   In SGOF table SATNUM = " + regionIdx + ", krg should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SGOF table SATNUM = " + regionIdx + ", krg should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check krow column.
         if (krog.front() > 1.0 || krog.back() < 0.0) {
-            const std::string msg = "-- Error:   In SGOF table SATNUM = " + regionIdx + ", krog should be in range [0, 1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SGOF table SATNUM = " + regionIdx + ", krog should be in range [0, 1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         //TODO check if run with water.
     }
@@ -306,36 +253,26 @@ namespace Opm{
         //Check sl column.
         //TODO first value means sl = swco + sor
         if (sl.front() < 0.0 || sl.back() > 1.0) {
-            const std::string msg = "-- Error:   In SLGOF table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SLGOF table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         if (sl.back() != 1.0) {
-            const std::string msg = "-- Error:   In SLGOF table SATNUM = " + regionIdx + ", last value of sl should be 1.";
-            messages_.push_back(msg);
+            const std::string msg = "In SLGOF table SATNUM = " + regionIdx + ", last value of sl should be 1.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         if (krg.front() > 1.0 || krg.back() < 0) {
-            const std::string msg = "-- Error:   In SLGOF table SATNUM = " + regionIdx + ", krg shoule be in range [0, 1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SLGOF table SATNUM = " + regionIdx + ", krg shoule be in range [0, 1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         if (krg.back() != 0.0) {
-            const std::string msg = "-- Error:   In SLGOF table SATNUM = " + regionIdx + ", last value of krg hould be 0.";
-            messages_.push_back(msg);
+            const std::string msg = "In SLGOF table SATNUM = " + regionIdx + ", last value of krg hould be 0.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         if (krog.front() < 0.0 || krog.back() > 1.0) {
-            const std::string msg = "-- Error:   In SLGOF table SATNUM = " + regionIdx + ", krog shoule be in range [0, 1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SLGOF table SATNUM = " + regionIdx + ", krog shoule be in range [0, 1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
     }
 
@@ -351,25 +288,19 @@ namespace Opm{
         const std::string regionIdx = std::to_string(satnumIdx);
         //Check sw column.
         if (sw.front() < 0.0 || sw.back() > 1.0) {
-            const std::string msg = "-- Error:   In SWFN table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SWFN table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         
         //Check krw column.
         if (krw.front() < 0.0 || krw.back() > 1.0) {
-            const std::string msg = "-- Error:   In SWFN table SATNUM = " + regionIdx + ", krw should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SWFN table SATNUM = " + regionIdx + ", krw should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         if (krw.front() != 0.0) {
-            const std::string msg = "-- Error:   In SWFN table SATNUM = " + regionIdx + ", first value of krw should be 0.";
-            messages_.push_back(msg);
+            const std::string msg = "In SWFN table SATNUM = " + regionIdx + ", first value of krw should be 0.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
     }
 
@@ -385,24 +316,18 @@ namespace Opm{
         const std::string regionIdx = std::to_string(satnumIdx);
         //Check sg column.
         if (sg.front() < 0.0 || sg.back() > 1.0) {
-            const std::string msg = "-- Error:   In SGFN table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SGFN table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         
         //Check krg column.
         if (krg.front() < 0.0 || krg.back() > 1.0) {
-            const std::string msg = "-- Error:   In SGFN table SATNUM = " + regionIdx + ", krg should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SGFN table SATNUM = " + regionIdx + ", krg should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         if (krg.front() != 0.0) {
-            const std::string msg = "-- Error:   In SGFN table SATNUM = " + regionIdx + ", first value of krg should be 0.";
-            messages_.push_back(msg);
+            const std::string msg = "In SGFN table SATNUM = " + regionIdx + ", first value of krg should be 0.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
     }
 
@@ -420,46 +345,34 @@ namespace Opm{
         //Check so column.
         //TODO: The max so = 1 - Swco
         if (so.front() < 0.0 || so.back() > 1.0) {
-            const std::string msg = "-- Error:   In SOF3 table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SOF3 table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check krow column.
         if (krow.front() < 0.0 || krow.back() > 1.0) {
-            const std::string msg = "-- Error:   In SOF3 table SATNUM = " + regionIdx + ", krow should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SOF3 table SATNUM = " + regionIdx + ", krow should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         if (krow.front() != 0.0) {
-            const std::string msg = "-- Error:   In SOF3 table SATNUM = " + regionIdx + ", first value of krow should be 0.";
-            messages_.push_back(msg);
+            const std::string msg = "In SOF3 table SATNUM = " + regionIdx + ", first value of krow should be 0.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check krog column.
         if (krog.front() < 0.0 || krog.back() > 1.0) {
-            const std::string msg = "-- Error:   In SOF3 table SATNUM = " + regionIdx + ", krog should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SOF3 table SATNUM = " + regionIdx + ", krog should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         if (krog.front() != 0.0) {
-            const std::string msg = "-- Error:   In SOF3 table SATNUM = " + regionIdx + ", first value of krog should be 0.";
-            messages_.push_back(msg);
+            const std::string msg = "In SOF3 table SATNUM = " + regionIdx + ", first value of krog should be 0.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
     
         if (krog.back() != krow.back()) {
-            const std::string msg = "-- Error:   In SOF3 table SATNUM = " + regionIdx + ", max value of krog and krow should be the same.";
-            messages_.push_back(msg);
+            const std::string msg = "In SOF3 table SATNUM = " + regionIdx + ", max value of krog and krow should be the same.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
     }
 
@@ -476,24 +389,18 @@ namespace Opm{
         //Check so column.
         //TODO: The max so = 1 - Swco
         if (so.front() < 0.0 || so.back() > 1.0) {
-            const std::string msg = "-- Error:   In SOF2 table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SOF2 table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check krow column.
         if (kro.front() < 0.0 || kro.back() > 1.0) {
-            const std::string msg = "-- Error:   In SOF2 table SATNUM = " + regionIdx + ", krow should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SOF2 table SATNUM = " + regionIdx + ", krow should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         if (kro.front() != 0.0) {
-            const std::string msg = "-- Error:   In SOF2 table SATNUM = " + regionIdx + ", first value of krow should be 0.";
-            messages_.push_back(msg);
+            const std::string msg = "In SOF2 table SATNUM = " + regionIdx + ", first value of krow should be 0.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
     }
 
@@ -510,39 +417,29 @@ namespace Opm{
         const std::string regionIdx = std::to_string(satnumIdx);
         //Check sg column.
         if (sg.front() < 0.0 || sg.back() > 1.0) {
-            const std::string msg = "-- Error:   In SGWFN table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SGWFN table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check krg column.
         if (krg.front() < 0.0 || krg.back() > 1.0) {
-            const std::string msg = "-- Error:   In SGWFN table SATNUM = " + regionIdx + ", krg should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SGWFN table SATNUM = " + regionIdx + ", krg should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         if (krg.front() != 0.0) {
-            const std::string msg = "-- Error:   In SGWFN table SATNUM = " + regionIdx + ", first value of krg should be 0.";
-            messages_.push_back(msg);
+            const std::string msg = "In SGWFN table SATNUM = " + regionIdx + ", first value of krg should be 0.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check krgw column.
         //TODO check saturation sw = 1. - sg
         if (krgw.front() > 1.0 || krgw.back() < 0.0) {
-            const std::string msg = "-- Error:   In SGWFN table SATNUM = " + regionIdx + ", krgw should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SGWFN table SATNUM = " + regionIdx + ", krgw should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
         if (krgw.back() != 0.0) {
-            const std::string msg = "-- Error:   In SGWFN table SATNUM = " + regionIdx + ", last value of krgw should be 0.";
-            messages_.push_back(msg);
+            const std::string msg = "In SGWFN table SATNUM = " + regionIdx + ", last value of krgw should be 0.";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
     }
 
@@ -556,18 +453,14 @@ namespace Opm{
         const std::string regionIdx = std::to_string(satnumIdx);
         //Check sw column.
         if (sw.front() < 0.0 || sw.back() > 1.0) {
-            const std::string msg = "-- Error:   In SGCWMIS table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SGCWMIS table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check critical gas column.
         if (sgc.front() < 0.0 || sgc.back() > 1.0) {
-            const std::string msg = "-- Error:   In SGCWMIS table SATNUM = " + regionIdx + ", critical gas saturation should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SGCWMIS table SATNUM = " + regionIdx + ", critical gas saturation should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
     }
 
@@ -583,18 +476,14 @@ namespace Opm{
         const std::string regionIdx = std::to_string(satnumIdx);
         //Check sw column.
         if (sw.front() < 0.0 || sw.back() > 1.0) {
-            const std::string msg = "-- Error:   In SORWMIS table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SORWMIS table SATNUM = " + regionIdx + ", saturation should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check critical oil column.
         if (sor.front() < 0.0 || sor.back() > 1.0) {
-            const std::string msg = "-- Error:   In SORWMIS table SATNUM = " + regionIdx + ", critical oil saturation should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SORWMIS table SATNUM = " + regionIdx + ", critical oil saturation should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
     }
 
@@ -610,26 +499,20 @@ namespace Opm{
         const std::string regionIdx = std::to_string(satnumIdx);
         //Check phase fraction column.
         if (frac.front() < 0.0 || frac.back() > 1.0) {
-            const std::string msg = "-- Error:   In SSFN table SATNUM = " + regionIdx + ", phase fraction should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SSFN table SATNUM = " + regionIdx + ", phase fraction should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check gas relperm multiplier column.
         if (krgm.front() < 0.0 || krgm.back() > 1.0) {
-            const std::string msg = "-- Error:   In SSFN table SATNUM = " + regionIdx + ", gas relative permeability multiplier should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SSFN table SATNUM = " + regionIdx + ", gas relative permeability multiplier should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check solvent relperm multiplier column.
         if (krsm.front() < 0.0 || krsm.back() > 1.0) {
-            const std::string msg = "-- Error:   In SSFN table SATNUM = " + regionIdx + ", solvent relative permeability multiplier should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In SSFN table SATNUM = " + regionIdx + ", solvent relative permeability multiplier should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
     }
 
@@ -647,18 +530,14 @@ namespace Opm{
         const std::string regionIdx = std::to_string(satnumIdx);
         //Check phase fraction column.
         if (frac.front() < 0.0 || frac.back() > 1.0) {
-            const std::string msg = "-- Error:   In MISC table SATNUM = " + regionIdx + ", phase fraction should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In MISC table SATNUM = " + regionIdx + ", phase fraction should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check miscibility column.
         if (misc.front() < 0.0 || misc.back() > 1.0) {
-            const std::string msg = "-- Error:   In MISC table SATNUM = " + regionIdx + ", miscibility should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In MISC table SATNUM = " + regionIdx + ", miscibility should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
     }
 
@@ -676,26 +555,20 @@ namespace Opm{
         const std::string regionIdx = std::to_string(satnumIdx);
         //Check phase fraction column.
         if (frac.front() < 0.0 || frac.back() > 1.0) {
-            const std::string msg = "-- Error:   In MSFN table SATNUM = " + regionIdx + ", total gas fraction should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In MSFN table SATNUM = " + regionIdx + ", total gas fraction should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check gas_solvent relperm multiplier column.
         if (krgsm.front() < 0.0 || krgsm.back() > 1.0) {
-            const std::string msg = "-- Error:   In MSFN table SATNUM = " + regionIdx + ", gas+solvent relative permeability multiplier should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In MSFN table SATNUM = " + regionIdx + ", gas+solvent relative permeability multiplier should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
 
         //Check oil relperm multiplier column.
         if (krom.front() > 1.0 || krom.back() < 0.0) {
-            const std::string msg = "-- Error:   In MSFN table SATNUM = " + regionIdx + ", oil relative permeability multiplier should be in range [0,1].";
-            messages_.push_back(msg);
+            const std::string msg = "In MSFN table SATNUM = " + regionIdx + ", oil relative permeability multiplier should be in range [0,1].";
             OpmLog::error(msg);
-            counter_.error += 1;
         }
     }
 
@@ -722,16 +595,12 @@ namespace Opm{
              const std::string regionIdx = std::to_string(satnumIdx + 1);
              ///Consistency check.
              if (unscaledEpsInfo_[satnumIdx].Sgu > (1. - unscaledEpsInfo_[satnumIdx].Swl)) {
-                const std::string msg = "-- Warning: In saturation table SATNUM = " + regionIdx + ", Sgmax should not exceed 1-Swco.";
-                messages_.push_back(msg);
+                const std::string msg = "In saturation table SATNUM = " + regionIdx + ", Sgmax should not exceed 1-Swco.";
                 OpmLog::warning(msg);
-                counter_.warning += 1;
              }
              if (unscaledEpsInfo_[satnumIdx].Sgl > (1. - unscaledEpsInfo_[satnumIdx].Swu)) {
-                const std::string msg = "-- Warning: In saturation table SATNUM = " + regionIdx + ", Sgco should not exceed 1-Swmax.";
-                messages_.push_back(msg);
+                const std::string msg = "In saturation table SATNUM = " + regionIdx + ", Sgco should not exceed 1-Swmax.";
                 OpmLog::warning(msg);
-                counter_.warning += 1;
              }
 
              //Krow(Sou) == Krog(Sou) for three-phase
@@ -762,25 +631,19 @@ namespace Opm{
                      krog_value = table.evaluate("KROG" , Sou);
                  }
                  if (krow_value != krog_value) {
-                     const std::string msg = "-- Warning: In saturation table SATNUM = " + regionIdx + ", Krow(Somax) should be equal to Krog(Somax).";
-                     messages_.push_back(msg);
+                     const std::string msg = "In saturation table SATNUM = " + regionIdx + ", Krow(Somax) should be equal to Krog(Somax).";
                      OpmLog::warning(msg);
-                     counter_.warning += 1;
                  }
              }
              //Krw(Sw=0)=Krg(Sg=0)=Krow(So=0)=Krog(So=0)=0.
              //Mobile fluid requirements
             if (((unscaledEpsInfo_[satnumIdx].Sowcr + unscaledEpsInfo_[satnumIdx].Swcr)-1) >= 0) {
-                const std::string msg = "-- Warning: In saturation table SATNUM = " + regionIdx + ", Sowcr + Swcr should be less than 1.";
-                messages_.push_back(msg);
+                const std::string msg = "In saturation table SATNUM = " + regionIdx + ", Sowcr + Swcr should be less than 1.";
                 OpmLog::warning(msg);
-                counter_.warning += 1;
             }
             if (((unscaledEpsInfo_[satnumIdx].Sogcr + unscaledEpsInfo_[satnumIdx].Sgcr + unscaledEpsInfo_[satnumIdx].Swl) - 1 ) > 0) {
-                const std::string msg = "-- Warning: In saturation table SATNUM = " + regionIdx + ", Sogcr + Sgcr + Swco should be less than 1.";
-                messages_.push_back(msg);
+                const std::string msg = "In saturation table SATNUM = " + regionIdx + ", Sogcr + Sgcr + Swco should be less than 1.";
                 OpmLog::warning(msg);
-                counter_.warning += 1;
             }
         }
     }
diff --git a/opm/core/props/satfunc/RelpermDiagnostics.hpp b/opm/core/props/satfunc/RelpermDiagnostics.hpp
index 191998582..cb5cf138f 100644
--- a/opm/core/props/satfunc/RelpermDiagnostics.hpp
+++ b/opm/core/props/satfunc/RelpermDiagnostics.hpp
@@ -62,8 +62,6 @@ namespace Opm {
                        DeckConstPtr deck,
                        const GridT& grid);
 
-        std::vector<std::string> getMessages() const;
-
     private:
         enum FluidSystem {
             OilWater,
@@ -83,22 +81,9 @@ namespace Opm {
   
         SaturationFunctionFamily satFamily_;
 
-        struct Counter {
-            Counter();
-            int error;
-            int warning;
-            int problem;
-            int bug;
-        };
-
-        Counter counter_;
-        
         std::vector<Opm::EclEpsScalingPointsInfo<double> > unscaledEpsInfo_;
         std::vector<Opm::EclEpsScalingPointsInfo<double> > scaledEpsInfo_;
 
-        std::vector<std::string> messages_;
-        ///Store scaled information.
-        std::vector<std::string> scaled_messages_;
 
         ///Check the phase that used.
         void phaseCheck_(DeckConstPtr deck);
diff --git a/opm/core/props/satfunc/RelpermDiagnostics_impl.hpp b/opm/core/props/satfunc/RelpermDiagnostics_impl.hpp
index 869e8cbb5..2f72a2fcd 100644
--- a/opm/core/props/satfunc/RelpermDiagnostics_impl.hpp
+++ b/opm/core/props/satfunc/RelpermDiagnostics_impl.hpp
@@ -33,41 +33,12 @@ namespace Opm {
                                        Opm::DeckConstPtr deck,
                                        const GridT& grid)
     {
-        std::cout << "\n\n***************Saturation Functions Diagnostics***************\n\n";
+        OpmLog::info("\n***************Saturation Functions Diagnostics***************");
         phaseCheck_(deck);
         satFamilyCheck_(eclState);
         tableCheck_(eclState, deck);
         unscaledEndPointsCheck_(deck, eclState);
         scaledEndPointsCheck_(deck, eclState, grid);
-        if (!messages_.empty()) {
-            std::sort(messages_.begin(), messages_.end());
-            auto it = std::unique(messages_.begin(), messages_.end());
-            messages_.erase(it, messages_.end());
-            std::cout << std::endl;
-            for (const auto& x : messages_) {
-                std::cout << x << std::endl;
-            }
-        }
-        int limits = 0;
-        if (!scaled_messages_.empty()) {
-            std::cout << std::endl;
-            for (const auto& x : scaled_messages_) {
-                if (limits < 10) {
-                    std::cout << x << std::endl;
-                    limits++;
-                } else {
-                    std::cout << "\nMore inconsistencies exist. Check saturation function input and LOGFILE!" << std::endl;
-                    break;
-                }
-            }
-        }
-
-        const std::string summary_msg = "\n\nError summary:" + 
-            std::string("\nWarnings          " + std::to_string(counter_.warning)) +
-            std::string("\nProblems          " + std::to_string(counter_.problem)) +
-            std::string("\nErrors            " + std::to_string(counter_.error)) + 
-            std::string("\nBugs              " + std::to_string(counter_.bug))+ "\n";
-        std::cout << summary_msg << std::endl;
     }
 
     template <class GridT>
@@ -83,7 +54,8 @@ namespace Opm {
         EclEpsGridProperties epsGridProperties;
         epsGridProperties.initFromDeck(deck, eclState, /*imbibition=*/false);       
         const auto& satnum = eclState->get3DProperties().getIntGridProperty("SATNUM");
-
+        
+        const std::string tag = "Scaled endpoints";
         for (int c = 0; c < nc; ++c) {
             const int cartIdx = compressedToCartesianIdx[c];
             const std::string satnumIdx = std::to_string(satnum.iget(cartIdx));
@@ -98,82 +70,62 @@ namespace Opm {
 
             // SGU <= 1.0 - SWL
             if (scaledEpsInfo_[c].Sgu > (1.0 - scaledEpsInfo_[c].Swl)) {
-                const std::string msg = "-- Warning: For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SGU exceed 1.0 - SWL";
-                scaled_messages_.push_back(msg);
-                OpmLog::warning(msg);
-                counter_.warning += 1;
+                const std::string msg = "For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SGU exceed 1.0 - SWL";
+                OpmLog::warning(tag, msg);
             }
             
             // SGL <= 1.0 - SWU
             if (scaledEpsInfo_[c].Sgl > (1.0 - scaledEpsInfo_[c].Swu)) {
-                const std::string msg = "-- Warning: For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SGL exceed 1.0 - SWU";
-                scaled_messages_.push_back(msg);
-                OpmLog::warning(msg);
-                counter_.warning += 1;
+                const std::string msg = "For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SGL exceed 1.0 - SWU";
+                OpmLog::warning(tag, msg);
             }
 
             if (deck->hasKeyword("SCALECRS") && fluidSystem_ == FluidSystem::BlackOil) {
                 // Mobilility check.
                 if ((scaledEpsInfo_[c].Sowcr + scaledEpsInfo_[c].Swcr) >= 1.0) {
-                    const std::string msg = "-- Warning: For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SOWCR + SWCR exceed 1.0";
-                    scaled_messages_.push_back(msg);
-                    OpmLog::warning(msg);
-                    counter_.warning += 1;
+                    const std::string msg = "For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SOWCR + SWCR exceed 1.0";
+                    OpmLog::warning(tag, msg);
                 }
 
                 if ((scaledEpsInfo_[c].Sogcr + scaledEpsInfo_[c].Sgcr + scaledEpsInfo_[c].Swl) >= 1.0) {
-                    const std::string msg = "-- Warning: For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SOGCR + SGCR + SWL exceed 1.0";
-                    scaled_messages_.push_back(msg);
-                    OpmLog::warning(msg);
-                    counter_.warning += 1;
+                    const std::string msg = "For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SOGCR + SGCR + SWL exceed 1.0";
+                    OpmLog::warning(tag, msg);
                 }
             }
             ///Following rules come from NEXUS.
             if (fluidSystem_ != FluidSystem::WaterGas) {
                 if (scaledEpsInfo_[c].Swl > scaledEpsInfo_[c].Swcr) {
-                    const std::string msg = "-- Warning: For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SWL > SWCR";
-                    scaled_messages_.push_back(msg);
-                    OpmLog::warning(msg);
-                    counter_.warning += 1;
+                    const std::string msg = "For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SWL > SWCR";
+                    OpmLog::warning(tag, msg);
                 }
 
                 if (scaledEpsInfo_[c].Swcr > scaledEpsInfo_[c].Sowcr) {
-                    const std::string msg = "-- Warning: For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SWCR > SOWCR";
-                    scaled_messages_.push_back(msg);
-                    OpmLog::warning(msg);
-                    counter_.warning += 1;
+                    const std::string msg = "For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SWCR > SOWCR";
+                    OpmLog::warning(tag, msg);
                 }
             
                 if (scaledEpsInfo_[c].Sowcr > scaledEpsInfo_[c].Swu) {
-                    const std::string msg = "-- Warning: For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SOWCR > SWU";
-                    scaled_messages_.push_back(msg);
-                    OpmLog::warning(msg);
-                    counter_.warning += 1;
+                    const std::string msg = "For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SOWCR > SWU";
+                    OpmLog::warning(tag, msg);
                 }
             }
 
             if (fluidSystem_ != FluidSystem::OilWater) {
                 if (scaledEpsInfo_[c].Sgl > scaledEpsInfo_[c].Sgcr) {
-                    const std::string msg = "-- Warning: For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SGL > SGCR";
-                    scaled_messages_.push_back(msg);
-                    OpmLog::warning(msg);
-                    counter_.warning += 1;
+                    const std::string msg = "For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SGL > SGCR";
+                    OpmLog::warning(tag, msg);
                 }
             }
 
             if (fluidSystem_ != FluidSystem::BlackOil) {
                 if (scaledEpsInfo_[c].Sgcr > scaledEpsInfo_[c].Sogcr) {
-                    const std::string msg = "-- Warning: For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SGCR > SOGCR";
-                    scaled_messages_.push_back(msg);
-                    OpmLog::warning(msg);
-                    counter_.warning += 1;
+                    const std::string msg = "For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SGCR > SOGCR";
+                    OpmLog::warning(tag, msg);
                 }
 
                 if (scaledEpsInfo_[c].Sogcr > scaledEpsInfo_[c].Sgu) {
-                    const std::string msg = "-- Warning: For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SOGCR > SGU";
-                    scaled_messages_.push_back(msg);
-                    OpmLog::warning(msg);
-                    counter_.warning += 1;
+                    const std::string msg = "For scaled endpoints input, cell" + cellIdx + " SATNUM = " + satnumIdx + ", SOGCR > SGU";
+                    OpmLog::warning(tag, msg);
                 }
             }
         } 
diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck.cpp b/opm/core/props/satfunc/SaturationPropsFromDeck.cpp
index db7ab6221..48b1feb38 100644
--- a/opm/core/props/satfunc/SaturationPropsFromDeck.cpp
+++ b/opm/core/props/satfunc/SaturationPropsFromDeck.cpp
@@ -136,9 +136,11 @@ namespace Opm
                                            double* pc,
                                            double* dpcds) const
     {
-        assert(cells != 0);
+        assert(cells != 0);        
+        assert(phaseUsage_.phase_used[BlackoilPhases::Liquid]);
 
         const int np = numPhases();
+
         if (dpcds) {
             ExplicitArraysSatDerivativesFluidState fluidState(phaseUsage_);
             typedef ExplicitArraysSatDerivativesFluidState::Evaluation Evaluation;
@@ -151,12 +153,26 @@ namespace Opm
                 MaterialLaw::capillaryPressures(capillaryPressures, params, fluidState);
 
                 // copy the values calculated using opm-material to the target arrays
-                for (int pcPhaseIdx = 0; pcPhaseIdx < np; ++pcPhaseIdx) {
-                    double sign = (pcPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0;
-                    pc[np*i + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx].value;
+                for (int canonicalPhaseIdx = 0; canonicalPhaseIdx < BlackoilPhases::MaxNumPhases; ++canonicalPhaseIdx) {
+                    // skip unused phases
+                    if ( ! phaseUsage_.phase_used[canonicalPhaseIdx]) {
+                        continue;
+                    }
+                    const int pcPhaseIdx = phaseUsage_.phase_pos[canonicalPhaseIdx];
 
-                    for (int satPhaseIdx = 0; satPhaseIdx < np; ++satPhaseIdx)
-                        dpcds[np*np*i + satPhaseIdx*np + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx].derivatives[satPhaseIdx];
+                    const double sign = (canonicalPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0;
+                    // in opm-material the wetting phase is the reference phase
+                    // for two-phase problems i.e water for oil-water system,
+                    // but for flow it is always oil. Add oil (liquid) capillary pressure value
+                    // to shift the reference phase to oil
+                    pc[np*i + pcPhaseIdx] = capillaryPressures[BlackoilPhases::Liquid].value + sign * capillaryPressures[canonicalPhaseIdx].value;
+                    for (int canonicalSatPhaseIdx = 0; canonicalSatPhaseIdx < BlackoilPhases::MaxNumPhases; ++canonicalSatPhaseIdx) {
+                        if ( ! phaseUsage_.phase_used[canonicalSatPhaseIdx])
+                            continue;
+
+                        const int satPhaseIdx = phaseUsage_.phase_pos[canonicalSatPhaseIdx];
+                        dpcds[np*np*i + satPhaseIdx*np + pcPhaseIdx] = capillaryPressures[BlackoilPhases::Liquid].derivatives[canonicalSatPhaseIdx] + sign * capillaryPressures[canonicalPhaseIdx].derivatives[canonicalSatPhaseIdx];
+                    }
                 }
             }
         } else {
@@ -170,9 +186,18 @@ namespace Opm
                 MaterialLaw::capillaryPressures(capillaryPressures, params, fluidState);
 
                 // copy the values calculated using opm-material to the target arrays
-                for (int pcPhaseIdx = 0; pcPhaseIdx < np; ++pcPhaseIdx) {
-                    double sign = (pcPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0;
-                    pc[np*i + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx];
+                for (int canonicalPhaseIdx = 0; canonicalPhaseIdx < BlackoilPhases::MaxNumPhases; ++canonicalPhaseIdx) {
+                    // skip unused phases
+                    if ( ! phaseUsage_.phase_used[canonicalPhaseIdx])
+                        continue;
+
+                    const int pcPhaseIdx = phaseUsage_.phase_pos[canonicalPhaseIdx];
+                    double sign = (canonicalPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0;
+                    // in opm-material the wetting phase is the reference phase
+                    // for two-phase problems i.e water for oil-water system,
+                    // but for flow it is always oil. Add oil (liquid) capillary pressure value
+                    // to shift the reference phase to oil
+                    pc[np*i + pcPhaseIdx] = capillaryPressures[BlackoilPhases::Liquid] + sign * capillaryPressures[canonicalPhaseIdx];
                 }
             }
         }
diff --git a/opm/core/simulator/BlackoilState.cpp b/opm/core/simulator/BlackoilState.cpp
index 1f51dabba..97d39f72b 100644
--- a/opm/core/simulator/BlackoilState.cpp
+++ b/opm/core/simulator/BlackoilState.cpp
@@ -1,6 +1,8 @@
 #include "BlackoilState.hpp"
 #include <opm/common/util/numeric/cmp.hpp>
 #include <opm/core/props/BlackoilPropertiesInterface.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/output/Wells.hpp>
 
 
 using namespace Opm;
@@ -21,15 +23,18 @@ BlackoilState::BlackoilState( size_t num_cells , size_t num_faces , size_t num_p
 }
 
 BlackoilState::BlackoilState( const BlackoilState& other )
-    : SimulationDataContainer(other)
+    : SimulationDataContainer(other),
+      hydrocarbonstate_(other.hydroCarbonState())
 {
     setBlackoilStateReferencePointers();
+
 }
 
 BlackoilState& BlackoilState::operator=( const BlackoilState& other )
 {
     SimulationDataContainer::operator=(other);
     setBlackoilStateReferencePointers();
+    hydrocarbonstate_ = other.hydroCarbonState();
     return *this;
 }
 
diff --git a/opm/core/simulator/BlackoilState.hpp b/opm/core/simulator/BlackoilState.hpp
index 7f45e029f..fe6c0dbc5 100644
--- a/opm/core/simulator/BlackoilState.hpp
+++ b/opm/core/simulator/BlackoilState.hpp
@@ -30,6 +30,12 @@
 namespace Opm
 {
 
+    enum HydroCarbonState {
+        GasOnly = 0,
+        GasAndOil = 1,
+        OilOnly = 2
+    };
+
     /// Simulator state for a blackoil simulator.
     class BlackoilState : public SimulationDataContainer
     {
@@ -62,10 +68,12 @@ namespace Opm
         std::vector<double>& surfacevol  () { return *surfacevol_ref_;  }
         std::vector<double>& gasoilratio () { return *gasoilratio_ref_; }
         std::vector<double>& rv ()          { return *rv_ref_;          }
+        std::vector<HydroCarbonState>& hydroCarbonState() { return hydrocarbonstate_;  }
 
         const std::vector<double>& surfacevol  () const { return *surfacevol_ref_;  }
         const std::vector<double>& gasoilratio () const { return *gasoilratio_ref_; }
         const std::vector<double>& rv ()          const { return *rv_ref_;          }
+        const std::vector<HydroCarbonState>& hydroCarbonState() const { return hydrocarbonstate_;  }
 
     private:
         void setBlackoilStateReferencePointers();
@@ -73,6 +81,9 @@ namespace Opm
         std::vector<double>* gasoilratio_ref_;
         std::vector<double>* rv_ref_;
 
+        // A vector storing the hydro carbon state.
+        std::vector<HydroCarbonState> hydrocarbonstate_;
+
 
     };
 } // namespace Opm
diff --git a/opm/core/simulator/ExplicitArraysSatDerivativesFluidState.hpp b/opm/core/simulator/ExplicitArraysSatDerivativesFluidState.hpp
index d85856ecd..c6b5aa26c 100644
--- a/opm/core/simulator/ExplicitArraysSatDerivativesFluidState.hpp
+++ b/opm/core/simulator/ExplicitArraysSatDerivativesFluidState.hpp
@@ -20,8 +20,8 @@
 #ifndef OPM_EXPLICIT_ARRAYS_SAT_DERIVATIVES_FLUID_STATE_HEADER_INCLUDED
 #define OPM_EXPLICIT_ARRAYS_SAT_DERIVATIVES_FLUID_STATE_HEADER_INCLUDED
 
-#include <opm/material/localad/Evaluation.hpp>
-#include <opm/material/localad/Math.hpp>
+#include <opm/material/densead/Evaluation.hpp>
+#include <opm/material/densead/Math.hpp>
 
 #include <opm/core/props/BlackoilPhases.hpp>
 
@@ -44,7 +44,7 @@ public:
     enum { numPhases = BlackoilPhases::MaxNumPhases };
     enum { numComponents = 3 };
 
-    typedef Opm::LocalAd::Evaluation<double, numPhases> Evaluation;
+    typedef Opm::DenseAd::Evaluation<double, numPhases> Evaluation;
     typedef Evaluation Scalar;
 
     ExplicitArraysSatDerivativesFluidState(const PhaseUsage& phaseUsage)
diff --git a/opm/core/simulator/SimulatorReport.cpp b/opm/core/simulator/SimulatorReport.cpp
index 2ddca2857..e860c946c 100644
--- a/opm/core/simulator/SimulatorReport.cpp
+++ b/opm/core/simulator/SimulatorReport.cpp
@@ -27,6 +27,7 @@ namespace Opm
         : pressure_time(0.0),
           transport_time(0.0),
           total_time(0.0),
+          total_linearizations( 0 ),
           total_newton_iterations( 0 ),
           total_linear_iterations( 0 ),
           verbose_(verbose)
@@ -61,6 +62,7 @@ namespace Opm
         {
             os << "Total time taken (seconds):  " << total_time
                << "\nSolver time (seconds):       " << pressure_time
+               << "\nOverall Linearizations:      " << total_linearizations
                << "\nOverall Newton Iterations:   " << total_newton_iterations
                << "\nOverall Linear Iterations:   " << total_linear_iterations
                << std::endl;
diff --git a/opm/core/simulator/SimulatorReport.hpp b/opm/core/simulator/SimulatorReport.hpp
index 22043918d..9e35905ad 100644
--- a/opm/core/simulator/SimulatorReport.hpp
+++ b/opm/core/simulator/SimulatorReport.hpp
@@ -32,6 +32,7 @@ namespace Opm
         double transport_time;
         double total_time;
 
+        unsigned int total_linearizations;
         unsigned int total_newton_iterations;
         unsigned int total_linear_iterations;
 
diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp
index bbf0694b2..7fc82cc12 100644
--- a/opm/core/simulator/WellState.hpp
+++ b/opm/core/simulator/WellState.hpp
@@ -22,8 +22,11 @@
 
 #include <opm/core/wells.h>
 #include <opm/core/well_controls.h>
+#include <opm/output/Wells.hpp>
+
 #include <array>
 #include <map>
+#include <memory>
 #include <string>
 #include <vector>
 #include <cassert>
@@ -49,6 +52,7 @@ namespace Opm
         {
             // clear old name mapping
             wellMap_.clear();
+            wells_.reset( clone_wells( wells ) );
 
             if (wells) {
                 const int nw = wells->number_of_wells;
@@ -208,6 +212,43 @@ namespace Opm
             return wellRates().size() / numWells();
         }
 
+        virtual data::Wells report() const
+        {
+            return { { /* WellState offers no completion data, so that has to be added later */ },
+                this->bhp(),
+                this->temperature(),
+                this->wellRates(),
+                this->perfPress(),
+                this->perfRates() };
+        }
+
+        virtual ~WellState() {}
+
+        WellState() = default;
+        WellState( const WellState& rhs ) :
+            bhp_( rhs.bhp_ ),
+            thp_( rhs.thp_ ),
+            temperature_( rhs.temperature_ ),
+            wellrates_( rhs.wellrates_ ),
+            perfrates_( rhs.perfrates_ ),
+            perfpress_( rhs.perfpress_ ),
+            wellMap_( rhs.wellMap_ ),
+            wells_( clone_wells( rhs.wells_.get() ) )
+        {}
+
+        WellState& operator=( const WellState& rhs ) {
+            this->bhp_ = rhs.bhp_;
+            this->thp_ = rhs.thp_;
+            this->temperature_ = rhs.temperature_;
+            this->wellrates_ = rhs.wellrates_;
+            this->perfrates_ = rhs.perfrates_;
+            this->perfpress_ = rhs.perfpress_;
+            this->wellMap_ = rhs.wellMap_;
+            this->wells_.reset( clone_wells( rhs.wells_.get() ) );
+
+            return *this;
+        }
+
     private:
         std::vector<double> bhp_;
         std::vector<double> thp_;
@@ -217,6 +258,12 @@ namespace Opm
         std::vector<double> perfpress_;
 
         WellMapType wellMap_;
+
+    protected:
+        struct wdel {
+            void operator()( Wells* w ) { destroy_wells( w ); }
+        };
+        std::unique_ptr< Wells, wdel > wells_;
     };
 
 } // namespace Opm
diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp
index 2e539ea49..cbda500ee 100644
--- a/opm/core/simulator/initStateEquil.hpp
+++ b/opm/core/simulator/initStateEquil.hpp
@@ -37,6 +37,7 @@
 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/RsvdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/RvvdTable.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
 
 #include <array>
 #include <cassert>
@@ -199,7 +200,7 @@ namespace Opm
             std::vector<EquilRecord>
             getEquil(const Opm::EclipseState& state)
             {
-                const auto& init = *state.getInitConfig();
+                const auto& init = state.getInitConfig();
 
                 if( !init.hasEquil() ) {
                     OPM_THROW(std::domain_error, "Deck does not provide equilibration data.");
@@ -264,6 +265,11 @@ namespace Opm
                     if (deck->hasKeyword("DISGAS")) {                    
                         const TableContainer& rsvdTables = tables.getRsvdTables();
                         for (size_t i = 0; i < rec.size(); ++i) {
+                            if (eqlmap.cells(i).empty())
+                            {
+                                rs_func_.push_back(std::shared_ptr<Miscibility::RsVD>());
+                                continue;
+                            }
                             const int cell = *(eqlmap.cells(i).begin());                   
                             if (!rec[i].liveOilInitConstantRs()) {
                                 if (rsvdTables.size() <= 0 ) {
@@ -297,6 +303,11 @@ namespace Opm
                     if (deck->hasKeyword("VAPOIL")) {                    
                         const TableContainer& rvvdTables = tables.getRvvdTables();
                         for (size_t i = 0; i < rec.size(); ++i) {
+                            if (eqlmap.cells(i).empty())
+                            {
+                                rv_func_.push_back(std::shared_ptr<Miscibility::RvVD>());
+                                continue;
+                            }
                             const int cell = *(eqlmap.cells(i).begin());                   
                             if (!rec[i].wetGasInitConstantRv()) {
                                 if (rvvdTables.size() <= 0) { 
@@ -381,6 +392,12 @@ namespace Opm
                 {
                     for (const auto& r : reg.activeRegions()) {
                         const auto& cells = reg.cells(r);
+                        if (cells.empty())
+                        {
+                            OpmLog::warning("Equilibration region " + std::to_string(r + 1) 
+                                            + " has no active cells");
+                            continue;
+                        }
                         const int repcell = *cells.begin();
 
                         const RhoCalc calc(props, repcell);
diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp
index b3d750c34..20f56c5a4 100644
--- a/opm/core/simulator/initStateEquil_impl.hpp
+++ b/opm/core/simulator/initStateEquil_impl.hpp
@@ -767,28 +767,34 @@ namespace Opm
                 double sat[BlackoilPhases::MaxNumPhases];
                 double threshold_sat = 1.0e-6;
 
-                sat[waterpos] = smax[waterpos];
-                sat[gaspos] = smax[gaspos];
-                sat[oilpos] = 1.0 - sat[waterpos] - sat[gaspos];
-                if (sw > smax[waterpos]-threshold_sat ) {
-                    sat[waterpos] = smax[waterpos];
-                    props.capPress(1, sat, &cell, pc, 0);                   
-                    phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
-                } else if (sg > smax[gaspos]-threshold_sat) {
-                    sat[gaspos] = smax[gaspos];
-                    props.capPress(1, sat, &cell, pc, 0);                   
-                    phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
-                }
-                if (sg < smin[gaspos]+threshold_sat) {
-                    sat[gaspos] = smin[gaspos];
-                    props.capPress(1, sat, &cell, pc, 0);
-                    phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
-                }
-                if (sw < smin[waterpos]+threshold_sat) {
-                    sat[waterpos] = smin[waterpos];
-                    props.capPress(1, sat, &cell, pc, 0);
-                    phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
-                }
+                sat[oilpos] = 1.0;
+                  if (water) {
+                      sat[waterpos] = smax[waterpos];
+                      sat[oilpos] -= sat[waterpos];
+                  }
+                  if (gas) {
+                      sat[gaspos] = smax[gaspos];
+                      sat[oilpos] -= sat[gaspos];
+                  }
+                  if (water && sw > smax[waterpos]-threshold_sat ) {
+                      sat[waterpos] = smax[waterpos];
+                      props.capPress(1, sat, &cell, pc, 0);
+                      phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
+                  } else if (gas && sg > smax[gaspos]-threshold_sat) {
+                      sat[gaspos] = smax[gaspos];
+                      props.capPress(1, sat, &cell, pc, 0);
+                      phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
+                  }
+                  if (gas && sg < smin[gaspos]+threshold_sat) {
+                      sat[gaspos] = smin[gaspos];
+                      props.capPress(1, sat, &cell, pc, 0);
+                      phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
+                  }
+                  if (water && sw < smin[waterpos]+threshold_sat) {
+                      sat[waterpos] = smin[waterpos];
+                      props.capPress(1, sat, &cell, pc, 0);
+                      phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
+                  }
             }
             return phase_saturations;
         }
diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp
index 8e8fc24f9..910262a2a 100644
--- a/opm/core/simulator/initState_impl.hpp
+++ b/opm/core/simulator/initState_impl.hpp
@@ -741,7 +741,7 @@ namespace Opm
     }
 
     /// Initialize surface volume from pressure and saturation by z = As.
-    /// Here  saturation is used as an intial guess for z in the
+    /// Here  saturation is used as an initial guess for z in the
     /// computation of A.
     template <class Props, class State>
     void initBlackoilSurfvol(const UnstructuredGrid& grid,
@@ -800,10 +800,6 @@ namespace Opm
                                         const Props& props,
                                         State& state)
     {
-
-        if (props.numPhases() != 3) {
-            OPM_THROW(std::runtime_error, "initBlackoilSurfvol() is only supported in three-phase simulations.");
-        }
         const std::vector<double>& rs = state.gasoilratio();
         const std::vector<double>& rv = state.rv();
 
@@ -828,15 +824,6 @@ namespace Opm
         std::vector<double> capPressures(number_of_cells*np);
         props.capPress(number_of_cells,&state.saturation()[0],&allcells[0],&capPressures[0],NULL);
 
-        std::vector<double> Pw(number_of_cells);
-        std::vector<double> Pg(number_of_cells);
-
-        for (int c = 0; c < number_of_cells; ++c){
-            Pw[c] = state.pressure()[c] + capPressures[c*np + BlackoilPhases::Aqua];
-            Pg[c] = state.pressure()[c] + capPressures[c*np + BlackoilPhases::Vapour];
-        }
-
-
         double z_tmp;
 
         // Water phase
diff --git a/opm/core/utility/initHydroCarbonState.hpp b/opm/core/utility/initHydroCarbonState.hpp
new file mode 100644
index 000000000..222946cf7
--- /dev/null
+++ b/opm/core/utility/initHydroCarbonState.hpp
@@ -0,0 +1,44 @@
+#ifndef INITHYDROCARBONSTATE_HPP
+#define INITHYDROCARBONSTATE_HPP
+
+#include "opm/core/simulator/BlackoilState.hpp"
+
+namespace Opm
+{
+
+void initHydroCarbonState(BlackoilState& state, const PhaseUsage& pu, const int num_cells, const bool has_disgas, const bool has_vapoil) {
+    enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour, Water = BlackoilPhases::Aqua };
+    // hydrocarbonstate is only used when gas and oil is present
+    assert(pu.phase_used[Oil]);
+    std::vector<HydroCarbonState>& hydroCarbonState = state.hydroCarbonState();
+    hydroCarbonState.resize(num_cells);
+    if (!pu.phase_used[Gas]) {
+        // hydroCarbonState should only be used when oil and gas is present. Return OilOnly to avoid potential trouble.
+        std::fill(hydroCarbonState.begin(), hydroCarbonState.end(), HydroCarbonState::OilOnly);
+        return;
+    }
+    const int np = pu.num_phases;
+    std::fill(hydroCarbonState.begin(), hydroCarbonState.end(), HydroCarbonState::GasAndOil);
+
+    // set hydrocarbon state
+    const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
+    const std::vector<double>& saturation = state.saturation();
+    for (int c = 0; c < num_cells; ++c) {
+        if (pu.phase_used[Water]) {
+            if ( saturation[c*np + pu.phase_pos[ Water ]] > (1.0 - epsilon)) {
+                continue; // cases (almost) filled with water is treated as GasAndOil case;
+            }
+        }
+        if ( saturation[c*np + pu.phase_pos[ Gas ]] == 0.0 && has_disgas) {
+            hydroCarbonState[c] = HydroCarbonState::OilOnly;
+            continue;
+        }
+        if ( saturation[c*np + pu.phase_pos[ Oil ]] == 0.0 && has_vapoil) {
+            hydroCarbonState[c] = HydroCarbonState::GasOnly;
+        }
+    }
+}
+
+
+} // namespace Opm
+#endif // INITHYDROCARBONSTATE_HPP
diff --git a/opm/core/wells/DynamicListEconLimited.hpp b/opm/core/wells/DynamicListEconLimited.hpp
new file mode 100644
index 000000000..22bf2355d
--- /dev/null
+++ b/opm/core/wells/DynamicListEconLimited.hpp
@@ -0,0 +1,92 @@
+/*
+  Copyright 2016 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_DYNAMICLISTECONLIMITED_HPP
+#define OPM_DYNAMICLISTECONLIMITED_HPP
+
+#include <vector>
+#include <string>
+#include <map>
+
+#include <cassert>
+
+namespace Opm
+{
+
+    /// to handle the wells and connections violating economic limits.
+    class DynamicListEconLimited
+    {
+    public:
+
+        DynamicListEconLimited() {
+        }
+
+        bool wellShutEconLimited(const std::string& well_name) const {
+            return std::find(m_shut_wells.begin(), m_shut_wells.end(), well_name) != m_shut_wells.end();
+        }
+
+        void addShutWell(const std::string& well_name) {
+            assert( !wellShutEconLimited(well_name) );
+            assert( !wellStoppedEconLimited(well_name) );
+
+            m_shut_wells.push_back(well_name);
+        }
+
+        bool wellStoppedEconLimited(const std::string& well_name) const {
+            return std::find(m_stopped_wells.begin(), m_stopped_wells.end(), well_name) != m_stopped_wells.end();
+        }
+
+        void addStoppedWell(const std::string& well_name) {
+            assert( !wellShutEconLimited(well_name) );
+            assert( !wellStoppedEconLimited(well_name) );
+
+            m_stopped_wells.push_back(well_name);
+        }
+
+
+        // TODO: maybe completion better here
+        bool anyConnectionClosedForWell(const std::string& well_name) const {
+            return (m_cells_closed_connections.find(well_name) != m_cells_closed_connections.end());
+        }
+
+        const std::vector<int>& getClosedConnectionsForWell(const std::string& well_name) const {
+            return (m_cells_closed_connections.find(well_name)->second);
+        }
+
+        void addClosedConnectionsForWell(const std::string& well_name,
+                                         const int cell_closed_connection) {
+            if (!anyConnectionClosedForWell(well_name)) {
+                // first time adding a connection for the well
+                std::vector<int> vector_cells = {cell_closed_connection};
+                m_cells_closed_connections[well_name] = vector_cells;
+            } else {
+                std::vector<int>& closed_connections = m_cells_closed_connections.find(well_name)->second;
+                closed_connections.push_back(cell_closed_connection);
+            }
+        }
+
+    private:
+        std::vector <std::string> m_shut_wells;
+        std::vector <std::string> m_stopped_wells;
+        // using grid cell number to indicate the location of the connections
+        std::map<std::string, std::vector<int>> m_cells_closed_connections;
+    };
+
+} // namespace Opm
+#endif  /* OPM_DYNAMICLISTECONLIMITED_HPP */
+
diff --git a/opm/core/wells/InjectionSpecification.cpp b/opm/core/wells/InjectionSpecification.cpp
index 342012467..09d11f60f 100644
--- a/opm/core/wells/InjectionSpecification.cpp
+++ b/opm/core/wells/InjectionSpecification.cpp
@@ -1,4 +1,26 @@
+/*
+  Copyright 2012  Sintef.
+  Copyright 2016  Statoil.
+
+  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 <stdexcept>
 #include "config.h"
+
+#include <opm/common/ErrorMacros.hpp>
 #include <opm/core/wells/InjectionSpecification.hpp>
 
 namespace Opm
@@ -18,4 +40,43 @@ namespace Opm
 
     }
 
+    std::string
+    InjectionSpecification::toString(const ControlMode& mode)
+    {
+        switch(mode) {
+        case ControlMode::NONE: return "NONE";
+        case ControlMode::RATE: return "RATE";
+        case ControlMode::RESV: return "RESV";
+        case ControlMode::BHP : return "BHP" ;
+        case ControlMode::THP : return "THP" ;
+        case ControlMode::REIN: return "REIN";
+        case ControlMode::VREP: return "VREP";
+        case ControlMode::GRUP: return "GRUP";
+        case ControlMode::FLD : return "FLD" ;
+        }
+        OPM_THROW(std::domain_error, "Unknown control mode " << mode << " encountered in injection specification");
+    }
+
+
+    std::string
+    InjectionSpecification::toString(const InjectorType& type)
+    {
+        switch(type) {
+        case InjectorType::WATER: return "WATER";
+        case InjectorType::OIL  : return "OIL"  ;
+        case InjectorType::GAS  : return "GAS"  ;
+        }
+        OPM_THROW(std::domain_error, "Unknown injector type " << type << " encountered in injection specification");
+    }
+
+
+    std::string 
+    InjectionSpecification::toString(const GuideRateType& type)
+    {
+        switch(type) {
+        case GuideRateType::RAT     : return "RAT"     ;
+        case GuideRateType::NONE_GRT: return "NONE_GRT";
+        }
+        OPM_THROW(std::domain_error, "Unknown guide rate type " << type << " encountered in injection specification");
+    }
 } // namespace Opm
diff --git a/opm/core/wells/InjectionSpecification.hpp b/opm/core/wells/InjectionSpecification.hpp
index 74d20f6a7..e47900b7d 100644
--- a/opm/core/wells/InjectionSpecification.hpp
+++ b/opm/core/wells/InjectionSpecification.hpp
@@ -2,6 +2,7 @@
 #define	OPM_INJECTORSPECIFICATION_HPP
 
 #include <opm/core/wells.h>
+#include <string>
 
 namespace Opm
 {
@@ -25,7 +26,9 @@ namespace Opm
         };
 
         InjectionSpecification();
-
+        static std::string toString(const ControlMode& mode);
+        static std::string toString(const InjectorType& type);
+        static std::string toString(const GuideRateType& type);
         InjectorType injector_type_;
         ControlMode control_mode_;
         double surface_flow_max_rate_;
diff --git a/opm/core/wells/ProductionSpecification.cpp b/opm/core/wells/ProductionSpecification.cpp
index 1c6d443a5..20abe7157 100644
--- a/opm/core/wells/ProductionSpecification.cpp
+++ b/opm/core/wells/ProductionSpecification.cpp
@@ -1,4 +1,26 @@
+/*
+  Copyright 2012  Sintef.
+  Copyright 2016  Statoil.
+
+  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 <stdexcept>
 #include "config.h"
+
+#include <opm/common/ErrorMacros.hpp>
 #include <opm/core/wells/ProductionSpecification.hpp>
 
 namespace Opm
@@ -19,4 +41,51 @@ namespace Opm
     {
     }
 
+
+    std::string
+    ProductionSpecification::toString(const ControlMode& mode)
+    {
+        switch(mode) {
+        case ControlMode::NONE: return "NONE";
+        case ControlMode::ORAT: return "ORAT";
+        case ControlMode::WRAT: return "WRAT";
+        case ControlMode::GRAT: return "GRAT";
+        case ControlMode::LRAT: return "LRAT";
+        case ControlMode::CRAT: return "CRAT";
+        case ControlMode::RESV: return "RESV";
+        case ControlMode::PRBL: return "RPBL";
+        case ControlMode::BHP : return "BHP" ;
+        case ControlMode::THP : return "THP" ;
+        case ControlMode::GRUP: return "GRUP";
+        case ControlMode::FLD : return "FLD" ;
+        }
+        OPM_THROW(std::domain_error, "Unknown control mode " << mode << " encountered in production specification");
+    }
+
+
+    std::string 
+    ProductionSpecification::toString(const Procedure& type)
+    {
+        switch(type) {
+        case Procedure::NONE_P: return "NONE_P";
+        case Procedure::RATE  : return "RATE"  ;
+        case Procedure::WELL  : return "WELL"  ;
+        }
+        OPM_THROW(std::domain_error, "Unknown procedure " << type << " encountered in production specification");
+    }
+
+
+    std::string 
+    ProductionSpecification::toString(const GuideRateType& type)
+    {
+        switch(type) {
+        case GuideRateType::OIL     : return "OIL"     ;
+        case GuideRateType::GAS     : return "GAS"     ;
+        case GuideRateType::WATER   : return "WATER"   ;
+        case GuideRateType::NONE_GRT: return "NONE_GRT";
+        }
+        OPM_THROW(std::domain_error, "Unknown guide rate type " << type << " encountered in production specification");
+    }
+
+
 }
diff --git a/opm/core/wells/ProductionSpecification.hpp b/opm/core/wells/ProductionSpecification.hpp
index 6d064a2f0..4bae4cae5 100644
--- a/opm/core/wells/ProductionSpecification.hpp
+++ b/opm/core/wells/ProductionSpecification.hpp
@@ -2,6 +2,7 @@
 #define	OPM_PRODUCTIONSPECIFICATION_HPP
 
 #include <opm/core/wells.h>
+#include <string>
 
 namespace Opm
 {
@@ -25,6 +26,9 @@ namespace Opm
         };
 
         ProductionSpecification();
+        static std::string toString(const ControlMode& mode);
+        static std::string toString(const Procedure& type);
+        static std::string toString(const GuideRateType& type);
 
         ControlMode control_mode_;
         Procedure procedure_;
diff --git a/opm/core/wells/WellCollection.cpp b/opm/core/wells/WellCollection.cpp
index ba94d0973..d1238e7c7 100644
--- a/opm/core/wells/WellCollection.cpp
+++ b/opm/core/wells/WellCollection.cpp
@@ -29,7 +29,7 @@
 
 namespace Opm
 {
-    void WellCollection::addField(GroupConstPtr fieldGroup, size_t timeStep, const PhaseUsage& phaseUsage) {
+    void WellCollection::addField(const Group* fieldGroup, size_t timeStep, const PhaseUsage& phaseUsage) {
         WellsGroupInterface* fieldNode = findNode(fieldGroup->name());
         if (fieldNode) {
             OPM_THROW(std::runtime_error, "Trying to add FIELD node, but this already exists. Can only have one FIELD node.");
@@ -38,7 +38,7 @@ namespace Opm
         roots_.push_back(createGroupWellsGroup(fieldGroup, timeStep, phaseUsage));
     }
 
-    void WellCollection::addGroup(GroupConstPtr groupChild, std::string parent_name,
+    void WellCollection::addGroup(const Group* groupChild, std::string parent_name,
                                   size_t timeStep, const PhaseUsage& phaseUsage) {
         WellsGroupInterface* parent = findNode(parent_name);
         if (!parent) {
@@ -60,7 +60,7 @@ namespace Opm
         child->setParent(parent);
     }
 
-    void WellCollection::addWell(WellConstPtr wellChild, size_t timeStep, const PhaseUsage& phaseUsage) {
+    void WellCollection::addWell(const Well* wellChild, size_t timeStep, const PhaseUsage& phaseUsage) {
         if (wellChild->getStatus(timeStep) == WellCommon::SHUT) {
             //SHUT wells are not added to the well collection
             return;
diff --git a/opm/core/wells/WellCollection.hpp b/opm/core/wells/WellCollection.hpp
index b597c4b64..b00827431 100644
--- a/opm/core/wells/WellCollection.hpp
+++ b/opm/core/wells/WellCollection.hpp
@@ -36,11 +36,11 @@ namespace Opm
     {
     public:
 
-        void addField(GroupConstPtr fieldGroup, size_t timeStep, const PhaseUsage& phaseUsage);
+        void addField(const Group* fieldGroup, size_t timeStep, const PhaseUsage& phaseUsage);
 
-        void addWell(WellConstPtr wellChild, size_t timeStep, const PhaseUsage& phaseUsage);
+        void addWell(const Well* wellChild, size_t timeStep, const PhaseUsage& phaseUsage);
 
-        void addGroup(GroupConstPtr groupChild, std::string parent_name,
+        void addGroup(const Group* groupChild, std::string parent_name,
                       size_t timeStep, const PhaseUsage& phaseUsage);
 
         /// Adds the child to the collection
diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp
index bf80e4640..558f39fc7 100644
--- a/opm/core/wells/WellsGroup.cpp
+++ b/opm/core/wells/WellsGroup.cpp
@@ -345,9 +345,10 @@ namespace Opm
                                             mode);
 
                 if (my_rate > target_rate) {
-                    std::cout << "Group " << mode<<" target not met for group " << name() << std::endl;
-                    std::cout << "target = " << target_rate << '\n'
-                              << "rate   = " << my_rate << std::endl;
+                    OpmLog::warning("Group " + InjectionSpecification::toString(mode) 
+                                    + " target not met for group " + name() + "\n"
+                                    + "target = " + std::to_string(target_rate) + "\n"
+                                    + "rate   = " + std::to_string(my_rate));
                     applyInjGroupControl(mode, target_rate, true);
                     injSpec().control_mode_ = mode;
                     return false;
@@ -377,9 +378,10 @@ namespace Opm
                                                   child_phases_summed.surf_prod_rates,
                                                   mode);
                 if (std::fabs(my_rate) > target_rate) {
-                    std::cout << "Group" << mode << " target not met for group " << name() << std::endl;
-                    std::cout << "target = " << target_rate << '\n'
-                              << "rate   = " << my_rate << std::endl;
+                    OpmLog::warning("Group" + ProductionSpecification::toString(mode) 
+                                    + " target not met for group " + name() + "\n"
+                                    + "target = " + std::to_string(target_rate) + '\n'
+                                    + "rate   = " + std::to_string(my_rate));
                     production_violated = true;
                     production_mode_violated = mode;
                     break;
@@ -677,9 +679,9 @@ namespace Opm
                 ctrl_violated = is_producer ? (my_target_bhp > my_well_bhp)
                     : (my_target_bhp < my_well_bhp);
                 if (ctrl_violated) {
-                    std::cout << "BHP limit violated for well " << name() << ":\n";
-                    std::cout << "BHP limit = " << my_target_bhp << std::endl;
-                    std::cout << "BHP       = " << my_well_bhp << std::endl;
+                    OpmLog::info("BHP limit violated for well " + name() + ":\n"
+                                 + "BHP limit = " + std::to_string(my_target_bhp)
+                                 + "BHP       = " + std::to_string(my_well_bhp));
                 }
                 break;
             }
@@ -698,9 +700,9 @@ namespace Opm
                 const double my_rate_target = well_controls_iget_target(ctrls , ctrl_index);
                 ctrl_violated = std::fabs(my_rate) - std::fabs(my_rate_target)> std::max(std::abs(my_rate), std::abs(my_rate_target))*1e-6;
                 if (ctrl_violated) {
-                    std::cout << "RESERVOIR_RATE limit violated for well " << name() << ":\n";
-                    std::cout << "rate limit = " << my_rate_target << std::endl;
-                    std::cout << "rate       = " << my_rate << std::endl;
+                    OpmLog::info("RESERVOIR_RATE limit violated for well " + name() + ":\n"
+                                 + "rate limit = " + std::to_string(my_rate_target)
+                                 + "rate       = " + std::to_string(my_rate));
                 }
                 break;
             }
@@ -714,9 +716,9 @@ namespace Opm
                 const double my_rate_target = well_controls_iget_target(ctrls , ctrl_index);
                 ctrl_violated = std::fabs(my_rate) > std::fabs(my_rate_target);
                 if (ctrl_violated) {
-                    std::cout << "SURFACE_RATE limit violated for well " << name() << ":\n";
-                    std::cout << "rate limit = " << my_rate_target << std::endl;
-                    std::cout << "rate       = " << my_rate << std::endl;
+                    OpmLog::info("SURFACE_RATE limit violated for well " + name() + ":\n"
+                                 + "rate limit = " + std::to_string(my_rate_target)
+                                 + "rate       = " + std::to_string(my_rate));
                 }
                 break;
             }
@@ -1062,7 +1064,7 @@ namespace Opm
         }
     } // anonymous namespace
 
-    std::shared_ptr<WellsGroupInterface> createGroupWellsGroup(GroupConstPtr group, size_t timeStep, const PhaseUsage& phase_usage )
+    std::shared_ptr<WellsGroupInterface> createGroupWellsGroup(const Group* group, size_t timeStep, const PhaseUsage& phase_usage )
     {
         InjectionSpecification injection_specification;
         ProductionSpecification production_specification;
@@ -1097,7 +1099,7 @@ namespace Opm
       'CMODE_UNDEFINED' - we do not carry that over the specification
       objects here.
      */
-    std::shared_ptr<WellsGroupInterface> createWellWellsGroup(WellConstPtr well, size_t timeStep, const PhaseUsage& phase_usage )
+    std::shared_ptr<WellsGroupInterface> createWellWellsGroup(const Well* well, size_t timeStep, const PhaseUsage& phase_usage )
     {
         InjectionSpecification injection_specification;
         ProductionSpecification production_specification;
diff --git a/opm/core/wells/WellsGroup.hpp b/opm/core/wells/WellsGroup.hpp
index 697a23cbc..c93fcdac4 100644
--- a/opm/core/wells/WellsGroup.hpp
+++ b/opm/core/wells/WellsGroup.hpp
@@ -406,14 +406,14 @@ namespace Opm
     /// \param[in] well the Well to construct object for
     /// \param[in] timeStep the time step in question
     /// \param[in] the phase usage
-    std::shared_ptr<WellsGroupInterface> createWellWellsGroup(WellConstPtr well, size_t timeStep,
+    std::shared_ptr<WellsGroupInterface> createWellWellsGroup(const Well* well, size_t timeStep,
                                                               const PhaseUsage& phase_usage );
 
     /// Creates the WellsGroupInterface for the given Group
     /// \param[in] group the Group to construct object for
     /// \param[in] timeStep the time step in question
     /// \param[in] the phase usage
-    std::shared_ptr<WellsGroupInterface> createGroupWellsGroup(GroupConstPtr group, size_t timeStep,
+    std::shared_ptr<WellsGroupInterface> createGroupWellsGroup(const Group* group, size_t timeStep,
                                                                const PhaseUsage& phase_usage );
 }
 #endif	/* OPM_WELLSGROUP_HPP */
diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp
index f9ad8569a..5a86c1139 100644
--- a/opm/core/wells/WellsManager.cpp
+++ b/opm/core/wells/WellsManager.cpp
@@ -333,11 +333,14 @@ namespace Opm
         : w_(0), is_parallel_run_(false)
     {
         std::vector<double> dummy_well_potentials;
+        // TODO: not sure about the usage of this WellsManager constructor
+        // TODO: not sure whether this is the correct thing to do here.
+        DynamicListEconLimited dummy_list_econ_limited;
         init(eclipseState, timeStep, UgGridHelpers::numCells(grid),
              UgGridHelpers::globalCell(grid), UgGridHelpers::cartDims(grid), 
              UgGridHelpers::dimensions(grid),
              UgGridHelpers::cell2Faces(grid), UgGridHelpers::beginFaceCentroids(grid),
-             permeability, dummy_well_potentials);
+             permeability, dummy_list_econ_limited, dummy_well_potentials);
 
     }
 
@@ -414,9 +417,10 @@ namespace Opm
 
 
 
-    void WellsManager::setupWellControls(std::vector<WellConstPtr>& wells, size_t timeStep,
+    void WellsManager::setupWellControls(std::vector< const Well* >& wells, size_t timeStep,
                                          std::vector<std::string>& well_names, const PhaseUsage& phaseUsage,
-                                         const std::vector<int>& wells_on_proc) {
+                                         const std::vector<int>& wells_on_proc,
+                                         const DynamicListEconLimited& list_econ_limited) {
         int well_index = 0;
         auto well_on_proc = wells_on_proc.begin();
 
@@ -427,18 +431,23 @@ namespace Opm
                 continue;
             }
 
-           WellConstPtr well = (*wellIter);
-
-            if (well->getStatus(timeStep) == WellCommon::STOP) {
-                // STOPed wells are kept in the well list but marked as stopped.
-                well_controls_stop_well(w_->ctrls[well_index]);
-            }
+            const auto* well = (*wellIter);
 
             if (well->getStatus(timeStep) == WellCommon::SHUT) {
                 //SHUT wells are not added to the well list
                 continue;
             }
 
+            if (list_econ_limited.wellShutEconLimited(well->name())) {
+                continue;
+            }
+
+            if (well->getStatus(timeStep) == WellCommon::STOP || list_econ_limited.wellStoppedEconLimited(well->name())) {
+                // Stopped wells are kept in the well list but marked as stopped.
+                well_controls_stop_well(w_->ctrls[well_index]);
+            }
+
+
             if (well->isInjector(timeStep)) {
                 const WellInjectionProperties& injectionProperties = well->getInjectionProperties(timeStep);
                 int ok = 1;
@@ -728,12 +737,12 @@ namespace Opm
     }
 
 
-    void WellsManager::setupGuideRates(std::vector<WellConstPtr>& wells, const size_t timeStep, std::vector<WellData>& well_data, std::map<std::string, int>& well_names_to_index,
+    void WellsManager::setupGuideRates(std::vector< const Well* >& wells, const size_t timeStep, std::vector<WellData>& well_data, std::map<std::string, int>& well_names_to_index,
                                        const PhaseUsage& phaseUsage, const std::vector<double>& well_potentials)
     {
         const int np = phaseUsage.num_phases;
         for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) {
-            WellConstPtr well = *wellIter;
+            const auto* well = *wellIter;
 
             if (well->getStatus(timeStep) == WellCommon::SHUT) {
                 //SHUT wells does not need guide rates
diff --git a/opm/core/wells/WellsManager.hpp b/opm/core/wells/WellsManager.hpp
index 022871abf..79d356bff 100644
--- a/opm/core/wells/WellsManager.hpp
+++ b/opm/core/wells/WellsManager.hpp
@@ -25,6 +25,7 @@
 
 #include <opm/core/wells/WellCollection.hpp>
 #include <opm/core/wells/WellsGroup.hpp>
+#include <opm/core/wells/DynamicListEconLimited.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/GroupTree.hpp>
 
 #include <opm/core/utility/CompressedPropertyAccess.hpp>
@@ -86,6 +87,7 @@ namespace Opm
                      const F2C& f2c,
                      FC begin_face_centroids,
                      const double* permeability,
+                     const DynamicListEconLimited& list_econ_limited,
                      bool is_parallel_run=false,
                      const std::vector<double>& well_potentials={});
 
@@ -155,17 +157,19 @@ namespace Opm
                   const C2F& cell_to_faces,
                   FC begin_face_centroids,
                   const double* permeability,
+                  const DynamicListEconLimited& list_econ_limited,
                   const std::vector<double>& well_potentials);
         // Disable copying and assignment.
         WellsManager(const WellsManager& other);
         WellsManager& operator=(const WellsManager& other);
         static void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed );
-        void setupWellControls(std::vector<WellConstPtr>& wells, size_t timeStep,
+        void setupWellControls(std::vector<const Well*>& wells, size_t timeStep,
                                std::vector<std::string>& well_names, const PhaseUsage& phaseUsage,
-                               const std::vector<int>& wells_on_proc);
+                               const std::vector<int>& wells_on_proc,
+                               const DynamicListEconLimited& list_econ_limited);
 
         template<class C2F, class FC, class NTG>
-        void createWellsFromSpecs( std::vector<WellConstPtr>& wells, size_t timeStep,
+        void createWellsFromSpecs( std::vector<const Well*>& wells, size_t timeStep,
                                    const C2F& cell_to_faces, 
                                    const int* cart_dims,
                                    FC begin_face_centroids, 
@@ -178,10 +182,11 @@ namespace Opm
                                    const std::map<int,int>& cartesian_to_compressed,
                                    const double* permeability,
                                    const NTG& ntg,
-                                   std::vector<int>& wells_on_proc);
+                                   std::vector<int>& wells_on_proc,
+                                   const DynamicListEconLimited& list_econ_limited);
 
         void addChildGroups(GroupTreeNodeConstPtr parentNode, std::shared_ptr< const Schedule > schedule, size_t timeStep, const PhaseUsage& phaseUsage);
-        void setupGuideRates(std::vector<WellConstPtr>& wells, const size_t timeStep, std::vector<WellData>& well_data, std::map<std::string, int>& well_names_to_index,
+        void setupGuideRates(std::vector<const Well*>& wells, const size_t timeStep, std::vector<WellData>& well_data, std::map<std::string, int>& well_names_to_index,
                              const PhaseUsage& phaseUsage, const std::vector<double>& well_potentials);
         // Data
         Wells* w_;
diff --git a/opm/core/wells/WellsManager_impl.hpp b/opm/core/wells/WellsManager_impl.hpp
index 1f44443a0..c637eeff4 100644
--- a/opm/core/wells/WellsManager_impl.hpp
+++ b/opm/core/wells/WellsManager_impl.hpp
@@ -2,6 +2,7 @@
 #include <opm/core/grid/GridHelpers.hpp>
 
 #include <opm/common/ErrorMacros.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
 #include <opm/core/utility/compressedToCartesian.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/CompletionSet.hpp>
@@ -106,7 +107,7 @@ getCubeDim(const C2F& c2f,
 namespace Opm
 {
 template<class C2F, class FC, class NTG>
-void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t timeStep,
+void WellsManager::createWellsFromSpecs(std::vector<const Well*>& wells, size_t timeStep,
                                         const C2F& c2f,
                                         const int* cart_dims,
                                         FC begin_face_centroids,
@@ -119,7 +120,8 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
                                         const std::map<int,int>& cartesian_to_compressed,
                                         const double* permeability,
                                         const NTG& ntg,
-                                        std::vector<int>& wells_on_proc)
+                                        std::vector<int>& wells_on_proc,
+                                        const DynamicListEconLimited& list_econ_limited)
 {
     if (dimensions != 3) {
         OPM_THROW(std::domain_error,
@@ -137,12 +139,21 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
     // the index of the well according to the eclipse state
     int well_index_on_proc = 0;
     for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
-        WellConstPtr well = (*wellIter);
+        const auto* well = (*wellIter);
 
         if (well->getStatus(timeStep) == WellCommon::SHUT) {
             continue;
         }
 
+        if (list_econ_limited.wellShutEconLimited(well->name())) {
+            continue;
+        }
+
+        std::vector<int> cells_connection_closed;
+        if (list_econ_limited.anyConnectionClosedForWell(well->name())) {
+            cells_connection_closed = list_econ_limited.getClosedConnectionsForWell(well->name());
+        }
+
         {   // COMPDAT handling
             auto completionSet = well->getCompletions(timeStep);
             // shut completions and open ones stored in this process will have 1 others 0.
@@ -174,6 +185,16 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
                     else
                     {
                         int cell = cgit->second;
+                        // check if the connection is closed due to economic limits
+                        if (!cells_connection_closed.empty()) {
+                            const bool connection_found = std::find(cells_connection_closed.begin(),
+                                                                    cells_connection_closed.end(), cell)
+                                                          != cells_connection_closed.end();
+                            if (connection_found) {
+                                continue;
+                            }
+                        }
+
                         PerfData pd;
                         pd.cell = cell;
                         {
@@ -233,9 +254,9 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
                     // Check that the complete well is on this process
                     if ( sum_completions_on_proc < completionSet->size() )
                     {
-                        std::cout<< "Well "<< well->name() << " does not seem to be"
-                                 << "completely in the disjoint partition of "
-                                 << "process. Therefore we  deactivate it here." << std::endl;
+                        OpmLog::warning("Well " + well->name() + " does not seem to be"
+                                        + "completely in the disjoint partition of "
+                                        + "process. Therefore we  deactivate it here.");
                         // Mark well as not existent on this process
                         wells_on_proc[wellIter-wells.begin()] = 0;
                         wellperf_data[well_index_on_proc].clear();
@@ -326,13 +347,14 @@ WellsManager(const Opm::EclipseStateConstPtr eclipseState,
              const C2F&                      cell_to_faces,
              FC                              begin_face_centroids,
              const double*                   permeability,
+             const DynamicListEconLimited&   list_econ_limited,
              bool                            is_parallel_run,
              const std::vector<double>&      well_potentials)
     : w_(0), is_parallel_run_(is_parallel_run)
 {
     init(eclipseState, timeStep, number_of_cells, global_cell,
          cart_dims, dimensions,
-         cell_to_faces, begin_face_centroids, permeability, well_potentials);
+         cell_to_faces, begin_face_centroids, permeability, list_econ_limited, well_potentials);
 }
 
 /// Construct wells from deck.
@@ -347,6 +369,7 @@ WellsManager::init(const Opm::EclipseStateConstPtr eclipseState,
                    const C2F&                      cell_to_faces,
                    FC                              begin_face_centroids,
                    const double*                   permeability,
+                   const DynamicListEconLimited&   list_econ_limited,
                    const std::vector<double>&      well_potentials)
 {
     if (dimensions != 3) {
@@ -378,8 +401,8 @@ WellsManager::init(const Opm::EclipseStateConstPtr eclipseState,
     std::map<std::string, int> well_names_to_index;
 
     auto schedule = eclipseState->getSchedule();
-    std::vector<WellConstPtr> wells    = schedule->getWells(timeStep);
-    std::vector<int>          wells_on_proc;
+    auto wells       = schedule->getWells(timeStep);
+    std::vector<int> wells_on_proc;
 
     well_names.reserve(wells.size());
     well_data.reserve(wells.size());
@@ -390,7 +413,7 @@ WellsManager::init(const Opm::EclipseStateConstPtr eclipseState,
     DoubleArray ntg_glob(eclipseState, "NTG", 1.0);
     NTGArray    ntg(ntg_glob, global_cell);
 
-    EclipseGridConstPtr eclGrid = eclipseState->getEclipseGrid();
+    EclipseGridConstPtr eclGrid = eclipseState->getInputGrid();
 
     // use cell thickness (dz) from eclGrid
     // dz overwrites values calculated by WellDetails::getCubeDim
@@ -409,16 +432,15 @@ WellsManager::init(const Opm::EclipseStateConstPtr eclipseState,
                          dz,
                          well_names, well_data, well_names_to_index,
                          pu, cartesian_to_compressed, permeability, ntg,
-                         wells_on_proc);
+                         wells_on_proc, list_econ_limited);
 
-    setupWellControls(wells, timeStep, well_names, pu, wells_on_proc);
+    setupWellControls(wells, timeStep, well_names, pu, wells_on_proc, list_econ_limited);
 
     {
         GroupTreeNodeConstPtr fieldNode =
             schedule->getGroupTree(timeStep)->getNode("FIELD");
 
-        GroupConstPtr fieldGroup =
-            schedule->getGroup(fieldNode->name());
+        const auto* fieldGroup = schedule->getGroup(fieldNode->name());
 
         well_collection_.addField(fieldGroup, timeStep, pu);
         addChildGroups(fieldNode, schedule, timeStep, pu);
diff --git a/opm/core/wells/wells.c b/opm/core/wells/wells.c
index 026b470e3..75827697f 100644
--- a/opm/core/wells/wells.c
+++ b/opm/core/wells/wells.c
@@ -545,6 +545,16 @@ bool
 wells_equal(const struct Wells *W1, const struct Wells *W2 , bool verbose)
 /* ---------------------------------------------------------------------- */
 {
+    // Cater the case where W1 and W2 are the same (null) pointers.
+    if( W1 == W2 )
+    {
+        return true;
+    }
+    if( W1 == NULL || W2 == NULL)
+    {
+        return false;
+    }
+    
     bool are_equal = true;
     are_equal = (W1->number_of_wells == W2->number_of_wells);
     are_equal = are_equal && (W1->number_of_phases == W2->number_of_phases);
diff --git a/tests/norne_pvt.data b/tests/norne_pvt.data
index bdb7af927..55192cdd6 100644
--- a/tests/norne_pvt.data
+++ b/tests/norne_pvt.data
@@ -15,6 +15,14 @@ DIMENS
 1 1 1 /
 
 GRID
+DX
+1 /
+DY
+1 /
+DZ
+1 /
+TOPS
+1 /
 
 PROPS
 
diff --git a/tests/test_blackoilstate.cpp b/tests/test_blackoilstate.cpp
index 0cff885a9..25d41c82c 100644
--- a/tests/test_blackoilstate.cpp
+++ b/tests/test_blackoilstate.cpp
@@ -28,25 +28,36 @@ using namespace Opm;
 using namespace std;
 
 
+// ACTNUM 1 998*2 3
+std::vector<int> get_testBlackoilStateActnum() {
+    std::vector<int> actnum(10 * 10 * 10, 2);
+    actnum.front() = 1;
+    actnum.back() = 3;
+    return actnum;
+}
 
 BOOST_AUTO_TEST_CASE(EqualsDifferentDeckReturnFalse) {
-
-    ParseContext parseContext;
     const string filename1 = "testBlackoilState1.DATA";
     const string filename2 = "testBlackoilState2.DATA";
-    Opm::ParserPtr parser(new Opm::Parser());
-    Opm::DeckConstPtr deck1(parser->parseFile(filename1, parseContext));
-    Opm::DeckConstPtr deck2(parser->parseFile(filename2, parseContext));
 
-    GridManager gridManager1(deck1);
+    const auto es1 = Opm::Parser::parse(filename1);
+    auto eg1 = es1.getInputGridCopy();
+    std::vector<int> actnum = get_testBlackoilStateActnum();
+    eg1->resetACTNUM(actnum.data());
+
+    const auto es2 = Opm::Parser::parse(filename2);
+    auto eg2 = es2.getInputGrid();
+
+    GridManager gridManager1(eg1);
     const UnstructuredGrid& grid1 = *gridManager1.c_grid();
-    GridManager gridManager2(deck2);
+
+    GridManager gridManager2(eg2);
     const UnstructuredGrid& grid2 = *gridManager2.c_grid();
 
     BlackoilState state1( UgGridHelpers::numCells( grid1 ) , UgGridHelpers::numFaces( grid1 ) , 3);
     BlackoilState state2( UgGridHelpers::numCells( grid2 ) , UgGridHelpers::numFaces( grid2 ) , 3);
 
-    BOOST_CHECK_EQUAL( false , state1.equal(state2) );
+    BOOST_CHECK( ! state1.equal(state2) );
 }
 
 
@@ -55,44 +66,47 @@ BOOST_AUTO_TEST_CASE(EqualsDifferentDeckReturnFalse) {
 BOOST_AUTO_TEST_CASE(EqualsNumericalDifferenceReturnFalse) {
 
     const string filename = "testBlackoilState1.DATA";
-    Opm::ParseContext parseContext;
-    Opm::ParserPtr parser(new Opm::Parser());
-    Opm::DeckConstPtr deck(parser->parseFile(filename , parseContext));
 
-    GridManager gridManager(deck);
+    const auto es = Opm::Parser::parse(filename);
+    auto eg = es.getInputGridCopy();
+
+    std::vector<int> actnum = get_testBlackoilStateActnum();
+    eg->resetACTNUM(actnum.data());
+
+    GridManager gridManager(eg);
     const UnstructuredGrid& grid = *gridManager.c_grid();
 
     BlackoilState state1( UgGridHelpers::numCells( grid ) , UgGridHelpers::numFaces( grid ) , 3);
     BlackoilState state2( UgGridHelpers::numCells( grid ) , UgGridHelpers::numFaces( grid ) , 3);
 
 
-    BOOST_CHECK_EQUAL( true , state1.equal(state2) );
+    BOOST_CHECK( state1.equal(state2) );
     {
         std::vector<double>& p1 = state1.pressure();
         std::vector<double>& p2 = state2.pressure();
         p1[0] = p1[0] * 2 + 1;
 
-        BOOST_CHECK_EQUAL( false , state1.equal(state2) );
+        BOOST_CHECK( ! state1.equal(state2) );
         p1[0] = p2[0];
-        BOOST_CHECK_EQUAL( true , state1.equal(state2) );
+        BOOST_CHECK(   state1.equal(state2) );
     }
     {
         std::vector<double>& gor1 = state1.gasoilratio();
         std::vector<double>& gor2 = state2.gasoilratio();
         gor1[0] = gor1[0] * 2 + 1;
 
-        BOOST_CHECK_EQUAL( false , state1.equal(state2) );
+        BOOST_CHECK( ! state1.equal(state2) );
         gor1[0] = gor2[0];
-        BOOST_CHECK_EQUAL( true , state1.equal(state2) );
+        BOOST_CHECK(   state1.equal(state2) );
     }
     {
         std::vector<double>& p1 = state1.facepressure();
         std::vector<double>& p2 = state2.facepressure();
         p1[0] = p1[0] * 2 + 1;
 
-        BOOST_CHECK_EQUAL( false , state1.equal(state2) );
+        BOOST_CHECK( ! state1.equal(state2) );
         p1[0] = p2[0];
-        BOOST_CHECK_EQUAL( true , state1.equal(state2) );
+        BOOST_CHECK(   state1.equal(state2) );
     }
 
     {
@@ -101,9 +115,9 @@ BOOST_AUTO_TEST_CASE(EqualsNumericalDifferenceReturnFalse) {
         if (f1.size() > 0 ) {
             f1[0] = f1[0] * 2 + 1;
 
-            BOOST_CHECK_EQUAL( false , state1.equal(state2) );
+            BOOST_CHECK( ! state1.equal(state2) );
             f1[0] = f2[0];
-            BOOST_CHECK_EQUAL( true , state1.equal(state2) );
+            BOOST_CHECK(   state1.equal(state2) );
         }
     }
     {
@@ -112,9 +126,9 @@ BOOST_AUTO_TEST_CASE(EqualsNumericalDifferenceReturnFalse) {
         if (sv1.size() > 0) {
             sv1[0] = sv1[0] * 2 + 1;
 
-            BOOST_CHECK_EQUAL( false , state1.equal(state2) );
+            BOOST_CHECK( ! state1.equal(state2) );
             sv1[0] = sv2[0];
-            BOOST_CHECK_EQUAL( true , state1.equal(state2) );
+            BOOST_CHECK(   state1.equal(state2) );
         }
     }
     {
@@ -122,8 +136,8 @@ BOOST_AUTO_TEST_CASE(EqualsNumericalDifferenceReturnFalse) {
         std::vector<double>& sat2 = state2.saturation();
         sat1[0] = sat1[0] * 2 + 1;
 
-        BOOST_CHECK_EQUAL( false , state1.equal(state2) );
+        BOOST_CHECK( ! state1.equal(state2) );
         sat1[0] = sat2[0];
-        BOOST_CHECK_EQUAL( true , state1.equal(state2) );
+        BOOST_CHECK(   state1.equal(state2) );
     }
 }
diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp
index ed7f74ffa..caedd0af9 100644
--- a/tests/test_equil.cpp
+++ b/tests/test_equil.cpp
@@ -366,7 +366,7 @@ BOOST_AUTO_TEST_CASE (DeckAllDead)
     Opm::ParseContext parseContext;
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::DeckConstPtr deck = parser->parseFile("deadfluids.DATA" , parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck, parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
     Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, *grid, false);
     Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, *grid, 10.0);
     const auto& pressures = comp.press();
@@ -394,7 +394,7 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::ParseContext parseContext;
     Opm::DeckConstPtr deck = parser->parseFile("capillary.DATA" , parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
     Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
 
     // Test the capillary inversion for oil-water.
@@ -448,7 +448,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary)
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::ParseContext parseContext;
     Opm::DeckConstPtr deck = parser->parseFile("capillary.DATA" , parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
     Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
 
     Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 10.0);
@@ -489,7 +489,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::ParseContext parseContext;
     Opm::DeckConstPtr deck = parser->parseFile("capillary_overlap.DATA" , parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
     Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
 
     Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665);
@@ -552,7 +552,7 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::ParseContext parseContext;
     Opm::DeckConstPtr deck = parser->parseFile("equil_liveoil.DATA" , parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
     Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
 
     Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665);
@@ -632,7 +632,7 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::ParseContext parseContext;
     Opm::DeckConstPtr deck = parser->parseFile("equil_livegas.DATA" , parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
     Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
 
     Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665);
@@ -715,7 +715,7 @@ BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::ParseContext parseContext;
     Opm::DeckConstPtr deck = parser->parseFile("equil_rsvd_and_rvvd.DATA", parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
     Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
 
     Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665);
diff --git a/tests/test_norne_pvt.cpp b/tests/test_norne_pvt.cpp
index 664698fb2..4fb5b6e81 100644
--- a/tests/test_norne_pvt.cpp
+++ b/tests/test_norne_pvt.cpp
@@ -281,7 +281,7 @@ BOOST_AUTO_TEST_CASE( Test_Norne_PVT) {
     std::shared_ptr<const Deck> deck;
     deck = parser->parseFile("norne_pvt.data", parseContext);
 
-    Opm::EclipseStateConstPtr eclState(new EclipseState(deck, parseContext));
+    Opm::EclipseStateConstPtr eclState(new EclipseState(*deck, parseContext));
 
     verify_norne_oil_pvt_region1( deck, eclState );
     verify_norne_oil_pvt_region2( deck, eclState );
diff --git a/tests/test_parallelistlinformation.cpp b/tests/test_parallelistlinformation.cpp
index 6bd927008..f0a07283e 100644
--- a/tests/test_parallelistlinformation.cpp
+++ b/tests/test_parallelistlinformation.cpp
@@ -33,6 +33,7 @@
 #include <functional>
 #ifdef HAVE_DUNE_ISTL
 
+
 template<typename T>
 void runSumMaxMinTest(const T offset)
 {
@@ -45,12 +46,13 @@ void runSumMaxMinTest(const T offset)
     assert(comm.indexSet()->size()==x.size());
     for(auto it=comm.indexSet()->begin(), itend=comm.indexSet()->end(); it!=itend; ++it)
         x[it->local()]=it->global()+offset;
-    auto containers = std::make_tuple(x, x, x, x);
+    auto containers = std::make_tuple(x, x, x, x, x);
     auto operators  = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<T>(),
                                       Opm::Reduction::makeGlobalMaxFunctor<T>(),
                                       Opm::Reduction::makeGlobalMinFunctor<T>(),
-                                      Opm::Reduction::makeInnerProductFunctor<T>());
-    auto values     = std::tuple<T,T,T,T>(0,0,100000, 0);
+                                      Opm::Reduction::makeInnerProductFunctor<T>(),
+                                      Opm::Reduction::makeLInfinityNormFunctor<T>());
+    auto values     = std::tuple<T,T,T,T,T>(0,0,100000, 0, 0);
     auto oldvalues  = values;
     start = offset;
     end   = start+N;
@@ -59,10 +61,14 @@ void runSumMaxMinTest(const T offset)
     BOOST_CHECK(std::get<1>(values)==std::max(N+offset-1, std::get<1>(oldvalues)));
     BOOST_CHECK(std::get<2>(values)==std::min(offset, std::get<2>(oldvalues)));
     BOOST_CHECK(std::get<3>(values)==((end-1)*end*(2*end-1)-(start-1)*start*(2*start-1))/6+std::get<3>(oldvalues));
+    // Must avoid std::abs() directly to prevent ambiguity with unsigned integers.
+    Opm::Reduction::detail::MaxAbsFunctor<T> maxabsfunc;
+    BOOST_CHECK(std::get<4>(values)==maxabsfunc(offset, N+offset-1));
 }
 
 BOOST_AUTO_TEST_CASE(tupleReductionTestInt)
 {
+    runSumMaxMinTest<int>(-200);
     runSumMaxMinTest<int>(0);
     runSumMaxMinTest<int>(20);
     runSumMaxMinTest<int>(-20);
@@ -75,6 +81,7 @@ BOOST_AUTO_TEST_CASE(tupleReductionTestUnsignedInt)
 }
 BOOST_AUTO_TEST_CASE(tupleReductionTestFloat)
 {
+    runSumMaxMinTest<float>(-200);
     runSumMaxMinTest<float>(0);
     runSumMaxMinTest<float>(20);
     runSumMaxMinTest<float>(-20);
diff --git a/tests/test_pinchprocessor.cpp b/tests/test_pinchprocessor.cpp
index 6a2711a5d..78659b8fe 100644
--- a/tests/test_pinchprocessor.cpp
+++ b/tests/test_pinchprocessor.cpp
@@ -54,9 +54,9 @@ BOOST_AUTO_TEST_CASE(Processing)
     Opm::ParserPtr parser(new Opm::Parser());
     Opm::ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }});
     Opm::DeckConstPtr deck = parser->parseFile(filename, parseContext);
-    std::shared_ptr<EclipseState> eclstate (new Opm::EclipseState(deck, parseContext));
+    std::shared_ptr<EclipseState> eclstate (new Opm::EclipseState(*deck, parseContext));
     const auto& porv = eclstate->get3DProperties().getDoubleGridProperty("PORV").getData();
-    EclipseGridConstPtr eclgrid = eclstate->getEclipseGrid();
+    EclipseGridConstPtr eclgrid = eclstate->getInputGrid();
 
     BOOST_CHECK_EQUAL(eclgrid->getMinpvMode(), MinpvMode::EclSTD);
 
@@ -92,10 +92,10 @@ BOOST_AUTO_TEST_CASE(Processing)
     std::vector<double> htrans(Opm::UgGridHelpers::numCellFaces(grid));
     Grid* ug = const_cast<Grid*>(& grid);
     tpfa_htrans_compute(ug, rock.permeability(), htrans.data());
-    auto transMult = eclstate->getTransMult();
+    const auto& transMult = eclstate->getTransMult();
     std::vector<double> multz(nc, 0.0);
     for (int i = 0; i < nc; ++i) {
-        multz[i] = transMult->getMultiplier(global_cell[i], Opm::FaceDir::ZPlus);
+        multz[i] = transMult.getMultiplier(global_cell[i], Opm::FaceDir::ZPlus);
     }
     Opm::NNC nnc(deck, eclgrid);
     pinch.process(grid, htrans, actnum, multz, porv, nnc);
diff --git a/tests/test_relpermdiagnostics.cpp b/tests/test_relpermdiagnostics.cpp
index 65862ae3d..814fc073c 100644
--- a/tests/test_relpermdiagnostics.cpp
+++ b/tests/test_relpermdiagnostics.cpp
@@ -34,7 +34,7 @@
 #include <boost/test/floating_point_comparison.hpp>
 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
 #include <opm/common/OpmLog/OpmLog.hpp>
-#include <opm/common/OpmLog/EclipsePRTLog.hpp>
+#include <opm/common/OpmLog/CounterLog.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/grid/cart_grid.h>
@@ -57,17 +57,13 @@ BOOST_AUTO_TEST_CASE(diagnosis)
                               { ParseContext::PARSE_RANDOM_TEXT, InputError::IGNORE}
                              });
     Opm::DeckConstPtr deck(parser->parseFile("../tests/relpermDiagnostics.DATA", parseContext));
-    eclState.reset(new EclipseState(deck, parseContext));
-
-    GridManager gm(deck);
+    eclState.reset(new EclipseState(*deck, parseContext));
+    GridManager gm(eclState->getInputGrid());
     const UnstructuredGrid& grid = *gm.c_grid();
-    std::string logFile = "LOGFILE.txt";
-    std::shared_ptr<EclipsePRTLog> prtLog = std::make_shared<EclipsePRTLog>(logFile, Log::DefaultMessageTypes);
-    OpmLog::addBackend( "ECLIPSEPRTLOG" , prtLog );
+    std::shared_ptr<CounterLog> counterLog = std::make_shared<CounterLog>(Log::DefaultMessageTypes);
+    OpmLog::addBackend( "COUNTERLOG" , counterLog );
     RelpermDiagnostics diagnostics;
     diagnostics.diagnosis(eclState, deck, grid);
-    auto msg = diagnostics.getMessages();
-    BOOST_CHECK(!msg.empty());
-    BOOST_CHECK_EQUAL(msg.size(), 1);
+    BOOST_CHECK_EQUAL(1, counterLog->numMessages(Log::MessageType::Warning));
 }
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/tests/test_satfunc.cpp b/tests/test_satfunc.cpp
index b5630589a..898acf12a 100644
--- a/tests/test_satfunc.cpp
+++ b/tests/test_satfunc.cpp
@@ -69,7 +69,7 @@ BOOST_AUTO_TEST_CASE (GwsegStandard)
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::ParseContext parseContext;
     Opm::DeckConstPtr deck = parser->parseFile("satfuncStandard.DATA", parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
     Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
     
     const int np = 3;
@@ -155,7 +155,7 @@ BOOST_AUTO_TEST_CASE (GwsegEPSBase)
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::ParseContext parseContext;
     Opm::DeckConstPtr deck = parser->parseFile("satfuncEPSBase.DATA" , parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
     Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
     
     const int np = 3;
@@ -241,7 +241,7 @@ BOOST_AUTO_TEST_CASE (GwsegEPS_A)
     Opm::ParseContext parseContext;
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::DeckConstPtr deck = parser->parseFile("satfuncEPS_A.DATA" , parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
     Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
     
     const int np = 3;
@@ -493,7 +493,7 @@ BOOST_AUTO_TEST_CASE (GwsegEPS_C)
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::ParseContext parseContext;
     Opm::DeckConstPtr deck = parser->parseFile("satfuncEPS_C.DATA", parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
     Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
     
     const int np = 3;
@@ -596,7 +596,7 @@ BOOST_AUTO_TEST_CASE (GwsegEPS_D)
     Opm::ParserPtr parser(new Opm::Parser() );
     Opm::ParseContext parseContext;
     Opm::DeckConstPtr deck = parser->parseFile("satfuncEPS_D.DATA" , parseContext);
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
     Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
     
     const int np = 3;
diff --git a/tests/test_stoppedwells.cpp b/tests/test_stoppedwells.cpp
index 213e59857..219377d7f 100644
--- a/tests/test_stoppedwells.cpp
+++ b/tests/test_stoppedwells.cpp
@@ -47,8 +47,8 @@ BOOST_AUTO_TEST_CASE(TestStoppedWells)
     Opm::ParseContext parseContext;
     Opm::DeckConstPtr deck(parser->parseFile(filename , parseContext));
 
-    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
-    Opm::GridManager gridManager(deck);
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
+    Opm::GridManager gridManager(eclipseState->getInputGrid());
 
     double target_surfacerate_inj;
     double target_surfacerate_prod;
diff --git a/tests/test_wellcollection.cpp b/tests/test_wellcollection.cpp
index 6870853d3..35860e9a2 100644
--- a/tests/test_wellcollection.cpp
+++ b/tests/test_wellcollection.cpp
@@ -44,7 +44,7 @@ BOOST_AUTO_TEST_CASE(AddWellsAndGroupToCollection) {
     std::string scheduleFile("wells_group.data");
     ParseContext parseContext;
     DeckConstPtr deck =  parser->parseFile(scheduleFile, parseContext);
-    EclipseStateConstPtr eclipseState(new EclipseState(deck, parseContext));
+    EclipseStateConstPtr eclipseState(new EclipseState(*deck, parseContext));
     PhaseUsage pu = phaseUsageFromDeck(eclipseState);
 
     GroupTreeNodePtr field=eclipseState->getSchedule()->getGroupTree(2)->getNode("FIELD");
@@ -54,24 +54,24 @@ BOOST_AUTO_TEST_CASE(AddWellsAndGroupToCollection) {
     WellCollection collection;
 
     // Add groups to WellCollection
-    GroupConstPtr fieldGroup =  eclipseState->getSchedule()->getGroup(field->name());
+    const auto* fieldGroup =  eclipseState->getSchedule()->getGroup(field->name());
     collection.addField(fieldGroup, 2, pu);
 
     for (auto iter = field->begin(); iter != field->end(); ++iter) {
-        GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name());
+        const auto* childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name());
         collection.addGroup(childGroupNode, fieldGroup->name(), 2, pu);
     }
 
-    GroupConstPtr g1Group =  eclipseState->getSchedule()->getGroup(g1->name());
+    const auto* g1Group = eclipseState->getSchedule()->getGroup(g1->name());
     for (auto iter = g1->begin(); iter != g1->end(); ++iter) {
-        GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name());
+        const auto* childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name());
         collection.addGroup(childGroupNode, g1Group->name(), 2, pu);
     }
 
 
-    GroupConstPtr g2Group =  eclipseState->getSchedule()->getGroup(g2->name());
+    const auto* g2Group =  eclipseState->getSchedule()->getGroup(g2->name());
     for (auto iter = g2->begin(); iter != g2->end(); ++iter) {
-        GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name());
+        const auto* childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name());
         collection.addGroup(childGroupNode, g2Group->name(), 2, pu);
     }
 
@@ -81,7 +81,7 @@ BOOST_AUTO_TEST_CASE(AddWellsAndGroupToCollection) {
 
     // Add wells to WellCollection
     WellCollection wellCollection;
-    std::vector<WellConstPtr> wells = eclipseState->getSchedule()->getWells();
+    auto wells = eclipseState->getSchedule()->getWells();
     for (size_t i=0; i<wells.size(); i++) {
         collection.addWell(wells[i], 2, pu);
     }
diff --git a/tests/test_wellsgroup.cpp b/tests/test_wellsgroup.cpp
index 0556bc71e..f6e24cd41 100644
--- a/tests/test_wellsgroup.cpp
+++ b/tests/test_wellsgroup.cpp
@@ -52,13 +52,13 @@ BOOST_AUTO_TEST_CASE(ConstructGroupFromWell) {
     std::string scheduleFile("wells_group.data");
     ParseContext parseContext;
     DeckConstPtr deck =  parser->parseFile(scheduleFile, parseContext);
-    EclipseStateConstPtr eclipseState(new EclipseState(deck , parseContext));
+    EclipseStateConstPtr eclipseState(new EclipseState(*deck , parseContext));
     PhaseUsage pu = phaseUsageFromDeck(eclipseState);
 
-    std::vector<WellConstPtr> wells = eclipseState->getSchedule()->getWells();
+    auto wells = eclipseState->getSchedule()->getWells();
 
     for (size_t i=0; i<wells.size(); i++) {
-        WellConstPtr well = wells[i];
+        const auto* well = wells[i];
         std::shared_ptr<WellsGroupInterface> wellsGroup = createWellWellsGroup(well, 2, pu);
         BOOST_CHECK_EQUAL(well->name(), wellsGroup->name());
         if (well->isInjector(2)) {
@@ -85,13 +85,13 @@ BOOST_AUTO_TEST_CASE(ConstructGroupFromGroup) {
     ParseContext parseContext;
     std::string scheduleFile("wells_group.data");
     DeckConstPtr deck =  parser->parseFile(scheduleFile, parseContext);
-    EclipseStateConstPtr eclipseState(new EclipseState(deck , parseContext));
+    EclipseStateConstPtr eclipseState(new EclipseState(*deck , parseContext));
     PhaseUsage pu = phaseUsageFromDeck(eclipseState);
 
     std::vector<GroupTreeNodeConstPtr> nodes = eclipseState->getSchedule()->getGroupTree(2)->getNodes();
 
     for (size_t i=0; i<nodes.size(); i++) {
-        GroupConstPtr group = eclipseState->getSchedule()->getGroup(nodes[i]->name());
+        const auto* group = eclipseState->getSchedule()->getGroup(nodes[i]->name());
         std::shared_ptr<WellsGroupInterface> wellsGroup = createGroupWellsGroup(group, 2, pu);
         BOOST_CHECK_EQUAL(group->name(), wellsGroup->name());
         if (group->isInjectionGroup(2)) {
diff --git a/tests/test_wellsmanager.cpp b/tests/test_wellsmanager.cpp
index da5995f20..cb8cd7ad0 100644
--- a/tests/test_wellsmanager.cpp
+++ b/tests/test_wellsmanager.cpp
@@ -29,273 +29,255 @@
 #include <boost/test/unit_test.hpp>
 
 #include <opm/parser/eclipse/Parser/Parser.hpp>
-        #include <opm/parser/eclipse/Parser/ParseContext.hpp>
-        #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
 
-        #include <opm/core/wells/WellsManager.hpp>
-        #include <opm/core/wells.h>
-        #include <opm/core/well_controls.h>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/wells.h>
+#include <opm/core/well_controls.h>
 
-        #include <opm/core/grid/GridManager.hpp>
+#include <opm/core/grid/GridManager.hpp>
 
+void wells_static_check(const Wells* wells) {
+    BOOST_CHECK_EQUAL(2, wells->number_of_wells);
+    BOOST_CHECK_EQUAL(3, wells->number_of_phases);
 
-        void wells_static_check(const Wells* wells) {
-            BOOST_CHECK_EQUAL(2 , wells->number_of_wells);
-            BOOST_CHECK_EQUAL(3 , wells->number_of_phases);
+    BOOST_CHECK_EQUAL("INJ1", wells->name[0]);
+    BOOST_CHECK_EQUAL("PROD1", wells->name[1]);
 
-            BOOST_CHECK_EQUAL("INJ1" , wells->name[0]);
-            BOOST_CHECK_EQUAL("PROD1" , wells->name[1]);
+    /* The mapping from well number into the wells->WI and wells->well_cells arrays. */
+    BOOST_CHECK_EQUAL(0, wells->well_connpos[0]);
+    BOOST_CHECK_EQUAL(1, wells->well_connpos[1]);
+    BOOST_CHECK_EQUAL(2, wells->well_connpos[2]);
 
-            /* The mapping from well number into the wells->WI and wells->well_cells arrays. */
-            BOOST_CHECK_EQUAL(0 , wells->well_connpos[0]);
-            BOOST_CHECK_EQUAL(1 , wells->well_connpos[1]);
-            BOOST_CHECK_EQUAL(2 , wells->well_connpos[2]);
+    /* Connection factor */
+    BOOST_CHECK_CLOSE(1.2279166666666664e-12, wells->WI[0], 0.001);
+    BOOST_CHECK_CLOSE(1.2279166666666664e-12, wells->WI[1], 0.001);
 
-            /* Connection factor */
-            BOOST_CHECK_CLOSE(1.2279166666666664e-12 , wells->WI[0] , 0.001);
-            BOOST_CHECK_CLOSE(1.2279166666666664e-12 , wells->WI[1] , 0.001);
+    /* Completed cells */
+    BOOST_CHECK_EQUAL(0, wells->well_cells[0]);
+    BOOST_CHECK_EQUAL(9 + 2 * 10 + 2 * 10 * 10, wells->well_cells[1]);
+}
 
-            /* Completed cells */
-            BOOST_CHECK_EQUAL(0 , wells->well_cells[0]);
-            BOOST_CHECK_EQUAL(9 + 2*10 + 2*10*10 , wells->well_cells[1]);
+/*
+ The number of controls is determined by looking at which elements
+ have been given explicit - non-default - values in the WCONxxxx
+ keyword. Is that at all interesting?
+ */
+
+void check_controls_epoch0(struct WellControls ** ctrls) {
+    // The injector
+    {
+        const struct WellControls * ctrls0 = ctrls[0];
+        BOOST_CHECK_EQUAL(3, well_controls_get_num(ctrls0));   // The number of controls for the injector == 3??
+
+        BOOST_CHECK_EQUAL(SURFACE_RATE, well_controls_iget_type(ctrls0, 0));
+        BOOST_CHECK_EQUAL(RESERVOIR_RATE, well_controls_iget_type(ctrls0, 1));
+        BOOST_CHECK_EQUAL(BHP, well_controls_iget_type(ctrls0, 2));
+
+        // The different targets
+        BOOST_CHECK_EQUAL(100.0 / 86400, well_controls_iget_target(ctrls0, 0));
+        BOOST_CHECK_EQUAL(200.0 / 86400, well_controls_iget_target(ctrls0, 1));
+        BOOST_CHECK_EQUAL(400 * 100000, well_controls_iget_target(ctrls0, 2));
+
+        // Which control is active
+        BOOST_CHECK_EQUAL(0, well_controls_get_current(ctrls0));
+
+        // The phase distribution in the active target
+        {
+            const double * distr = well_controls_iget_distr(ctrls0, 0);
+            BOOST_CHECK_EQUAL(0, distr[0]);  // Water
+            BOOST_CHECK_EQUAL(0, distr[1]);  // Oil
+            BOOST_CHECK_EQUAL(1, distr[2]);  // Gas
         }
+    }
 
+    // The producer
+    {
+        const struct WellControls * ctrls1 = ctrls[1];
+        BOOST_CHECK_EQUAL(2, well_controls_get_num(ctrls1));   // The number of controls for the producer == 2??
+        BOOST_CHECK_EQUAL(SURFACE_RATE, well_controls_iget_type(ctrls1, 0));
+        BOOST_CHECK_EQUAL(BHP, well_controls_iget_type(ctrls1, 1));
 
-        /*
-           The number of controls is determined by looking at which elements
-           have been given explicit - non-default - values in the WCONxxxx
-           keyword. Is that at all interesting?
-        */
+        // The different targets
+        BOOST_CHECK_EQUAL(-20000.0 / 86400, well_controls_iget_target(ctrls1, 0));
+        BOOST_CHECK_EQUAL(1000 * 100000, well_controls_iget_target(ctrls1, 1));
 
+        // Which control is active
+        BOOST_CHECK_EQUAL(0, well_controls_get_current(ctrls1));
 
-        void check_controls_epoch0( struct WellControls ** ctrls) {
-            // The injector
-            {
-                const struct WellControls * ctrls0 = ctrls[0];
-                BOOST_CHECK_EQUAL( 3 , well_controls_get_num(ctrls0));   // The number of controls for the injector == 3??
-
-                BOOST_CHECK_EQUAL( SURFACE_RATE   , well_controls_iget_type(ctrls0 , 0) );
-                BOOST_CHECK_EQUAL( RESERVOIR_RATE , well_controls_iget_type(ctrls0 , 1) );
-                BOOST_CHECK_EQUAL( BHP            , well_controls_iget_type(ctrls0 , 2) );
-
-                // The different targets
-                BOOST_CHECK_EQUAL( 100.0 / 86400 , well_controls_iget_target(ctrls0,0));
-                BOOST_CHECK_EQUAL( 200.0 / 86400 , well_controls_iget_target(ctrls0,1));
-                BOOST_CHECK_EQUAL( 400 * 100000  , well_controls_iget_target(ctrls0,2));
-
-                // Which control is active
-                BOOST_CHECK_EQUAL( 0 , well_controls_get_current(ctrls0) );
-
-                // The phase distribution in the active target
-                {
-                     const double * distr = well_controls_iget_distr( ctrls0 , 0 );
-                     BOOST_CHECK_EQUAL( 0 , distr[0] );  // Water
-                     BOOST_CHECK_EQUAL( 0 , distr[1] );  // Oil
-                     BOOST_CHECK_EQUAL( 1 , distr[2] );  // Gas
-                }
-            }
-
-            // The producer
-            {
-                const struct WellControls * ctrls1 = ctrls[1];
-                BOOST_CHECK_EQUAL( 2 , well_controls_get_num( ctrls1 ));   // The number of controls for the producer == 2??
-                BOOST_CHECK_EQUAL( SURFACE_RATE   , well_controls_iget_type(ctrls1 , 0) );
-                BOOST_CHECK_EQUAL( BHP            , well_controls_iget_type(ctrls1 , 1) );
-
-                // The different targets
-                BOOST_CHECK_EQUAL( -20000.0 / 86400 , well_controls_iget_target(ctrls1,0));
-                BOOST_CHECK_EQUAL(  1000 * 100000   , well_controls_iget_target(ctrls1,1));
-
-                // Which control is active
-                BOOST_CHECK_EQUAL( 0 , well_controls_get_current(ctrls1));
-
-                // The phase distribution in the active target
-               {
-                    const double * distr = well_controls_iget_distr( ctrls1 , 0 );
-                    BOOST_CHECK_EQUAL( 0 , distr[0] );  // Water
-                    BOOST_CHECK_EQUAL( 1 , distr[1] );  // Oil
-                    BOOST_CHECK_EQUAL( 0 , distr[2] );  // Gas
-                }
-            }
+        // The phase distribution in the active target
+        {
+            const double * distr = well_controls_iget_distr(ctrls1, 0);
+            BOOST_CHECK_EQUAL(0, distr[0]);  // Water
+            BOOST_CHECK_EQUAL(1, distr[1]);  // Oil
+            BOOST_CHECK_EQUAL(0, distr[2]);  // Gas
         }
+    }
+}
 
+void check_controls_epoch1(struct WellControls ** ctrls) {
+    // The injector
+    {
+        const struct WellControls * ctrls0 = ctrls[0];
+        BOOST_CHECK_EQUAL(3, well_controls_get_num(ctrls0));   // The number of controls for the injector == 3??
 
+        BOOST_CHECK_EQUAL(SURFACE_RATE, well_controls_iget_type(ctrls0, 0));
+        BOOST_CHECK_EQUAL(RESERVOIR_RATE, well_controls_iget_type(ctrls0, 1));
+        BOOST_CHECK_EQUAL(BHP, well_controls_iget_type(ctrls0, 2));
 
+        // The different targets
+        BOOST_CHECK_CLOSE(10.0 / 86400, well_controls_iget_target(ctrls0, 0), 0.001);
+        BOOST_CHECK_CLOSE(20.0 / 86400, well_controls_iget_target(ctrls0, 1), 0.001);
+        BOOST_CHECK_CLOSE(40 * 100000, well_controls_iget_target(ctrls0, 2), 0.001);
 
-        void check_controls_epoch1( struct WellControls ** ctrls) {
-            // The injector
-            {
-                const struct WellControls * ctrls0 = ctrls[0];
-                BOOST_CHECK_EQUAL( 3 , well_controls_get_num(ctrls0));   // The number of controls for the injector == 3??
+        // Which control is active
+        BOOST_CHECK_EQUAL(1, well_controls_get_current(ctrls0));
 
-                BOOST_CHECK_EQUAL( SURFACE_RATE   , well_controls_iget_type(ctrls0 , 0 ));
-                BOOST_CHECK_EQUAL( RESERVOIR_RATE , well_controls_iget_type(ctrls0 , 1 ));
-                BOOST_CHECK_EQUAL( BHP            , well_controls_iget_type(ctrls0 , 2 ));
-
-                // The different targets
-                BOOST_CHECK_CLOSE( 10.0 / 86400 , well_controls_iget_target(ctrls0 , 0) , 0.001);
-                BOOST_CHECK_CLOSE( 20.0 / 86400 , well_controls_iget_target(ctrls0 , 1) , 0.001);
-                BOOST_CHECK_CLOSE( 40 * 100000  , well_controls_iget_target(ctrls0 , 2) , 0.001);
-
-                // Which control is active
-                BOOST_CHECK_EQUAL( 1 , well_controls_get_current(ctrls0));
-
-                {
-                    const double * distr = well_controls_iget_distr( ctrls0 , 1 );
-                    BOOST_CHECK_EQUAL( 1 , distr[0] );  // Water
-                    BOOST_CHECK_EQUAL( 0 , distr[1] );  // Oil
-                    BOOST_CHECK_EQUAL( 0 , distr[2] );  // Gas
-                }
-            }
-
-            // The producer
-            {
-                const struct WellControls * ctrls1 = ctrls[1];
-                BOOST_CHECK_EQUAL( 3 , well_controls_get_num(ctrls1));   // The number of controls for the producer - now 3.
-                BOOST_CHECK_EQUAL( SURFACE_RATE   , well_controls_iget_type(ctrls1 , 0) );
-                BOOST_CHECK_EQUAL( RESERVOIR_RATE , well_controls_iget_type(ctrls1 , 1) );
-                BOOST_CHECK_EQUAL( BHP            , well_controls_iget_type(ctrls1 , 2) );
-
-                // The different targets
-                BOOST_CHECK_CLOSE( -999.0 / 86400 , well_controls_iget_target(ctrls1 , 0), 0.001);
-                BOOST_CHECK_CLOSE( -123.0 / 86400 , well_controls_iget_target(ctrls1 , 1), 0.001);
-                BOOST_CHECK_CLOSE(  100 * 100000  , well_controls_iget_target(ctrls1 , 2), 0.001);
-
-                // Which control is active
-                BOOST_CHECK_EQUAL( 1 , well_controls_get_current(ctrls1) );
-
-                {
-                    const double * distr = well_controls_iget_distr( ctrls1 , 1 );
-                    BOOST_CHECK_EQUAL( 1 , distr[0] );  // Water
-                    BOOST_CHECK_EQUAL( 1 , distr[1] );  // Oil
-                    BOOST_CHECK_EQUAL( 1 , distr[2] );  // Gas
-                }
-            }
+        {
+            const double * distr = well_controls_iget_distr(ctrls0, 1);
+            BOOST_CHECK_EQUAL(1, distr[0]);  // Water
+            BOOST_CHECK_EQUAL(0, distr[1]);  // Oil
+            BOOST_CHECK_EQUAL(0, distr[2]);  // Gas
         }
+    }
 
+    // The producer
+    {
+        const struct WellControls * ctrls1 = ctrls[1];
+        BOOST_CHECK_EQUAL(3, well_controls_get_num(ctrls1));   // The number of controls for the producer - now 3.
+        BOOST_CHECK_EQUAL(SURFACE_RATE, well_controls_iget_type(ctrls1, 0));
+        BOOST_CHECK_EQUAL(RESERVOIR_RATE, well_controls_iget_type(ctrls1, 1));
+        BOOST_CHECK_EQUAL(BHP, well_controls_iget_type(ctrls1, 2));
 
-        void check_controls_epoch3( struct WellControls ** ctrls) {
-            // The new producer
-            const struct WellControls * ctrls1 = ctrls[1];
-            // Note: controls include default (1 atm) BHP control.
-            BOOST_CHECK_EQUAL( 6 , well_controls_get_num(ctrls1));
+        // The different targets
+        BOOST_CHECK_CLOSE(-999.0 / 86400, well_controls_iget_target(ctrls1, 0), 0.001);
+        BOOST_CHECK_CLOSE(-123.0 / 86400, well_controls_iget_target(ctrls1, 1), 0.001);
+        BOOST_CHECK_CLOSE(100 * 100000, well_controls_iget_target(ctrls1, 2), 0.001);
+
+        // Which control is active
+        BOOST_CHECK_EQUAL(1, well_controls_get_current(ctrls1));
+
+        {
+            const double * distr = well_controls_iget_distr(ctrls1, 1);
+            BOOST_CHECK_EQUAL(1, distr[0]);  // Water
+            BOOST_CHECK_EQUAL(1, distr[1]);  // Oil
+            BOOST_CHECK_EQUAL(1, distr[2]);  // Gas
         }
+    }
+}
 
+void check_controls_epoch3(struct WellControls ** ctrls) {
+    // The new producer
+    const struct WellControls * ctrls1 = ctrls[1];
+    // Note: controls include default (1 atm) BHP control.
+    BOOST_CHECK_EQUAL(6, well_controls_get_num(ctrls1));
+}
 
+BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
 
+    const std::string filename = "wells_manager_data.data";
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::ParseContext parseContext;
+    Opm::DeckConstPtr deck(parser->parseFile(filename, parseContext));
 
-        BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
+    Opm::GridManager gridManager(eclipseState->getInputGrid());
 
-            const std::string filename = "wells_manager_data.data";
-            Opm::ParserPtr parser(new Opm::Parser());
-            Opm::ParseContext parseContext;
-            Opm::DeckConstPtr deck(parser->parseFile(filename , parseContext));
+    {
+        Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
+        wells_static_check(wellsManager.c_wells());
+        check_controls_epoch0(wellsManager.c_wells()->ctrls);
+    }
 
-            Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
-            Opm::GridManager gridManager(deck);
+    {
+        Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
+        wells_static_check(wellsManager.c_wells());
+        check_controls_epoch1(wellsManager.c_wells()->ctrls);
+    }
 
-            {
-                Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
-                wells_static_check( wellsManager.c_wells() );
-                check_controls_epoch0( wellsManager.c_wells()->ctrls );
-            }
+    {
+        Opm::WellsManager wellsManager(eclipseState, 3, *gridManager.c_grid(), NULL);
+        const Wells* wells = wellsManager.c_wells();
 
-            {
-                Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
-                wells_static_check( wellsManager.c_wells() );
-                check_controls_epoch1( wellsManager.c_wells()->ctrls );
-            }
+        // There is 3 wells in total in the deck at the 3rd schedule step.
+        // PROD1 is shut and should therefore not be counted.
+        // The new well is therefore the secound well.
+        BOOST_CHECK_EQUAL(2, wells->number_of_wells);
+        BOOST_CHECK_EQUAL(wells->name[0], "INJ1");
+        BOOST_CHECK_EQUAL(wells->name[1], "NEW");
 
+        check_controls_epoch3(wellsManager.c_wells()->ctrls);
+    }
 
-            {
-                Opm::WellsManager wellsManager(eclipseState, 3, *gridManager.c_grid(), NULL);
-                const Wells* wells = wellsManager.c_wells();
+}
 
-                // There is 3 wells in total in the deck at the 3rd schedule step.
-                // PROD1 is shut and should therefore not be counted.
-                // The new well is therefore the secound well.
-                BOOST_CHECK_EQUAL(2 , wells->number_of_wells);
-                BOOST_CHECK_EQUAL( wells->name[0] , "INJ1");
-                BOOST_CHECK_EQUAL( wells->name[1] , "NEW");
+BOOST_AUTO_TEST_CASE(WellsEqual) {
+    const std::string filename = "wells_manager_data.data";
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::ParseContext parseContext;
+    Opm::DeckConstPtr deck(parser->parseFile(filename, parseContext));
 
-                check_controls_epoch3( wellsManager.c_wells()->ctrls );
-            }
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
+    Opm::GridManager gridManager(eclipseState->getInputGrid());
 
-        }
+    Opm::WellsManager wellsManager0(eclipseState, 0, *gridManager.c_grid(), NULL);
+    Opm::WellsManager wellsManager1(eclipseState, 1, *gridManager.c_grid(), NULL);
 
+    BOOST_CHECK(wells_equal( wellsManager0.c_wells() , wellsManager0.c_wells(),false));
+    BOOST_CHECK(!wells_equal( wellsManager0.c_wells() , wellsManager1.c_wells(),false));
+}
 
+BOOST_AUTO_TEST_CASE(ControlsEqual) {
+    const std::string filename = "wells_manager_data.data";
+    Opm::ParseContext parseContext;
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::DeckConstPtr deck(parser->parseFile(filename, parseContext));
 
-        BOOST_AUTO_TEST_CASE(WellsEqual) {
-            const std::string filename = "wells_manager_data.data";
-            Opm::ParserPtr parser(new Opm::Parser());
-            Opm::ParseContext parseContext;
-            Opm::DeckConstPtr deck(parser->parseFile(filename , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
+    Opm::GridManager gridManager(eclipseState->getInputGrid());
 
-            Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
-            Opm::GridManager gridManager(deck);
+    Opm::WellsManager wellsManager0(eclipseState, 0, *gridManager.c_grid(), NULL);
+    Opm::WellsManager wellsManager1(eclipseState, 1, *gridManager.c_grid(), NULL);
 
-            Opm::WellsManager wellsManager0(eclipseState , 0 , *gridManager.c_grid(), NULL);
-            Opm::WellsManager wellsManager1(eclipseState , 1 , *gridManager.c_grid(), NULL);
+    BOOST_CHECK(well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0] , false));
+    BOOST_CHECK(well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1] , false));
+    BOOST_CHECK(well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager1.c_wells()->ctrls[0] , false));
+    BOOST_CHECK(well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager1.c_wells()->ctrls[1] , false));
 
-            BOOST_CHECK(  wells_equal( wellsManager0.c_wells() , wellsManager0.c_wells(),false) );
-            BOOST_CHECK( !wells_equal( wellsManager0.c_wells() , wellsManager1.c_wells(),false) );
-        }
+    BOOST_CHECK(!well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[1] , false));
+    BOOST_CHECK(!well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[0] , false));
+    BOOST_CHECK(!well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0] , false));
+    BOOST_CHECK(!well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1] , false));
+}
 
+BOOST_AUTO_TEST_CASE(WellShutOK) {
+    const std::string filename = "wells_manager_data.data";
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::ParseContext parseContext;
+    Opm::DeckConstPtr deck(parser->parseFile(filename, parseContext));
 
-        BOOST_AUTO_TEST_CASE(ControlsEqual) {
-            const std::string filename = "wells_manager_data.data";
-            Opm::ParseContext parseContext;
-            Opm::ParserPtr parser(new Opm::Parser());
-            Opm::DeckConstPtr deck(parser->parseFile(filename , parseContext));
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
+    Opm::GridManager gridManager(eclipseState->getInputGrid());
 
-            Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
-            Opm::GridManager gridManager(deck);
+    Opm::WellsManager wellsManager2(eclipseState, 2, *gridManager.c_grid(), NULL);
 
-            Opm::WellsManager wellsManager0(eclipseState , 0 , *gridManager.c_grid(), NULL);
-            Opm::WellsManager wellsManager1(eclipseState , 1 , *gridManager.c_grid(), NULL);
+    // Shut wells are not added to the deck. i.e number of wells should be 2-1
+    BOOST_CHECK(wellsManager2.c_wells()->number_of_wells == 1);
 
-            BOOST_CHECK(  well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0] , false));
-            BOOST_CHECK(  well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1] , false));
-            BOOST_CHECK(  well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager1.c_wells()->ctrls[0] , false));
-            BOOST_CHECK(  well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager1.c_wells()->ctrls[1] , false));
+    //BOOST_CHECK_NO_THROW( Opm::WellsManager wellsManager2(eclipseState , 2 , *gridManager.c_grid(), NULL));
+}
 
-            BOOST_CHECK(  !well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[1] , false));
-            BOOST_CHECK(  !well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[0] , false));
-            BOOST_CHECK(  !well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0] , false));
-            BOOST_CHECK(  !well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1] , false));
-        }
+BOOST_AUTO_TEST_CASE(WellSTOPOK) {
+    const std::string filename = "wells_manager_data_wellSTOP.data";
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::ParseContext parseContext;
+    Opm::DeckConstPtr deck(parser->parseFile(filename, parseContext));
 
+    Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
+    Opm::GridManager gridManager(eclipseState->getInputGrid());
 
-
-        BOOST_AUTO_TEST_CASE(WellShutOK) {
-            const std::string filename = "wells_manager_data.data";
-            Opm::ParserPtr parser(new Opm::Parser());
-            Opm::ParseContext parseContext;
-            Opm::DeckConstPtr deck(parser->parseFile(filename , parseContext));
-
-            Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
-            Opm::GridManager gridManager(deck);
-
-            Opm::WellsManager wellsManager2(eclipseState , 2 , *gridManager.c_grid(), NULL);
-
-            // Shut wells are not added to the deck. i.e number of wells should be 2-1
-            BOOST_CHECK( wellsManager2.c_wells()->number_of_wells == 1);
-
-            //BOOST_CHECK_NO_THROW( Opm::WellsManager wellsManager2(eclipseState , 2 , *gridManager.c_grid(), NULL));
-        }
-
-
-
-        BOOST_AUTO_TEST_CASE(WellSTOPOK) {
-            const std::string filename = "wells_manager_data_wellSTOP.data";
-            Opm::ParserPtr parser(new Opm::Parser());
-            Opm::ParseContext parseContext;
-            Opm::DeckConstPtr deck(parser->parseFile(filename , parseContext));
-
-            Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck , parseContext));
-    Opm::GridManager gridManager(deck);
-
-    Opm::WellsManager wellsManager(eclipseState , 0 , *gridManager.c_grid(), NULL);
+    Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
 
     const Wells* wells = wellsManager.c_wells();
     const struct WellControls* ctrls0 = wells->ctrls[0];