From b1de0bf3382cd39152338ecb602b442d8c4c6d1d Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Mon, 1 Dec 2014 18:15:20 +0100
Subject: [PATCH 1/4] unit system: add temperature

this requires the possibility of specifying an offset for the SI
conversion because Eclipse in its eternal wisdom chooses to specify
temperatures using degrees Celsius and degrees Fahrenheit instead of
using Kelvin an Rankine...
---
 opm/parser/eclipse/Deck/DeckDoubleItem.cpp    | 13 +++-----
 opm/parser/eclipse/Deck/DeckFloatItem.cpp     | 13 +++-----
 .../eclipse/Units/ConversionFactors.hpp       | 17 ++++++++++
 opm/parser/eclipse/Units/Dimension.cpp        | 32 ++++++++++++++++---
 opm/parser/eclipse/Units/Dimension.hpp        | 12 +++++--
 opm/parser/eclipse/Units/UnitSystem.cpp       | 16 ++++++++--
 opm/parser/eclipse/Units/UnitSystem.hpp       |  2 +-
 opm/parser/eclipse/Units/tests/UnitTests.cpp  | 10 ++++++
 8 files changed, 88 insertions(+), 27 deletions(-)

diff --git a/opm/parser/eclipse/Deck/DeckDoubleItem.cpp b/opm/parser/eclipse/Deck/DeckDoubleItem.cpp
index 9445fe0bb..7ca21a9f6 100644
--- a/opm/parser/eclipse/Deck/DeckDoubleItem.cpp
+++ b/opm/parser/eclipse/Deck/DeckDoubleItem.cpp
@@ -46,15 +46,10 @@ namespace Opm {
                 return;
             }
             m_SIdata.resize( m_data.size() );
-            if (m_dimensions.size() == 1) {
-                double SIfactor = m_dimensions[0]->getSIScaling();
-                std::transform( m_data.begin() , m_data.end() , m_SIdata.begin() , std::bind1st(std::multiplies<double>(),SIfactor));
-            } else {
-                for (size_t index=0; index < m_data.size(); index++) {
-                    size_t dimIndex = (index % m_dimensions.size());
-                    double SIfactor = m_dimensions[dimIndex]->getSIScaling();
-                    m_SIdata[index] = m_data[index] * SIfactor;
-                }
+
+            for (size_t index=0; index < m_data.size(); index++) {
+                size_t dimIndex = (index % m_dimensions.size());
+                m_SIdata[index] = m_dimensions[dimIndex]->convertRawToSi(m_data[index]);
             }
         } else
             throw std::invalid_argument("No dimension has been set for item:" + name() + " can not ask for SI data");
diff --git a/opm/parser/eclipse/Deck/DeckFloatItem.cpp b/opm/parser/eclipse/Deck/DeckFloatItem.cpp
index d68a0f7c3..517e145e9 100644
--- a/opm/parser/eclipse/Deck/DeckFloatItem.cpp
+++ b/opm/parser/eclipse/Deck/DeckFloatItem.cpp
@@ -59,15 +59,10 @@ namespace Opm {
                 return;
             }
             m_SIdata.resize( m_data.size() );
-            if (m_dimensions.size() == 1) {
-                float SIfactor = m_dimensions[0]->getSIScaling();
-                std::transform( m_data.begin() , m_data.end() , m_SIdata.begin() , std::bind1st(std::multiplies<float>(),SIfactor));
-            } else {
-                for (size_t index=0; index < m_data.size(); index++) {
-                    size_t dimIndex = (index % m_dimensions.size());
-                    float SIfactor = m_dimensions[dimIndex]->getSIScaling();
-                    m_SIdata[index] = m_data[index] * SIfactor;
-                }
+
+            for (size_t index=0; index < m_data.size(); index++) {
+                size_t dimIndex = (index % m_dimensions.size());
+                m_SIdata[index] = m_dimensions[dimIndex]->convertRawToSi(m_data[index]);
             }
         } else
             throw std::invalid_argument("No dimension has been set for item:" + name() + " can not ask for SI data");
