diff --git a/ApplicationCode/ProjectDataModel/RimEclipseInputCase.cpp b/ApplicationCode/ProjectDataModel/RimEclipseInputCase.cpp
index 6af6cc2056..4bbefbf730 100644
--- a/ApplicationCode/ProjectDataModel/RimEclipseInputCase.cpp
+++ b/ApplicationCode/ProjectDataModel/RimEclipseInputCase.cpp
@@ -249,9 +249,9 @@ void RimEclipseInputCase::loadAndSyncronizeInputProperties()
     // Then read the properties from all the files referenced by the InputReservoir
 
     std::vector<QString> filenames;
-    for (const RimEclipseInputProperty* inputProperty : m_inputPropertyCollection()->inputProperties())
+    for (const QString& fileName : additionalFiles())
     {
-        filenames.push_back(inputProperty->fileName);
+        filenames.push_back(fileName);
     }
     filenames.push_back(m_gridFileName);
 
diff --git a/ApplicationCode/ReservoirDataModel/RigNumberOfFloodedPoreVolumesCalculator.cpp b/ApplicationCode/ReservoirDataModel/RigNumberOfFloodedPoreVolumesCalculator.cpp
index 065270c68c..09f43e0821 100644
--- a/ApplicationCode/ReservoirDataModel/RigNumberOfFloodedPoreVolumesCalculator.cpp
+++ b/ApplicationCode/ReservoirDataModel/RigNumberOfFloodedPoreVolumesCalculator.cpp
@@ -236,7 +236,7 @@ void RigNumberOfFloodedPoreVolumesCalculator::calculate(RigMainGrid* mainGrid,
 
         const std::vector<double>* flowrateNNC = flowrateNNCatAllTimeSteps[timeStep-1];
 
-        if (flowrateNNC->size() > 0)
+        if (flowrateNNC && flowrateNNC->size() > 0)
         {
             distributeNNCflow(connections,
                               caseToApply,
diff --git a/ApplicationCode/UserInterface/RiuResultQwtPlot.cpp b/ApplicationCode/UserInterface/RiuResultQwtPlot.cpp
index 0cef8c26c6..f6bec832ac 100644
--- a/ApplicationCode/UserInterface/RiuResultQwtPlot.cpp
+++ b/ApplicationCode/UserInterface/RiuResultQwtPlot.cpp
@@ -66,6 +66,11 @@ RiuResultQwtPlot::~RiuResultQwtPlot()
 //--------------------------------------------------------------------------------------------------
 void RiuResultQwtPlot::addCurve(const QString& curveName, const cvf::Color3f& curveColor, const std::vector<QDateTime>& dateTimes, const std::vector<double>& timeHistoryValues)
 {
+    if (dateTimes.empty() || timeHistoryValues.empty())
+    {
+        return;
+    }
+
     RiuLineSegmentQwtPlotCurve* plotCurve = new RiuLineSegmentQwtPlotCurve("Curve 1");
 
     plotCurve->setSamplesFromDatesAndYValues(dateTimes, timeHistoryValues, false);
diff --git a/ResInsightVersion.cmake b/ResInsightVersion.cmake
index efce40c205..6ec04eeef5 100644
--- a/ResInsightVersion.cmake
+++ b/ResInsightVersion.cmake
@@ -10,7 +10,7 @@ set(RESINSIGHT_VERSION_TEXT "-dev")
 # Must be unique and increasing within one combination of major/minor/patch version 
 # The uniqueness of this text is independent of RESINSIGHT_VERSION_TEXT 
 # Format of text must be ".xx"
-set(RESINSIGHT_DEV_VERSION ".02")
+set(RESINSIGHT_DEV_VERSION ".103")
 
 # https://github.com/CRAVA/crava/tree/master/libs/nrlib
 set(NRLIB_GITHUB_SHA "ba35d4359882f1c6f5e9dc30eb95fe52af50fd6f") 
@@ -22,7 +22,7 @@ set(ERT_GITHUB_SHA "2e36798b43daf18c112b91aa3febbf2fccd4a95f")
 set(OPM_FLOWDIAGNOSTICS_SHA "7e2be931d430796ed42efcfb5c6b67a8d5962f7f")
 
 # https://github.com/OPM/opm-flowdiagnostics-applications
-set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "44f7e47ecdc87ba566ab4146629de49039a73b2e")
+set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "5bcd6d99259a63f5cd820db541b45c4f07aec808")
 
 # https://github.com/OPM/opm-parser/blob/master/opm/parser/eclipse/Units/Units.hpp
 # This file was moved from opm-core to opm-parser october 2016
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists.txt b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists.txt
index 3f086ceb51..565b1aafce 100644
--- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists.txt
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists.txt
@@ -20,7 +20,7 @@ cmake_minimum_required (VERSION 2.8)
 option(SIBLING_SEARCH "Search for other modules in sibling directories?" ON)
 
 # Mandatory call to project
-project(opm-flowdiagnostics-applications CXX)
+project(opm-flowdiagnostics-applications C CXX)
 
 if(SIBLING_SEARCH AND NOT opm-common_DIR)
   # guess the sibling dir
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/dune.module b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/dune.module
index 1285983fed..b55f213284 100644
--- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/dune.module
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/dune.module
@@ -5,8 +5,8 @@
 
 Module: opm-flowdiagnostics-applications
 Description: flow diagnostics applications and examples
-Version: 2016.10-pre
-Label: 2016.10-pre
+Version: 2018.04-pre
+Label: 2018.04-pre
 Maintainer: bard.skaflestad@sintef.no
 MaintainerName: Bård Skaflestad
 Url: http://opm-project.org
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm-flowdiagnostics-applications-prereqs.cmake b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm-flowdiagnostics-applications-prereqs.cmake
new file mode 100644
index 0000000000..5461c566c8
--- /dev/null
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm-flowdiagnostics-applications-prereqs.cmake
@@ -0,0 +1,26 @@
+# defines that must be present in config.h for our headers
+set (opm-flowdiagnostics-applications_CONFIG_VAR
+  )
+
+# Build prerequisites
+set (opm-flowdiagnostics-applications_DEPS
+  # This module requires C++11 support, including std::regex
+  "CXX11Features REQUIRED"
+  # We need Boost.Filesystem for advanced file handling
+  #   filesystem::path
+  #   filesystem::directory_iterator
+  #   filesystem::last_write_time()
+  "Boost 1.44.0
+    COMPONENTS filesystem system unit_test_framework REQUIRED"
+  # We need LibECL to handle ECL result sets.
+  "ecl REQUIRED"
+  # Prerequisite OPM modules
+  #   common -> Parameter System
+  #   fdiag  -> Solver
+  #   parser -> Unit Conversions
+  "opm-common REQUIRED"
+  "opm-flowdiagnostics REQUIRED"
+  "opm-parser REQUIRED"
+  )
+
+find_package_deps(opm-flowdiagnostics-applications)
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.cpp
index 5b03b0690f..0c9a99bc65 100644
--- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.cpp
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.cpp
@@ -81,6 +81,16 @@ namespace {
         return tep;
     }
 
+    double defaultedScaledSaturation(const double s, const double dflt)
+    {
+        // Use input scaled saturation ('s') if not defaulted (|s| < 1e20),
+        // default scaled saturation otherwise.  The sentinel value 1e20 is
+        // the common value used to represent unset/defaulted values in ECL
+        // result sets.
+
+        return (std::abs(s) < 1.0e+20) ? s : dflt;
+    }
+
     bool validSaturation(const double s)
     {
         return (! (s < 0.0)) && (! (s > 1.0));
@@ -135,6 +145,24 @@ private:
 
     void handleInvalidEndpoint(const SaturationAssoc& sp,
                                std::vector<double>&   effsat) const;
+
+    double sMin(const std::vector<int>::size_type cell,
+                const TableEndPoints&             tep) const
+    {
+        // Use model's scaled connate saturation if defined, otherwise fall
+        // back to table's unscaled connate saturation.
+
+        return defaultedScaledSaturation(this->smin_[cell], tep.low);
+    }
+
+    double sMax(const std::vector<int>::size_type cell,
+                const TableEndPoints&             tep) const
+    {
+        // Use model's scaled maximum saturation if defined, otherwise fall
+        // back to table's unscaled maximum saturation.
+
+        return defaultedScaledSaturation(this->smax_[cell], tep.high);
+    }
 };
 
 std::vector<double>
@@ -150,8 +178,8 @@ Impl::eval(const TableEndPoints&   tep,
     for (const auto& eval_pt : sp) {
         const auto cell = eval_pt.cell;
 
-        const auto sLO = this->smin_[cell];
-        const auto sHI = this->smax_[cell];
+        const auto sLO = this->sMin(cell, tep);
+        const auto sHI = this->sMax(cell, tep);
 
         if (! validSaturations({ sLO, sHI })) {
             this->handleInvalidEndpoint(eval_pt, effsat);
@@ -195,8 +223,8 @@ Impl::reverse(const TableEndPoints&   tep,
     for (const auto& eval_pt : sp) {
         const auto cell = eval_pt.cell;
 
-        const auto sLO = this->smin_[cell];
-        const auto sHI = this->smax_[cell];
+        const auto sLO = this->sMin(cell, tep);
+        const auto sHI = this->sMax(cell, tep);
 
         if (! validSaturations({ sLO, sHI })) {
             this->handleInvalidEndpoint(eval_pt, unscaledsat);
@@ -293,6 +321,33 @@ private:
 
     void handleInvalidEndpoint(const SaturationAssoc& sp,
                                std::vector<double>&   effsat) const;
+
+    double sMin(const std::vector<int>::size_type cell,
+                const TableEndPoints&             tep) const
+    {
+        // Use model's scaled connate saturation if defined, otherwise fall
+        // back to table's unscaled connate saturation.
+
+        return defaultedScaledSaturation(this->smin_[cell], tep.low);
+    }
+
+    double sDisp(const std::vector<int>::size_type cell,
+                 const TableEndPoints&             tep) const
+    {
+        // Use model's scaled displacing saturation if defined, otherwise
+        // fall back to table's unscaled displacing saturation.
+
+        return defaultedScaledSaturation(this->sdisp_[cell], tep.disp);
+    }
+
+    double sMax(const std::vector<int>::size_type cell,
+                const TableEndPoints&             tep) const
+    {
+        // Use model's scaled maximum saturation if defined, otherwise fall
+        // back to table's unscaled maximum saturation.
+
+        return defaultedScaledSaturation(this->smax_[cell], tep.high);
+    }
 };
 
 std::vector<double>
@@ -306,9 +361,9 @@ Impl::eval(const TableEndPoints&   tep,
     for (const auto& eval_pt : sp) {
         const auto cell = eval_pt.cell;
 
-        const auto sLO = this->smin_ [cell];
-        const auto sR  = this->sdisp_[cell];
-        const auto sHI = this->smax_ [cell];
+        const auto sLO = this->sMin (cell, tep);
+        const auto sR  = this->sDisp(cell, tep);
+        const auto sHI = this->sMax (cell, tep);
 
         if (! validSaturations({ sLO, sR, sHI })) {
             this->handleInvalidEndpoint(eval_pt, effsat);
@@ -355,9 +410,9 @@ Impl::reverse(const TableEndPoints&   tep,
     for (const auto& eval_pt : sp) {
         const auto cell = eval_pt.cell;
 
-        const auto sLO = this->smin_ [cell];
-        const auto sR  = this->sdisp_[cell];
-        const auto sHI = this->smax_ [cell];
+        const auto sLO = this->sMin (cell, tep);
+        const auto sR  = this->sDisp(cell, tep);
+        const auto sHI = this->sMax (cell, tep);
 
         if (! validSaturations({ sLO, sR, sHI })) {
             this->handleInvalidEndpoint(eval_pt, unscaledsat);
@@ -678,7 +733,13 @@ Create::TwoPoint::Pc::GO(const ::Opm::ECLGraph&        G,
                          const ::Opm::ECLInitFileData& init,
                          const InvBeh                  handle_invalid)
 {
-    auto sgl = G.rawLinearisedCellData<double>(init, "SGL");
+    // Try dedicated scaled Sg_conn for Pc first
+    auto sgl = G.rawLinearisedCellData<double>(init, "SGLPC");
+    if (sgl.empty()) {
+        // Fall back to general scaled Sg_conn if not available.
+        sgl = G.rawLinearisedCellData<double>(init, "SGL");
+    }
+
     auto sgu = G.rawLinearisedCellData<double>(init, "SGU");
 
     if ((sgl.size() != sgu.size()) ||
@@ -702,7 +763,13 @@ Create::TwoPoint::Pc::OW(const ::Opm::ECLGraph&        G,
                          const ::Opm::ECLInitFileData& init,
                          const InvBeh                  handle_invalid)
 {
-    auto swl = G.rawLinearisedCellData<double>(init, "SWL");
+    // Try dedicated scaled Sw_conn for Pc first
+    auto swl = G.rawLinearisedCellData<double>(init, "SWLPC");
+    if (swl.empty()) {
+        // Fall back to general scaled Sw_conn if not available.
+        swl = G.rawLinearisedCellData<double>(init, "SWL");
+    }
+
     auto swu = G.rawLinearisedCellData<double>(init, "SWU");
 
     if ((swl.size() != swu.size()) ||
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp
index 473635ce15..5550d376c3 100644
--- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp
@@ -20,7 +20,7 @@
 #include <opm/utility/ECLFluxCalc.hpp>
 #include <opm/utility/ECLPvtCommon.hpp>
 #include <opm/utility/ECLUnitHandling.hpp>
-
+#include <opm/utility/ECLSaturationFunc.hpp>
 #include <opm/parser/eclipse/Units/Units.hpp>
 
 #include <algorithm>
@@ -221,6 +221,21 @@ namespace Opm
     }
 
 
+    std::vector<double>
+    ECLFluxCalc::massflux(const ECLRestartData& rstrt,
+                          const ECLPhaseIndex   phase) const
+    {
+        // Obtain dynamic data.
+        const auto dyn_data = this->phaseProperties(rstrt, phase);
+
+        // Compute fluxes per connection.
+        const int num_conn = transmissibility_.size();
+        std::vector<double> fluxvec(num_conn);
+        for (int conn = 0; conn < num_conn; ++conn) {
+            fluxvec[conn] = singleMassFlux(conn, dyn_data);
+        }
+        return fluxvec;
+    }
 
 
 
@@ -251,7 +266,35 @@ namespace Opm
         return mob * T * dh;
     }
 
+    double ECLFluxCalc::singleMassFlux(const int connection,
+                                       const DynamicData& dyn_data) const
+    {
+        const int c1 = neighbours_[2*connection];
+        const int c2 = neighbours_[2*connection + 1];
 
+        // Phase pressure in connecting cells.
+        const auto p1 = dyn_data.pressure[c1];
+        const auto p2 = dyn_data.pressure[c2];
+
+        // Phase density at interface: Arith. avg. of cell values.
+        const auto rho =
+            (dyn_data.density[c1] + dyn_data.density[c2]) / 2.0;
+
+        // Phase potential drop across interface.
+        const auto dh = p1 - p2 + rho*this->gravDz_[connection];
+
+        // Phase mobility at interface: Upstream weighting (phase pot).
+        const auto ucell = (dh < 0.0) ? c2 : c1;
+        const auto mob   = dyn_data.mobility[ucell];
+
+        // Background (static) transmissibility.
+        const auto T = this->transmissibility_[connection];
+
+        // Upstream weighted phase density.
+        const auto urho = dyn_data.density[ucell];
+
+        return urho * mob * T * dh;
+    }
 
 
 
@@ -278,12 +321,39 @@ namespace Opm
 
         switch (phase) {
         case ECLPhaseIndex::Aqua:
+            dyn_data.saturation = this->graph_.rawLinearisedCellData<double>(rstrt, "SWAT");
             return this->watPVT(std::move(dyn_data));
 
         case ECLPhaseIndex::Liquid:
-            return this->oilPVT(rstrt, std::move(dyn_data));
+            dyn_data.saturation = this->graph_.rawLinearisedCellData<double>(rstrt, "SOIL");
+            if (!dyn_data.saturation.empty()) {
+                return this->oilPVT(rstrt, std::move(dyn_data));
+            } else {
+                // SOIL vector not provided. Compute from SWAT and/or SGAS.
+                // may read two times
+                auto sw = this->graph_.rawLinearisedCellData<double>(rstrt, "SWAT");
+                auto sg = this->graph_.rawLinearisedCellData<double>(rstrt, "SGAS");
+                std::vector<double>& so = dyn_data.saturation;
+                so.assign(this->graph_.numCells(), 1.0);
+                auto adjust_So_for_other_phase =
+                    [&so](const std::vector<double>& s)
+                {
+                    std::transform(std::begin(so), std::end(so),
+                                   std::begin(s) ,
+                                   std::begin(so), std::minus<double>());
+                };
+                if (sg.size() == this->graph_.numCells()) {
+                    adjust_So_for_other_phase(sg);
+                }
+
+                if (sw.size() == this->graph_.numCells()) {
+                    adjust_So_for_other_phase(sw);
+                }
+                return this->oilPVT(rstrt, std::move(dyn_data));
+            }
 
         case ECLPhaseIndex::Vapour:
+            dyn_data.saturation = this->graph_.rawLinearisedCellData<double>(rstrt, "SGAS");
             return this->gasPVT(rstrt, std::move(dyn_data));
         }
 
@@ -292,7 +362,18 @@ namespace Opm
         };
     }
 
+    double ECLFluxCalc::surfaceDensity(const ECLPhaseIndex   phase) const{
+        switch (phase) {
+        case ECLPhaseIndex::Aqua:
+            return this->pvtWat_->surfaceMassDensity(0);
 
+        case ECLPhaseIndex::Liquid:
+            return this->pvtOil_->surfaceMassDensity(0);
+
+        case ECLPhaseIndex::Vapour:
+            return this->pvtGas_->surfaceMassDensity(0);
+        }
+    }
 
 
 
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp
index a229983dfe..341f91326a 100644
--- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp
@@ -33,6 +33,7 @@
 
 namespace Opm
 {
+
     class ECLRestartData;
     class ECLInitFileData;
 
@@ -73,19 +74,62 @@ namespace Opm
         flux(const ECLRestartData& rstrt,
              const ECLPhaseIndex   phase) const;
 
-    private:
+        /// Retrive phase mass flux on all connections defined by \code
+        /// graph.neighbours() \endcode.
+        ///
+        /// \param[in] rstrt ECL Restart data set from which to extract
+        ///            relevant data per cell.
+        ///
+        /// \param[in] phase Canonical phase for which to retrive flux.
+        ///
+        /// \return Mass flux values corresponding to selected phase.
+        ///         Empty if required data is missing.
+        ///         Numerical values in SI units (kg/s).
+        std::vector<double>
+        massflux(const ECLRestartData& rstrt,
+                 const ECLPhaseIndex   phase) const;
+
+        /// Return type for the phaseProperties() method,
+        /// encapsulates dynamic properties for a single
+        /// phase.
         struct DynamicData
         {
             std::vector<double> pressure;
             std::vector<double> mobility;
             std::vector<double> density;
+            std::vector<double> saturation;
         };
 
+        /// Retrive dynamical properties of a single phase on all cells.
+        ///
+        /// \param[in] rstrt ECL Restart data set from which to extract
+        ///            relevant data per cell.
+        ///
+        /// \param[in] phase Canonical phase for which to retrive properties.
+        ///
+        /// \return DynamicData struct containing cell-values for phase properties.
+        ///         Numerical values in SI units (kg/s).
+        DynamicData phaseProperties(const ECLRestartData& rstrt,
+                                    const ECLPhaseIndex   phase) const;
+
+        /// Retrive the constant surface density of a phase.
+        ///
+        /// \param[in] phase Canonical phase for which to retrive the surface density.
+        ///
+        /// \return Density of given phase at surface conditions.
+        ///         Numerical value in SI units (kg/m^3).
+        double surfaceDensity(const ECLPhaseIndex phase) const;
+
+
+    private:
+
+
         double singleFlux(const int connection,
                           const DynamicData& dyn_data) const;
 
-        DynamicData phaseProperties(const ECLRestartData& rstrt,
-                                    const ECLPhaseIndex   phase) const;
+        double singleMassFlux(const int connection,
+                              const DynamicData& dyn_data) const;
+
 
         DynamicData gasPVT(const ECLRestartData& rstrt,
                            DynamicData&&         dyn_data) const;
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp
index f56f880c35..8783275b98 100644
--- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp
@@ -195,6 +195,8 @@ namespace {
             /// \endcode.
             const std::vector<int>& neighbours() const;
 
+            const std::vector<int>& getLocalCellID() const;
+
             /// Retrive static pore-volume values on active cells only.
             ///
             /// Corresponds to the \c PORV vector in the INIT file, possibly
@@ -349,6 +351,8 @@ namespace {
                 std::size_t
                 getGlobalCell(const int i, const int j, const int k) const;
 
+                const std::vector<int>& getLocalCellID() const;
+
                 /// Retrieve active cell ID of particular global cell's
                 /// neighbour in given Cartesian direction.
                 ///
@@ -872,6 +876,13 @@ CartesianCells::getGlobalCell(const int i, const int j, const int k) const
     return this->globIdx(ijk);
 }
 
+const std::vector<int>&
+ECL::CartesianGridData::
+CartesianCells::getLocalCellID() const
+{
+    return this->active_ID_;
+}
+
 int
 ECL::CartesianGridData::
 CartesianCells::getNeighbour(const std::size_t globalCell,
@@ -1006,6 +1017,12 @@ ECL::CartesianGridData::neighbours() const
     return this->neigh_;
 }
 
+const std::vector<int>&
+ECL::CartesianGridData::getLocalCellID() const
+{
+    return this->cells_.getLocalCellID();
+}
+
 const std::vector<double>&
 ECL::CartesianGridData::activePoreVolume() const
 {
@@ -1292,6 +1309,8 @@ public:
     /// \endcode.
     std::vector<int> neighbours() const;
 
+    std::vector<int> localCellID() const;
+
     /// Retrieve static pore-volume values on active cells only.
     ///
     /// Corresponds to the \c PORV vector in the INIT file, possibly
@@ -1372,6 +1391,12 @@ public:
                        const std::string&    vector,
                        UnitConvention        unit) const;
 
+
+    const ECLGeometryGrid& getGeometryGrid() const
+    {
+        return this->geomGrid_;
+    }
+
 private:
     /// Collection of non-Cartesian neighbourship relations attributed to a
     /// particular ECL keyword set (i.e., one of NNC{1,2}, NNC{G,L}, NNCLL).
@@ -1614,8 +1639,12 @@ private:
     /// Set of active grids in result set.
     std::vector<std::string> activeGrids_;
 
+    /// Map grid names to numeric grid IDs
     std::unordered_map<std::string, int> gridID_;
 
+    /// Geometry for main grid.
+    ECLGeometryGrid geomGrid_;
+
     /// Extract explicit non-neighbouring connections from ECL output.
     ///
     /// Writes to \c neigh_ and \c nncID_.
@@ -1885,6 +1914,46 @@ Opm::ECLGraph::Impl::Impl(const boost::filesystem::path& grid,
 
     this->defineNNCs(G.get(), init);
     this->defineActivePhases(init);
+
+    // Extract geometry of main grid.
+    {
+        // Size
+        {
+            this->geomGrid_.size[0] = ecl_grid_get_nx(G.get());
+            this->geomGrid_.size[1] = ecl_grid_get_ny(G.get());
+            this->geomGrid_.size[2] = ecl_grid_get_nz(G.get());
+        }
+
+        // COORD
+        {
+            const auto ncoord =
+                    static_cast<std::size_t>(ecl_grid_get_coord_size(G.get()));
+            this->geomGrid_.coord.resize(ncoord, 0.0);
+
+            ecl_grid_init_coord_data_double(G.get(), this->geomGrid_.coord.data());
+        }
+
+        // ZCORN
+        {
+            const auto nzcorn =
+                    static_cast<std::size_t>(ecl_grid_get_zcorn_size(G.get()));
+            this->geomGrid_.zcorn.resize(nzcorn, 0.0);
+
+            ecl_grid_init_zcorn_data_double(G.get(), this->geomGrid_.zcorn.data());
+        }
+
+        // ACTNUM
+        {
+            const auto nglob =
+                    static_cast<std::size_t>(this->geomGrid_.size[0]) *
+                    static_cast<std::size_t>(this->geomGrid_.size[1]) *
+                    static_cast<std::size_t>(this->geomGrid_.size[2]);
+
+            this->geomGrid_.actnum.assign(nglob, 1);
+
+            ecl_grid_init_actnum_data(G.get(), this->geomGrid_.actnum.data());
+        }
+    }
 }
 
 int
@@ -1983,6 +2052,28 @@ Opm::ECLGraph::Impl::neighbours() const
     return N;
 }
 
+std::vector<int>
+Opm::ECLGraph::Impl::localCellID() const
+{
+    auto ret = std::vector<int>{};
+    ret.reserve(this->numCells());
+
+    auto off = this->activeOffset_.begin();
+    for (const auto& G : this->grid_) {
+        for (const auto& loc : G.getLocalCellID()) {
+            const auto locID = (loc >= 0)
+                    ? static_cast<int>(*off) + loc
+                    : -1;
+
+            ret.push_back(locID);
+        }
+
+        ++off;
+    }
+
+    return ret;
+}
+
 std::vector<double>
 Opm::ECLGraph::Impl::activePoreVolume() const
 {
@@ -2322,6 +2413,11 @@ std::size_t Opm::ECLGraph::numConnections() const
     return this->pImpl_->numConnections();
 }
 
+std::vector<int> Opm::ECLGraph::localCellID() const
+{
+    return this->pImpl_->localCellID();
+}
+
 const std::vector<Opm::ECLPhaseIndex        >&
 Opm::ECLGraph::activePhases() const
 {
@@ -2393,3 +2489,9 @@ Opm::ECLGraph::linearisedCellData(const ECLRestartData& rstrt,
 {
     return this->pImpl_->linearisedCellData(rstrt, vector, unit);
 }
+
+const Opm::ECLGraph::ECLGeometryGrid&
+Opm::ECLGraph::getGeometryGrid() const
+{
+    return this->pImpl_->getGeometryGrid();
+}
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp
index 4ad8bc7490..e8114d7127 100644
--- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp
@@ -119,6 +119,9 @@ namespace Opm {
         /// Retrieve number of connections in graph.
         std::size_t numConnections() const;
 
+        /// Retrieve active cell indices for all global cells in all grids.
+        std::vector<int> localCellID() const;
+
         /// Retrieve the simulation scenario's active phases.
         ///
         /// Mostly useful to determine the set of \c PhaseIndex values for
@@ -211,6 +214,18 @@ namespace Opm {
                            const std::string&    vector,
                            UnitConvention        unit) const;
 
+        /// Raw grid data (geometric grid) for corner-point grids.
+        struct ECLGeometryGrid {
+            std::array<int, 3> size;
+            std::vector<double> zcorn;
+            std::vector<double> coord;
+            std::vector<int> actnum;
+        };
+
+        /// Retrieve corner-point grid definition,
+        /// only for main grid in case of local grid refinement (LGR).
+        const ECLGeometryGrid& getGeometryGrid() const;
+
     private:
         /// Implementation class.
         class Impl;
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.cpp
index 77d5db5942..e84cba83b9 100644
--- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.cpp
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.cpp
@@ -593,10 +593,34 @@ fromECLOutput(const ECLInitFileData& init)
     const auto& tabdims = init.keywordData<int>("TABDIMS");
     const auto& tab     = init.keywordData<double>("TAB");
 
-    raw.numPrimary = tabdims[ TABDIMS_NRPVTG_ITEM ]; // #Composition nodes
-    raw.numRows    = tabdims[ TABDIMS_NPPVTG_ITEM ]; // #Rv (or Pg) nodes
-    raw.numCols    = 5; // [ Rv, 1/B, 1/(B*mu), d(1/B)/dRv, d(1/(B*mu))/dRv ]
-    raw.numTables  = tabdims[ TABDIMS_NTPVTG_ITEM ]; // # PVTG tables
+    const auto numRv = tabdims[ TABDIMS_NRPVTG_ITEM ]; // #Composition (Rv) nodes
+    const auto numPg = tabdims[ TABDIMS_NPPVTG_ITEM ]; // #Pg nodes
+
+    raw.numCols   = 5; // [ x, 1/B, 1/(B*mu), d(1/B)/dx, d(1/(B*mu))/dx ]
+    raw.numTables = tabdims[ TABDIMS_NTPVTG_ITEM ]; // # PV{D,T}G tables
+
+    const auto& lh = init.keywordData<bool>(LOGIHEAD_KW);
+
+    if (lh[ LOGIHEAD_RV_INDEX ]) {
+        // Vaporised oil flag set => Wet Gas.
+        //
+        // Inner dimension (number of rows per sub-table) is number of
+        // composition nodes.  Number of primary keys (outer dimension,
+        // number of sub-tables per PVTG table) is number of pressure nodes.
+
+        raw.numPrimary = numPg;
+        raw.numRows    = numRv;
+    }
+    else {
+        // Vaporised oil flag NOT set => Dry Gas.
+        //
+        // Inner dimension (number of rows per sub-table) is number of
+        // pressure nodes.  Number of primary keys (outer dimension, number
+        // of sub-tables per PVDG table) is one.
+
+        raw.numPrimary = 1;
+        raw.numRows    = numPg;
+    }
 
     // Extract Primary Key (Pg)
     {
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.cpp
index 5161784815..b572ccb35a 100644
--- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.cpp
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.cpp
@@ -40,6 +40,48 @@
 #include <ert/ecl/ecl_kw_magic.h>
 
 namespace {
+    std::function<double(double)>
+    fvf(const Opm::ECLUnits::UnitSystem& usys)
+    {
+        const auto scale = usys.reservoirVolume()
+                / usys.surfaceVolumeLiquid();
+
+        return [scale](const double x) -> double
+        {
+            return Opm::unit::convert::from(x, scale);
+        };
+    }
+
+    std::function<double(double)>
+    viscosity(const Opm::ECLUnits::UnitSystem& usys)
+    {
+        const auto scale = usys.viscosity();
+
+        return [scale](const double x) -> double
+        {
+            return Opm::unit::convert::from(x, scale);
+        };
+    }
+
+    ::Opm::ECLPVT::ConvertUnits pvcdoUnitConverter(const int usys)
+    {
+        using ToSI = ::Opm::ECLPVT::CreateUnitConverter::ToSI;
+
+        const auto u = ::Opm::ECLUnits::createUnitSystem(usys);
+
+        // [ Pref, Bo, Co, mu_o, Cv ]
+
+        return ::Opm::ECLPVT::ConvertUnits {
+            ToSI::pressure(*u),
+            {
+                fvf(*u),
+                ToSI::compressibility(*u),
+                viscosity(*u),
+                ToSI::compressibility(*u)
+            }
+        };
+    }
+
     ::Opm::ECLPVT::ConvertUnits
     deadOilUnitConverter(const ::Opm::ECLUnits::UnitSystem& usys)
     {
@@ -101,6 +143,142 @@ public:
     virtual std::unique_ptr<PVxOBase> clone() const = 0;
 };
 
+// =====================================================================
+
+class DeadOilConstCompr : public PVxOBase
+{
+public:
+    using ElemIt = std::vector<double>::const_iterator;
+    using ConvertUnits = ::Opm::ECLPVT::ConvertUnits;
+
+    DeadOilConstCompr(ElemIt               xBegin,
+              ElemIt               xEnd,
+              const ConvertUnits&  convert,
+              std::vector<ElemIt>& colIt);
+
+    virtual std::vector<double>
+    formationVolumeFactor(const std::vector<double>& /* rs */,
+                          const std::vector<double>& po) const
+    {
+        return this->evaluate(po, [this](const double p) -> double
+            {
+                // 1 / (1 / B)
+                return 1.0 / this->recipFvf(p);
+            });
+    }
+
+    virtual std::vector<double>
+    viscosity(const std::vector<double>& /*rs*/,
+              const std::vector<double>& po) const
+    {
+        return this->evaluate(po, [this](const double p) -> double
+            {
+                // (1 / B) / (1 / (B * mu))
+                return this->recipFvf(p) / this->recipFvfVisc(p);
+            });
+    }
+
+    virtual std::vector<Opm::ECLPVT::PVTGraph>
+    getPvtCurve(const Opm::ECLPVT::RawCurve /*curve*/) const
+    {
+        return {};
+    }
+
+    virtual std::unique_ptr<PVxOBase> clone() const
+    {
+        return std::unique_ptr<PVxOBase>(new DeadOilConstCompr(*this));
+    }
+
+    double& surfaceMassDensity()
+    {
+        return this->rhoS_;
+    }
+
+    double surfaceMassDensity() const
+    {
+        return this->rhoS_;
+    }
+
+private:
+    double po_ref_       { 1.0 };
+    double fvf_     { 1.0 }; // B
+    double visc_ { 1.0 }; // mu
+    double Co_           { 1.0 };
+    double cv_     { 0.0 }; // Cv
+    double rhoS_         { 0.0 };
+
+    double recipFvf(const double po) const
+    {
+        const auto x = this->Co_ * (po - this->po_ref_);
+
+        return this->exp(x) / this->fvf_ ;
+    }
+
+    double recipFvfVisc(const double po) const
+    {
+        const auto y = (this->Co_ - this->cv_) * (po - this->po_ref_);
+
+        return  this->exp(y) / (this->fvf_ * this->visc_);
+    }
+
+    double exp(const double x) const
+    {
+        return 1.0 + x*(1.0 + x/2.0);
+    }
+
+    template <class CalcQuant>
+    std::vector<double>
+    evaluate(const std::vector<double>& po,
+             CalcQuant&&                calculate) const
+    {
+        auto q = std::vector<double>{};
+        q.reserve(po.size());
+
+        for (const auto& poi : po) {
+            q.push_back(calculate(poi));
+        }
+
+        return q;
+    }
+};
+
+DeadOilConstCompr::DeadOilConstCompr(ElemIt               xBegin,
+                     ElemIt               xEnd,
+                     const ConvertUnits&  convert,
+                     std::vector<ElemIt>& colIt)
+{
+    // Recall: Table is
+    //
+    //    [ Po, Bo, Co, mu_o, Cv ]
+    //
+    // xBegin is Pw, colIt is remaining four columns.
+
+    this->fvf_     = convert.column[0](*colIt[0]); // Bo
+    this->Co_           = convert.column[1](*colIt[1]); // Co
+    this->visc_ = convert.column[2](*colIt[2]); // mu_o
+    this->cv_     = convert.column[3](*colIt[3]); // Cw - Cv
+
+    // Honour requirement that constructor advances column iterators.
+    const auto N = std::distance(xBegin, xEnd);
+    for (auto& it : colIt) { it += N; }
+
+    if (! (std::abs(*xBegin) < 1.0e20)) {
+          throw std::invalid_argument {
+            "Invalid Input PVCDO Table"
+        };
+    }
+
+    this->po_ref_ = convert.indep(*xBegin);
+
+    ++xBegin;
+    if (std::abs(*xBegin) < 1.0e20) {
+        throw std::invalid_argument {
+          "Probably Invalid Input PVCDO Table"
+      };
+    }
+}
+
+
 // =====================================================================
 
 class DeadOil : public PVxOBase
@@ -227,6 +405,7 @@ private:
 namespace {
     std::vector<std::unique_ptr<PVxOBase>>
     createDeadOil(const ::Opm::ECLPropTableRawData& raw,
+                  const bool   const_compr,
                   const int                         usys)
     {
         using PVTInterp = std::unique_ptr<PVxOBase>;
@@ -235,6 +414,16 @@ namespace {
         assert ((raw.numPrimary == 1) &&
                 "Can't Create Dead Oil Function From Live Oil Table");
 
+        if (const_compr) {
+            const auto cvrt = pvcdoUnitConverter(usys);
+
+            return ::Opm::MakeInterpolants<PVTInterp>::fromRawData(raw,
+                [&cvrt](ElmIt xBegin, ElmIt xEnd, std::vector<ElmIt>& colIt)
+            {
+                return PVTInterp{ new DeadOilConstCompr(xBegin, xEnd, cvrt, colIt) };
+            });
+        }
+
         const auto cvrt = deadOilUnitConverter(usys);
 
         return ::Opm::MakeInterpolants<PVTInterp>::fromRawData(raw,
@@ -313,6 +502,7 @@ namespace {
 
     std::vector<std::unique_ptr<PVxOBase>>
     createPVTFunction(const ::Opm::ECLPropTableRawData& raw,
+                      const bool   const_compr,
                       const int                         usys)
     {
         if (raw.numPrimary == 0) {
@@ -344,7 +534,7 @@ namespace {
         }
 
         if (raw.numPrimary == 1) {
-            return createDeadOil(raw, usys);
+            return createDeadOil(raw, const_compr, usys);
         }
 
         return createLiveOil(raw, usys);
@@ -379,6 +569,7 @@ private:
 public:
     Impl(const ECLPropTableRawData& raw,
          const int                  usys,
+         const bool     const_compr,
          std::vector<double>        rhoS);
 
     Impl(const Impl& rhs);
@@ -420,8 +611,9 @@ private:
 Opm::ECLPVT::Oil::Impl::
 Impl(const ECLPropTableRawData& raw,
      const int                  usys,
+     const bool     const_compr,
      std::vector<double>        rhoS)
-    : eval_(createPVTFunction(raw, usys))
+    : eval_(createPVTFunction(raw, const_compr, usys))
     , rhoS_(std::move(rhoS))
 {}
 
@@ -504,8 +696,9 @@ Opm::ECLPVT::Oil::Impl::validateRegIdx(const RegIdx region) const
 
 Opm::ECLPVT::Oil::Oil(const ECLPropTableRawData& raw,
                       const int                  usys,
+                      const bool const_compr,
                       std::vector<double>        rhoS)
-    : pImpl_(new Impl(raw, usys, std::move(rhoS)))
+    : pImpl_(new Impl(raw, usys, const_compr, std::move(rhoS)))
 {}
 
 Opm::ECLPVT::Oil::~Oil()
@@ -583,15 +776,39 @@ fromECLOutput(const ECLInitFileData& init)
         return OPtr{};
     }
 
+    const auto& lh = init.keywordData<bool>(LOGIHEAD_KW);
+    const int LOGIHEAD_CONST_COMPR_INDEX = 38;
+    const bool is_const_compr = lh[LOGIHEAD_CONST_COMPR_INDEX];
+
     auto raw = ::Opm::ECLPropTableRawData{};
 
     const auto& tabdims = init.keywordData<int>("TABDIMS");
     const auto& tab     = init.keywordData<double>("TAB");
 
-    raw.numPrimary = tabdims[ TABDIMS_NRPVTO_ITEM ]; // #Rs nodes/full table
-    raw.numRows    = tabdims[ TABDIMS_NPPVTO_ITEM ]; // #Po nodes/sub-table
-    raw.numCols    = 5; // [ Po, 1/B, 1/(B*mu), d(1/B)/dPo, d(1/(B*mu))/dPo ]
-    raw.numTables  = tabdims[ TABDIMS_NTPVTO_ITEM ]; // # PVTO tables
+    const auto numRs = tabdims[ TABDIMS_NRPVTO_ITEM ]; // #Rs nodes/full table
+
+    // Inner dimension (number of rows per sub-table) is number of pressure
+    // nodes for both live oil and dead oil cases.
+    raw.numRows   = tabdims[ TABDIMS_NPPVTO_ITEM ]; // #Po nodes/sub-table
+    raw.numCols   = 5; // [ Po, 1/B, 1/(B*mu), d(1/B)/dPo, d(1/(B*mu))/dPo ]
+    raw.numTables = tabdims[ TABDIMS_NTPVTO_ITEM ]; // # PVTO tables
+
+    if (lh[ LOGIHEAD_RS_INDEX ]) {
+        // Dissolved gas flag set => Live Oil.
+        //
+        // Number of primary keys (outer dimension, number of sub-tables per
+        // PVTO table) is number of composition nodes.
+
+        raw.numPrimary = numRs;
+    }
+    else {
+        // Dissolved gas flag NOT set => Dead Oil.
+        //
+        // Number of primary keys (outer dimension, number of sub-tables per
+        // PVDO table) is one.
+
+        raw.numPrimary = 1;
+    }
 
     // Extract Primary Key (Rs)
     {
@@ -617,6 +834,6 @@ fromECLOutput(const ECLInitFileData& init)
     auto rhoS = surfaceMassDensity(init, ECLPhaseIndex::Liquid);
 
     return OPtr{
-        new Oil(raw, ih[ INTEHEAD_UNIT_INDEX ], std::move(rhoS))
+        new Oil(raw, ih[ INTEHEAD_UNIT_INDEX ], is_const_compr, std::move(rhoS))
     };
 }
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.hpp
index 78235de9b0..da8adf52b2 100644
--- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.hpp
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.hpp
@@ -40,8 +40,10 @@ namespace Opm { namespace ECLPVT {
     public:
         /// Constructor.
         ///
-        /// \param[in] raw Raw tabulated data.  Must correspond to the PVTO
-        ///    vector of an ECL INIT file..
+        /// \param[in] raw Raw tabulated data.  Must correspond to either
+        ///    the PVTO vector of an ECL INIT file (tabulated live oil),
+        ///    to PVDO data (tabulated dead oil) or to PVCDO data (constant
+        ///    compressibility dead oil).
         ///
         /// \param[in] usys Unit system convention of the result set from
         ///    which \p raw was extracted.  Must correspond to item 3 of the
@@ -52,6 +54,7 @@ namespace Opm { namespace ECLPVT {
         ///    \endcode.
         Oil(const ECLPropTableRawData& raw,
             const int                  usys,
+            const bool                 const_compr,
             std::vector<double>        rhoS);
 
         /// Destructor.
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.cpp
index 73eff37197..78a477fddd 100644
--- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.cpp
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.cpp
@@ -112,9 +112,22 @@ namespace Opm
             return ECLUnits::createUnitSystem(unit_code)->reservoirRate();
         }
 
+        double surfaceRateGasUnit(const int unit_code)
+        {
+            return ( ECLUnits::createUnitSystem(unit_code)->surfaceVolumeGas() )/
+                    ( ECLUnits::createUnitSystem(unit_code)->time() );
+        }
 
+        double pressureUnit(const int unit_code)
+        {
+             return ECLUnits::createUnitSystem(unit_code)->pressure();
+        }
 
-
+        double surfaceRateLiquidUnit(const int unit_code)
+        {
+            return ( ECLUnits::createUnitSystem(unit_code)->surfaceVolumeLiquid() )/
+                    ( ECLUnits::createUnitSystem(unit_code)->time() );
+        }
         // Return input string with spaces stripped of the right end.
         std::string trimSpacesRight(const std::string& s)
         {
@@ -187,6 +200,9 @@ namespace Opm
             return {};
         }
         const double qr_unit = resRateUnit(ih.unit);
+        const double qGs_unit = surfaceRateGasUnit(ih.unit);
+        const double qLs_unit = surfaceRateLiquidUnit(ih.unit);
+        const double pressure_unit = pressureUnit(ih.unit);
 
         // Load well topology and flow rates.
         auto zwel = restart.keywordData<std::string>(ZWEL_KW, gridName);
@@ -214,6 +230,16 @@ namespace Opm
             }
             // Otherwise: add data for this well.
             WellData wd;
+            // add well rates
+            int xwel_offset = well * ih.nxwel;
+            wd.qOs  = -unit::convert::from(xwel[ 0 + xwel_offset], qLs_unit);
+            wd.qWs  = -unit::convert::from(xwel[ 1 + xwel_offset], qLs_unit);
+            wd.qGs  = -unit::convert::from(xwel[ 2 + xwel_offset], qGs_unit);
+            wd.lrat = -unit::convert::from(xwel[ 3 + xwel_offset], qLs_unit);
+            wd.qr   = -unit::convert::from(xwel[ 4 + xwel_offset], qr_unit);
+            wd.bhp  = unit::convert::from(xwel[ 6 + xwel_offset], pressure_unit);
+
+
             wd.name = trimSpacesRight(zwel[well * ih.nzwel]);
             const bool is_producer = (iwel[well * ih.niwel + IWEL_TYPE_INDEX] == IWEL_TYPE_PRODUCER);
             wd.is_injector_well = !is_producer;
@@ -229,6 +255,9 @@ namespace Opm
                                    icon[icon_offset + ICON_K_INDEX] - 1 };
                 // Note: taking the negative input, to get inflow rate.
                 completion.reservoir_inflow_rate = -unit::convert::from(xcon[xcon_offset + XCON_QR_INDEX], qr_unit);
+                completion.qOs  =   -unit::convert::from(xcon[xcon_offset + 0 ], qLs_unit);
+                completion.qWs  =   -unit::convert::from(xcon[xcon_offset + 1] , qLs_unit);
+                completion.qGs  =   -unit::convert::from(xcon[xcon_offset + 2 ], qGs_unit);
                 if (disallow_crossflow_) {
                     // Add completion only if not cross-flowing (injecting producer or producing injector).
                     if ((completion.reservoir_inflow_rate < 0.0) == is_producer) {
diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.hpp
index 94e24c65ec..25f584d29c 100644
--- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.hpp
+++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLWellSolution.hpp
@@ -45,11 +45,23 @@ namespace Opm
         {
             std::string name;
             bool is_injector_well;
+
+            double qOs;   // Well oil surface volume rate.
+            double qWs;   // Well water surface volume rate.
+            double qGs;   // Well gas surface volume rate.
+            double lrat;  // Well liquid (oil + water) surface volume rate.
+            double bhp;   // Well bottom hole pressure.
+            double qr;    // Well total reservoir volume rate.
+
             struct Completion
             {
-                std::string gridName;          // empty for main grid, otherwise LGR grid.
+                std::string gridName;          // Empty for main grid, otherwise LGR grid.
                 std::array<int, 3> ijk;        // Cartesian location in grid.
-                double reservoir_inflow_rate;  // Total fluid rate in SI (m^3/s).
+                double reservoir_inflow_rate;  // Total reservoir volume fluid rate in SI (m^3/s).
+                double qOs;                    // Completion oil surface volute rate.
+                double qWs;                    // Completion water surface volute rate.
+                double qGs;                    // Completion gas surface volute rate.
+
             };
             std::vector<Completion> completions;
         };