diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index ea1783aaa..6ccb66117 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -211,3 +211,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/BlackoilWellModel.hpp
   opm/simulators/wells/BlackoilWellModel_impl.hpp
   )
+
+list (APPEND EXAMPLE_SOURCE_FILES
+  examples/printvfp.cpp
+  )
diff --git a/examples/printvfp.cpp b/examples/printvfp.cpp
new file mode 100644
index 000000000..eb8ce71a5
--- /dev/null
+++ b/examples/printvfp.cpp
@@ -0,0 +1,129 @@
+/*
+  Copyright 2020 Equinor.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Python/Python.hpp>
+#include <opm/parser/eclipse/Units/Units.hpp>
+#include <opm/simulators/wells/VFPProperties.hpp>
+#include <opm/simulators/wells/VFPInjProperties.hpp>
+#include <opm/simulators/wells/VFPProdProperties.hpp>
+
+#include <iostream>
+#include <iomanip>
+
+using namespace Opm;
+
+struct Setup
+{
+    using VFP = VFPProperties<VFPInjProperties, VFPProdProperties>;
+    std::unique_ptr<const EclipseState> ecl_state;
+    std::shared_ptr<Python> python;
+    std::unique_ptr<const Schedule> schedule;
+    std::unique_ptr<SummaryState> summary_state;
+    std::unique_ptr<VFP> vfp_properties;
+
+    Setup(const std::string& file)
+    {
+        Parser parser;
+        auto deck = parser.parseFile(file);
+        ecl_state.reset(new EclipseState(deck) );
+        {
+          const TableManager table( deck );
+          const Runspec runspec(deck);
+          python = std::make_shared<Python>();
+          schedule.reset( new Schedule(deck, *ecl_state, python));
+          summary_state.reset( new SummaryState(std::chrono::system_clock::from_time_t(schedule->getStartTime())));
+        }
+        const int step = 0;
+        vfp_properties = std::make_unique<VFP>(schedule->getVFPInjTables(step), schedule->getVFPProdTables(step));
+    };
+};
+
+
+
+
+double computeBhp(const VFPProdTable& table,
+                  const double flo,
+                  const double thp,
+                  const double wfr,
+                  const double gfr,
+                  const double alq)
+{
+
+    // First, find the values to interpolate between.
+    // Assuming positive flo here!
+    assert(flo > 0.0);
+    auto flo_i = detail::findInterpData(flo, table.getFloAxis());
+    auto thp_i = detail::findInterpData(thp, table.getTHPAxis()); // assume constant
+    auto wfr_i = detail::findInterpData(wfr, table.getWFRAxis());
+    auto gfr_i = detail::findInterpData(gfr, table.getGFRAxis());
+    auto alq_i = detail::findInterpData(alq, table.getALQAxis()); //assume constant
+
+    return detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value;
+}
+
+
+
+
+int main(int argc, char** argv)
+{
+    if (argc < 2) {
+        return EXIT_FAILURE;
+    }
+    Setup setup(argv[1]);
+
+//    const int table_id = 1;
+    const int table_id = 4;
+    const double wct = 0.0;
+    const double gor = 35.2743;
+    const double alq = 0.0;
+    const int n = 51;
+    const double m3pd = unit::cubic(unit::meter)/unit::day;
+    const double rate_min = 20.0 * m3pd;
+    const double rate_max = 2000.0 * m3pd;
+    // const double thp = 32.1744 * unit::barsa;
+    // const double thp = 10.0 * unit::barsa;
+    const double thp_min = 10.0 * unit::barsa;
+    const double thp_max = 35.0 * unit::barsa;
+    std::vector<double> rates(n);
+    std::vector<double> thps(n);
+    for (int ii = 0; ii < n; ++ii) {
+        const double q = double(ii) / double(n-1);
+        rates[ii] = (1.0 - q) * rate_min + q * rate_max;
+        thps[ii] = (1.0 - q) * thp_min + q * thp_max;
+    }
+
+    const VFPProdTable& table = *(setup.vfp_properties->getProd()->getTable(table_id));
+    std::cout.precision(12);
+    for (double rate : rates) {
+        for (double thp : thps) {
+            const double bhp = computeBhp(table, rate, thp, wct, gor, alq);
+            std::cout //<< std::setw(18) << unit::convert::to(rate, m3pd)
+                      //<< std::setw(18) << unit::convert::to(thp, unit::barsa)
+                      << std::setw(18) << unit::convert::to(bhp, unit::barsa) << '\n';
+        }
+    }
+
+    return EXIT_SUCCESS;
+}