diff --git a/opm/parser/eclipse/Units/ConversionFactors.hpp b/opm/parser/eclipse/Units/ConversionFactors.hpp
index 6e1d94c9c..55d6471c4 100644
--- a/opm/parser/eclipse/Units/ConversionFactors.hpp
+++ b/opm/parser/eclipse/Units/ConversionFactors.hpp
@@ -143,6 +143,19 @@ namespace Opm {
             const double psia   = lbf / square(inch);
             /// @}
 
+            /// \name Temperature. This one is more complicated
+            /// because the unit systems used by Eclipse (i.e. degrees
+            /// Celsius and degrees Fahrenheit require to add or
+            /// subtract an offset for the conversion between from/to
+            /// Kelvin
+            /// @{
+            const double degCelsius = 1.0; // scaling factor °C -> K
+            const double degCelsiusOffset = 273.15; // offset for the °C -> K conversion
+
+            const double degFahrenheit = 5.0/9; // scaling factor °F -> K
+            const double degFahrenheitOffset = 255.37; // offset for the °C -> K conversion
+            /// @}
+
             /// \name Viscosity
             /// @{
             const double Pas   = Pascal * second; // == 1
@@ -179,6 +192,8 @@ namespace Opm {
         using namespace details::prefix;
         using namespace details::unit;
         const double Pressure             = barsa;
+        const double Temperature          = degCelsius;
+        const double TemperatureOffset    = degCelsiusOffset;
         const double Length               = meter;
         const double Time                 = day;
         const double Mass                 = kilogram;
@@ -198,6 +213,8 @@ namespace Opm {
         using namespace details::prefix;
         using namespace details::unit;
         const double Pressure             = psia;
+        const double Temperature          = degFahrenheit;
+        const double TemperatureOffset    = degFahrenheitOffset;
         const double Length               = feet;
         const double Time                 = day;
         const double Mass                 = pound;
diff --git a/opm/parser/eclipse/Units/Dimension.cpp b/opm/parser/eclipse/Units/Dimension.cpp
index fe325741c..0e54799e5 100644
--- a/opm/parser/eclipse/Units/Dimension.cpp
+++ b/opm/parser/eclipse/Units/Dimension.cpp
@@ -29,7 +29,7 @@ namespace Opm {
         
     }
 
-    Dimension::Dimension(const std::string& name, double SIfactor) 
+    Dimension::Dimension(const std::string& name, double SIfactor, double SIoffset)
     {
         for (auto iter = name.begin(); iter != name.end(); ++iter) {
             if (!isalpha(*iter) && (*iter) != '1')
@@ -37,9 +37,9 @@ namespace Opm {
         }
         m_name = name;
         m_SIfactor = SIfactor;
+        m_SIoffset = SIoffset;
     }
 
-    
     double Dimension::getSIScaling() const {
         if (!std::isfinite(m_SIfactor))
             throw std::logic_error("The DeckItem contains a field with a context dependent unit. "
@@ -47,15 +47,39 @@ namespace Opm {
         return m_SIfactor;
     }
 
+    double Dimension::getSIOffset() const {
+        return m_SIoffset;
+    }
+
+    double Dimension::convertRawToSi(double rawValue) const {
+        if (!std::isfinite(m_SIfactor))
+            throw std::logic_error("The DeckItem contains a field with a context dependent unit. "
+                                   "Use getRawDoubleData() and convert the returned value manually!");
+
+        return rawValue*m_SIfactor + m_SIoffset;
+    }
+
+    double Dimension::convertSiToRaw(double siValue) const {
+        if (!std::isfinite(m_SIfactor))
+            throw std::logic_error("The DeckItem contains a field with a context dependent unit. "
+                                   "Use getRawDoubleData() and convert the returned value manually!");
+
+        return (siValue - m_SIoffset)/m_SIfactor;
+    }
+
     const std::string& Dimension::getName() const {
         return m_name;
     }
 
+    // only dimensions with zero offset are compositable...
+    bool Dimension::isCompositable() const
+    { return m_SIoffset == 0.0; }
 
-    Dimension * Dimension::newComposite(const std::string& dim , double SIfactor) {
+    Dimension * Dimension::newComposite(const std::string& dim , double SIfactor, double SIoffset) {
         Dimension * dimension = new Dimension();
         dimension->m_name = dim;
         dimension->m_SIfactor = SIfactor;
+        dimension->m_SIoffset = SIoffset;
 
         return dimension;
     }
@@ -64,7 +88,7 @@ namespace Opm {
     bool Dimension::equal(const Dimension& other) const {
         if (m_name != other.m_name)
             return false;
-        if (m_SIfactor == other.m_SIfactor)
+        if (m_SIfactor == other.m_SIfactor && m_SIoffset == other.m_SIoffset)
             return true;
         if (std::isnan(m_SIfactor) && std::isnan(other.m_SIfactor))
             return true;
diff --git a/opm/parser/eclipse/Units/Dimension.hpp b/opm/parser/eclipse/Units/Dimension.hpp
index 382dc3c13..acc93094f 100644
--- a/opm/parser/eclipse/Units/Dimension.hpp
+++ b/opm/parser/eclipse/Units/Dimension.hpp
@@ -26,16 +26,24 @@ namespace Opm {
 
     class Dimension {
     public:
-        Dimension(const std::string& name, double SI_factor);
+        Dimension(const std::string& name, double SIfactor, double SIoffset = 0.0);
+
         double getSIScaling() const;
+        double getSIOffset() const;
+
+        double convertRawToSi(double rawValue) const;
+        double convertSiToRaw(double siValue) const;
+
         bool equal(const Dimension& other) const;
         const std::string& getName() const;
-        static Dimension * newComposite(const std::string& dim , double SIfactor);
+        bool isCompositable() const;
+        static Dimension * newComposite(const std::string& dim, double SIfactor, double SIoffset = 0.0);
 
     private:
         Dimension();
         std::string m_name;
         double m_SIfactor;
+        double m_SIoffset;
     };
 }
 
diff --git a/opm/parser/eclipse/Units/UnitSystem.cpp b/opm/parser/eclipse/Units/UnitSystem.cpp
index 87b8d6dc4..6ad85b233 100644
--- a/opm/parser/eclipse/Units/UnitSystem.cpp
+++ b/opm/parser/eclipse/Units/UnitSystem.cpp
@@ -66,8 +66,8 @@ namespace Opm {
         m_dimensions.insert( std::make_pair(dimension->getName() , dimension));
     }
 
-    void UnitSystem::addDimension(const std::string& dimension , double SI_factor) {
-        std::shared_ptr<const Dimension> dim( new Dimension(dimension , SI_factor) );
+    void UnitSystem::addDimension(const std::string& dimension , double SIfactor, double SIoffset) {
+        std::shared_ptr<const Dimension> dim( new Dimension(dimension , SIfactor, SIoffset) );
         addDimension(dim);
     }
 
@@ -82,6 +82,13 @@ namespace Opm {
         double SIfactor = 1.0;
         for (auto iter = dimensionList.begin(); iter != dimensionList.end(); ++iter) {
             std::shared_ptr<const Dimension> dim = getDimension( *iter );
+
+            // all constituing dimension must be compositable. The
+            // only exception is if there is the "composite" dimension
+            // consists of exactly a single atomic dimension...
+            if (dimensionList.size() > 1 && !dim->isCompositable())
+                throw std::invalid_argument("Composite dimensions currently cannot require a conversion offset");
+
             SIfactor *= dim->getSIScaling();
         }
         return std::shared_ptr<Dimension>(Dimension::newComposite( dimension , SIfactor ));
@@ -106,6 +113,9 @@ namespace Opm {
             boost::split(parts , dimension , boost::is_any_of("/"));
             std::shared_ptr<const Dimension> dividend = parseFactor( parts[0] );
             std::shared_ptr<const Dimension> divisor = parseFactor( parts[1] );
+
+            if (dividend->getSIOffset() != 0.0 || divisor->getSIOffset() != 0.0)
+                throw std::invalid_argument("Composite dimensions cannot currently require a conversion offset");
         
             return std::shared_ptr<Dimension>( Dimension::newComposite( dimension , dividend->getSIScaling() / divisor->getSIScaling() ));
         } else {
@@ -139,6 +149,7 @@ namespace Opm {
         
         system->addDimension("1"         , 1.0);
         system->addDimension("Pressure"  , Metric::Pressure );
+        system->addDimension("Temperature", Metric::Temperature, Metric::TemperatureOffset);
         system->addDimension("Length"    , Metric::Length);
         system->addDimension("Time"      , Metric::Time );
         system->addDimension("Mass"         , Metric::Mass );
@@ -162,6 +173,7 @@ namespace Opm {
         
         system->addDimension("1"    , 1.0);
         system->addDimension("Pressure", Field::Pressure );
+        system->addDimension("Temperature", Field::Temperature, Field::TemperatureOffset);
         system->addDimension("Length", Field::Length);
         system->addDimension("Time" , Field::Time);
         system->addDimension("Mass", Field::Mass);
diff --git a/opm/parser/eclipse/Units/UnitSystem.hpp b/opm/parser/eclipse/Units/UnitSystem.hpp
index 1908d4b2d..d71abaf3e 100644
--- a/opm/parser/eclipse/Units/UnitSystem.hpp
+++ b/opm/parser/eclipse/Units/UnitSystem.hpp
@@ -33,7 +33,7 @@ namespace Opm {
         UnitSystem(const std::string& unitSystem);
         const std::string& getName() const;
 
-        void addDimension(const std::string& dimension , double SI_factor);
+        void addDimension(const std::string& dimension, double SIfactor, double SIoffset = 0.0);
         void addDimension(std::shared_ptr<const Dimension> dimension);
         std::shared_ptr<const Dimension> getNewDimension(const std::string& dimension);
         std::shared_ptr<const Dimension> getDimension(const std::string& dimension) const;
diff --git a/opm/parser/eclipse/Units/tests/UnitTests.cpp b/opm/parser/eclipse/Units/tests/UnitTests.cpp
index aeb691650..51ec0ca28 100644
--- a/opm/parser/eclipse/Units/tests/UnitTests.cpp
+++ b/opm/parser/eclipse/Units/tests/UnitTests.cpp
@@ -41,6 +41,11 @@ BOOST_AUTO_TEST_CASE(makeComposite) {
     BOOST_CHECK_EQUAL(100 , composite->getSIScaling());
 }
 
+BOOST_AUTO_TEST_CASE(MakeCompositeInvalid) {
+    // conversion to SI temperatures requires an offset, but such
+    // composite units are (currently?) invalid
+    BOOST_CHECK_THROW(Dimension::newComposite("Length*Temperature", 100), std::logic_error);
+}
 
 BOOST_AUTO_TEST_CASE(CreateDimensionInvalidNameThrows) {
     BOOST_CHECK_THROW(Dimension(" " , 1) , std::invalid_argument);
@@ -89,11 +94,15 @@ BOOST_AUTO_TEST_CASE(UnitSystemAddDimensions) {
     UnitSystem system("Metric");
     system.addDimension("Length" , 1 );
     system.addDimension("Time" , 86400 );
+    system.addDimension("Temperature", 1.0, 273.15);
 
     std::shared_ptr<const Dimension> length = system.getDimension("Length");
     std::shared_ptr<const Dimension> time = system.getDimension("Time");
+    std::shared_ptr<const Dimension> temperature = system.getDimension("Temperature");
     BOOST_CHECK_EQUAL(1     , length->getSIScaling());
     BOOST_CHECK_EQUAL(86400 , time->getSIScaling());
+    BOOST_CHECK_EQUAL(1.0   , temperature->getSIScaling());
+    BOOST_CHECK_EQUAL(273.15, temperature->getSIOffset());
 
     system.addDimension("Length" , 0.3048);
     length = system.getDimension("Length");
@@ -123,6 +132,7 @@ static void checkSystemHasRequiredDimensions(std::shared_ptr<const UnitSystem> s
     BOOST_CHECK( system->hasDimension("Time"));
     BOOST_CHECK( system->hasDimension("Permeability"));
     BOOST_CHECK( system->hasDimension("Pressure"));
+    BOOST_CHECK( system->hasDimension("Temperature"));
 }
 
 

From fdd66d911cc873f4e3cec82222856fa54fc706f5 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Mon, 1 Dec 2014 18:24:47 +0100
Subject: [PATCH 2/4] add the keywords required to specify the initial
 reservoir temperature

E100 does not know about TEMPI, but it is wise not to use the
geothermal gradient directly, but to go over the "temperature at
cell-center" detour as specified by the E300 TEMPI keyword. The depth
vs. temperature table is given by the RTEMPVD keyword. The intend is
that if TEMPI stays unspecified the temperature is calculated at each
cell center and the result is put into TEMPI from where it gets picked
up by the simulator.

(Note that E300 says that TEMPVD is an alias for RTEMPVD even though
the possible sections of these two do not exhibit any overlap, WTF?)
---
 .../share/keywords/000_Eclipse100/R/RTEMPVD    | 18 ++++++++++++++++++
 .../share/keywords/000_Eclipse100/T/TEMP       |  1 +
 .../share/keywords/001_Eclipse300/T/TEMPI      |  8 ++++++++
 .../share/keywords/001_Eclipse300/T/TEMPVD     | 18 ++++++++++++++++++
 4 files changed, 45 insertions(+)
 create mode 100644 opm/parser/share/keywords/000_Eclipse100/R/RTEMPVD
 create mode 100644 opm/parser/share/keywords/000_Eclipse100/T/TEMP
 create mode 100644 opm/parser/share/keywords/001_Eclipse300/T/TEMPI
 create mode 100644 opm/parser/share/keywords/001_Eclipse300/T/TEMPVD

diff --git a/opm/parser/share/keywords/000_Eclipse100/R/RTEMPVD b/opm/parser/share/keywords/000_Eclipse100/R/RTEMPVD
new file mode 100644
index 000000000..b98089dcd
--- /dev/null
+++ b/opm/parser/share/keywords/000_Eclipse100/R/RTEMPVD
@@ -0,0 +1,18 @@
+{
+    "name" : "RTEMPVD",
+    "sections" : [ "PROPS", "SOLUTION" ],
+    "size" : {
+       "keyword" : "EQLDIMS" ,
+       "item" : "NTEQUL"
+    },
+    "items" : [
+      {
+         "name":"DATA",
+         "value_type":"DOUBLE",
+          "size_type" : "ALL",
+          "dimension" : ["Length","Temperature"]
+      }
+    ]
+}
+
+
diff --git a/opm/parser/share/keywords/000_Eclipse100/T/TEMP b/opm/parser/share/keywords/000_Eclipse100/T/TEMP
new file mode 100644
index 000000000..4da2ae3c6
--- /dev/null
+++ b/opm/parser/share/keywords/000_Eclipse100/T/TEMP
@@ -0,0 +1 @@
+{"name" : "TEMP", "sections" : ["RUNSPEC"]}
diff --git a/opm/parser/share/keywords/001_Eclipse300/T/TEMPI b/opm/parser/share/keywords/001_Eclipse300/T/TEMPI
new file mode 100644
index 000000000..55717c1f3
--- /dev/null
+++ b/opm/parser/share/keywords/001_Eclipse300/T/TEMPI
@@ -0,0 +1,8 @@
+{
+ "name": "TEMPI",
+ "sections": ["SOLUTION"],
+ "data": {
+   "value_type": "DOUBLE",
+   "dimension" : "Temperature"
+ }
+}
diff --git a/opm/parser/share/keywords/001_Eclipse300/T/TEMPVD b/opm/parser/share/keywords/001_Eclipse300/T/TEMPVD
new file mode 100644
index 000000000..24dc66ddc
--- /dev/null
+++ b/opm/parser/share/keywords/001_Eclipse300/T/TEMPVD
@@ -0,0 +1,18 @@
+{
+    "name" : "TEMPVD",
+    "sections" : [ "PROPS", "SOLUTION" ],
+    "size" : {
+       "keyword" : "EQLDIMS" ,
+       "item" : "NTEQUL"
+    },
+    "items" : [
+      {
+         "name":"DATA",
+         "value_type":"DOUBLE",
+          "size_type" : "ALL",
+          "dimension" : ["Length","Temperature"]
+      }
+    ]
+}
+
+

From 9ab48975872d59c10dd973735a0a3fb98dfbbe28 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Mon, 1 Dec 2014 18:33:56 +0100
Subject: [PATCH 3/4] add a table for RTEMPVD

---
 opm/parser/eclipse/CMakeLists.txt             |  1 +
 .../eclipse/EclipseState/EclipseState.cpp     | 15 ++++
 .../eclipse/EclipseState/EclipseState.hpp     |  3 +
 .../EclipseState/Tables/RtempvdTable.hpp      | 68 +++++++++++++++++++
 4 files changed, 87 insertions(+)
 create mode 100644 opm/parser/eclipse/EclipseState/Tables/RtempvdTable.hpp

diff --git a/opm/parser/eclipse/CMakeLists.txt b/opm/parser/eclipse/CMakeLists.txt
index f2a300605..417273eb5 100644
--- a/opm/parser/eclipse/CMakeLists.txt
+++ b/opm/parser/eclipse/CMakeLists.txt
@@ -182,6 +182,7 @@ EclipseState/Tables/ImptvdTable.hpp
 EclipseState/Tables/RsvdTable.hpp
 EclipseState/Tables/PvtoInnerTable.hpp
 EclipseState/Tables/RvvdTable.hpp
+EclipseState/Tables/RtempvdTable.hpp
 EclipseState/Tables/PvtgOuterTable.hpp
 #
 Utility/WconinjeWrapper.hpp
diff --git a/opm/parser/eclipse/EclipseState/EclipseState.cpp b/opm/parser/eclipse/EclipseState/EclipseState.cpp
index 2fc856b2b..1ba63e8bd 100644
--- a/opm/parser/eclipse/EclipseState/EclipseState.cpp
+++ b/opm/parser/eclipse/EclipseState/EclipseState.cpp
@@ -200,6 +200,10 @@ namespace Opm {
         return m_rvvdTables;
     }
 
+    const std::vector<RtempvdTable>& EclipseState::getRtempvdTables() const {
+        return m_rtempvdTables;
+    }
+
     const std::vector<SgofTable>& EclipseState::getSgofTables() const {
         return m_sgofTables;
     }
@@ -244,6 +248,17 @@ namespace Opm {
         //depends on the presence of the RKTRMDIR keyword...
         initRocktabTables(deck, parserLog);
 
+        // the temperature vs depth table. the problem here is that
+        // the TEMPVD (E300) and RTEMPVD (E300 + E100) keywords are
+        // synonymous, but we want to provide only a single cannonical
+        // API here, so we jump through some small hoops...
+        if (deck->hasKeyword("TEMPVD") && deck->hasKeyword("RTEMPVD"))
+            throw std::invalid_argument("The TEMPVD and RTEMPVD tables are mutually exclusive!");
+        else if (deck->hasKeyword("TEMPVD"))
+            initSimpleTables(deck, parserLog, "TEMPVD", m_rtempvdTables);
+        else if (deck->hasKeyword("RTEMPVD"))
+            initSimpleTables(deck, parserLog, "RTEMPVD", m_rtempvdTables);
+
         initFullTables(deck, parserLog, "PVTG", m_pvtgTables);
         initFullTables(deck, parserLog, "PVTO", m_pvtoTables);
    }
diff --git a/opm/parser/eclipse/EclipseState/EclipseState.hpp b/opm/parser/eclipse/EclipseState/EclipseState.hpp
index a6586ae4c..9ec740d0c 100644
--- a/opm/parser/eclipse/EclipseState/EclipseState.hpp
+++ b/opm/parser/eclipse/EclipseState/EclipseState.hpp
@@ -47,6 +47,7 @@
 #include <opm/parser/eclipse/EclipseState/Tables/RocktabTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/RsvdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/RvvdTable.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/RtempvdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/SgofTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/SwofTable.hpp>
 
@@ -104,6 +105,7 @@ namespace Opm {
         const std::vector<RocktabTable>& getRocktabTables() const;
         const std::vector<RsvdTable>& getRsvdTables() const;
         const std::vector<RvvdTable>& getRvvdTables() const;
+        const std::vector<RtempvdTable>& getRtempvdTables() const;
         const std::vector<SgofTable>& getSgofTables() const;
         const std::vector<SwofTable>& getSwofTables() const;
 
@@ -219,6 +221,7 @@ namespace Opm {
         std::vector<RocktabTable> m_rocktabTables;
         std::vector<RsvdTable> m_rsvdTables;
         std::vector<RvvdTable> m_rvvdTables;
+        std::vector<RtempvdTable> m_rtempvdTables;
         std::vector<SgofTable> m_sgofTables;
         std::vector<SwofTable> m_swofTables;
 
diff --git a/opm/parser/eclipse/EclipseState/Tables/RtempvdTable.hpp b/opm/parser/eclipse/EclipseState/Tables/RtempvdTable.hpp
new file mode 100644
index 000000000..a5db388bc
--- /dev/null
+++ b/opm/parser/eclipse/EclipseState/Tables/RtempvdTable.hpp
@@ -0,0 +1,68 @@
+/*
+  Copyright (C) 2014 by Andreas Lauser
+
+  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_PARSER_RTEMPVD_TABLE_HPP
+#define	OPM_PARSER_RTEMPVD_TABLE_HPP
+
+#include "SingleRecordTable.hpp"
+
+namespace Opm {
+    // forward declaration
+    class EclipseState;
+
+    class RtempvdTable : protected SingleRecordTable {
+        typedef SingleRecordTable ParentType;
+
+        friend class EclipseState;
+        RtempvdTable() = default;
+
+        /*!
+         * \brief Read the RTEMPVD keyword and provide some convenience
+         *        methods for it.
+         */
+        void init(Opm::DeckKeywordConstPtr keyword, int recordIdx)
+        {
+            ParentType::init(keyword,
+                             std::vector<std::string>{
+                                 "Depth",
+                                 "Temperature"
+                             },
+                             recordIdx,
+                             /*firstEntityOffset=*/0);
+
+            ParentType::checkNonDefaultable("Depth");
+            ParentType::checkMonotonic("Depth", /*isAscending=*/true);
+
+            ParentType::checkNonDefaultable("Temperature");
+        }
+
+    public:
+        using ParentType::numTables;
+        using ParentType::numRows;
+        using ParentType::numColumns;
+        using ParentType::evaluate;
+
+        const std::vector<double> &getDepthColumn() const
+        { return ParentType::getColumn(0); }
+
+        const std::vector<double> &getTemperatureColumn() const
+        { return ParentType::getColumn(1); }
+    };
+}
+
+#endif

From ad7f915a3ee1965020e9887ab4d4ac09061e245d Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Mon, 1 Dec 2014 18:34:24 +0100
Subject: [PATCH 4/4] add a grid property for TEMPI

by default it is initialized using the temperature vs depth table...
---
 .../eclipse/EclipseState/EclipseState.cpp     |  4 ++
 .../Grid/GridPropertyInitializers.hpp         | 46 +++++++++++++++++++
 2 files changed, 50 insertions(+)

diff --git a/opm/parser/eclipse/EclipseState/EclipseState.cpp b/opm/parser/eclipse/EclipseState/EclipseState.cpp
index 1ba63e8bd..01a618cfb 100644
--- a/opm/parser/eclipse/EclipseState/EclipseState.cpp
+++ b/opm/parser/eclipse/EclipseState/EclipseState.cpp
@@ -537,6 +537,7 @@ namespace Opm {
 
         double nan = std::numeric_limits<double>::quiet_NaN();
         const auto eptLookup = std::make_shared<GridPropertyEndpointTableLookupInitializer<>>(*deck, *this);
+        const auto tempLookup = std::make_shared<GridPropertyTemperatureLookupInitializer<>>(*deck, *this);
         const auto distributeTopLayer = std::make_shared<const GridPropertyPostProcessor::DistributeTopLayer>(*this);
         const auto initPORV = std::make_shared<GridPropertyPostProcessor::InitPORV>(*this);
 
@@ -678,6 +679,9 @@ namespace Opm {
             SupportedDoubleKeywordInfo( "ISWCRZ"   , eptLookup, "1" ),
             SupportedDoubleKeywordInfo( "ISWCRZ-"  , eptLookup, "1" ),
 
+            // cell temperature (E300 only, but makes a lot of sense for E100, too)
+            SupportedDoubleKeywordInfo( "TEMPI"    , tempLookup, "Temperature" ),
+
             // porosity
             SupportedDoubleKeywordInfo( "PORO"  , nan, distributeTopLayer , "1" ),
 
diff --git a/opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp b/opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp
index a05b476da..50ea60464 100644
--- a/opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp
+++ b/opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp
@@ -21,6 +21,7 @@
 
 #include <opm/parser/eclipse/EclipseState/Tables/SwofTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/SgofTable.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/RtempvdTable.hpp>
 
 #include <vector>
 #include <string>
@@ -349,6 +350,51 @@ private:
     const Deck& m_deck;
     const EclipseState& m_eclipseState;
 };
+
+// initialize the TEMPI grid property using the temperature vs depth
+// table (stemming from the TEMPVD or the RTEMPVD keyword)
+template <class EclipseState=Opm::EclipseState,
+          class Deck=Opm::Deck>
+class GridPropertyTemperatureLookupInitializer
+    : public GridPropertyBaseInitializer<double>
+{
+public:
+    GridPropertyTemperatureLookupInitializer(const Deck& deck, const EclipseState& eclipseState)
+        : m_deck(deck)
+        , m_eclipseState(eclipseState)
+    { }
+
+    void apply(std::vector<double>& values,
+               const std::string& propertyName) const
+    {
+        if (propertyName != "TEMPI")
+            throw std::logic_error("The TemperatureLookupInitializer can only be used for the initial temperature!");
+
+        if (!m_deck.hasKeyword("EQLNUM")) {
+            // if values are defaulted in the TEMPI keyword, but no
+            // EQLNUM is specified, you will get NaNs...
+            double nan = std::numeric_limits<double>::quiet_NaN();
+            std::fill(values.begin(), values.end(), nan);
+            return;
+        }
+
+        auto eclipseGrid = m_eclipseState.getEclipseGrid();
+        const std::vector<RtempvdTable>& rtempvdTables = m_eclipseState.getRtempvdTables();
+        const std::vector<int>& eqlNum = m_eclipseState.getIntGridProperty("EQLNUM")->getData();
+
+        for (int cellIdx = 0; cellIdx < eqlNum.size(); ++ cellIdx) {
+            int cellEquilNum = eqlNum[cellIdx];
+            const RtempvdTable& rtempvdTable = rtempvdTables[cellEquilNum];
+            double cellDepth = std::get<2>(eclipseGrid->getCellCenter(cellIdx));
+            values[cellIdx] = rtempvdTable.evaluate("Temperature", cellDepth);
+        }
+    }
+
+private:
+    const Deck& m_deck;
+    const EclipseState& m_eclipseState;
+};
+
 }
 
 #endif