diff --git a/ewoms/common/parametersystem.hh b/ewoms/common/parametersystem.hh
index 69123a08b0..3bb5671557 100644
--- a/ewoms/common/parametersystem.hh
+++ b/ewoms/common/parametersystem.hh
@@ -961,7 +961,6 @@ public:
         check_(Dune::className<ParamType>(), propTagName, paramName);
 #endif
 
-        typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta;
         if (errorIfNotRegistered) {
             if (ParamsMeta::registrationOpen())
                 throw std::runtime_error("Parameters can only checked after _all_ of them have "
@@ -1038,7 +1037,6 @@ private:
         check_(Dune::className<ParamType>(), propTagName, paramName);
 #endif
 
-        typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta;
         if (errorIfNotRegistered) {
             if (ParamsMeta::registrationOpen())
                 throw std::runtime_error("Parameters can only retieved after _all_ of them have "
diff --git a/ewoms/common/simulator.hh b/ewoms/common/simulator.hh
index b5c0d628af..cd342d6362 100644
--- a/ewoms/common/simulator.hh
+++ b/ewoms/common/simulator.hh
@@ -56,6 +56,7 @@ NEW_PROP_TAG(EndTime);
 NEW_PROP_TAG(RestartTime);
 NEW_PROP_TAG(InitialTimeStepSize);
 NEW_PROP_TAG(PredeterminedTimeStepsFile);
+NEW_PROP_TAG(SimulatorManageTimeStep);
 
 END_PROPERTIES
 
@@ -107,6 +108,8 @@ class Simulator
     typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
 
+    enum { simulatorManageTimeStep = GET_PROP_VALUE(TypeTag, SimulatorManageTimeStep) };
+
 public:
     // do not allow to copy simulators around
     Simulator(const Simulator& ) = delete;
@@ -866,11 +869,11 @@ public:
      * name and uses the extension <tt>.ers</tt>. (Ewoms ReStart
      * file.)  See Ewoms::Restart for details.
      */
-    void serialize()
+    void serialize(const bool atEndOfStep = true)
     {
         typedef Ewoms::Restart Restarter;
         Restarter res;
-        res.serializeBegin(*this);
+        res.serializeBegin(*this, atEndOfStep);
         if (gridView().comm().rank() == 0)
             std::cout << "Serialize to file '" << res.fileName() << "'"
                       << ", next time step size: " << timeStepSize()
@@ -882,6 +885,23 @@ public:
         res.serializeEnd();
     }
 
+    void deserializeAll(Scalar t,bool only_reservoir)
+    {
+        typedef Ewoms::Restart Restarter;
+        Restarter res;
+        res.deserializeBegin(*this, t);
+        if (gridView().comm().rank() == 0)
+            std::cout << "Deserialize file '" << res.fileName() << "'"
+                      << ", next time step size: " << timeStepSize()
+                      << "\n" << std::flush;
+        this->deserialize(res);
+        if( not(only_reservoir) ){
+            problem_->deserialize(res, false);
+        }
+         model_->deserialize(res);
+
+        res.deserializeEnd();
+    }
     /*!
      * \brief Write the time manager's state to a restart file.
      *
@@ -899,6 +919,7 @@ public:
             << episodeLength_ << " "
             << startTime_ << " "
             << time_ << " "
+            << timeStepSize_ << " "
             << timeStepIdx_ << " ";
         restarter.serializeSectionEnd();
     }
@@ -920,6 +941,7 @@ public:
             >> episodeLength_
             >> startTime_
             >> time_
+            >> timeStepSize_
             >> timeStepIdx_;
         restarter.deserializeSectionEnd();
     }
diff --git a/ewoms/disc/common/fvbasediscretization.hh b/ewoms/disc/common/fvbasediscretization.hh
index 4b8dbaa77a..28de030e5e 100644
--- a/ewoms/disc/common/fvbasediscretization.hh
+++ b/ewoms/disc/common/fvbasediscretization.hh
@@ -268,6 +268,7 @@ SET_BOOL_PROP(FvBaseDiscretization, UseVolumetricResidual, true);
 
 //! eWoms is mainly targeted at research, so experimental features are enabled by
 //! default.
+SET_BOOL_PROP(FvBaseDiscretization, SimulatorManageTimeStep, true);
 SET_BOOL_PROP(FvBaseDiscretization, EnableExperiments, true);
 
 END_PROPERTIES
@@ -759,6 +760,15 @@ public:
     bool enableStorageCache() const
     { return enableStorageCache_; }
 
+    /*!
+     * \brief Set the value of enable storage cache  
+     *
+     * Be aware that calling the *CachedStorage() methods if the storage cache is
+     * disabled will crash the program.
+     */
+    void setEnableStorageCache(bool enableStorageCache) 
+    { enableStorageCache_= enableStorageCache; }
+    
     /*!
      * \brief Retrieve an entry of the cache for the storage term.
      *
diff --git a/ewoms/disc/common/fvbaseelementcontext.hh b/ewoms/disc/common/fvbaseelementcontext.hh
index 5769493c41..e8c691903b 100644
--- a/ewoms/disc/common/fvbaseelementcontext.hh
+++ b/ewoms/disc/common/fvbaseelementcontext.hh
@@ -97,9 +97,10 @@ public:
     {
         // remember the simulator object
         simulatorPtr_ = &simulator;
-        enableStorageCache_ = EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache);
+        enableStorageCache_ = simulator.model().enableStorageCache();//EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache);
         stashedDofIdx_ = -1;
         focusDofIdx_ = -1;
+        focusTimeIdx_ = 0;
     }
 
     static void *operator new(size_t size)
@@ -260,6 +261,17 @@ public:
     void setFocusDofIndex(unsigned dofIdx)
     { focusDofIdx_ = dofIdx; }
 
+    /*!
+     * \brief Sets the time index on which the simulator is currently "focused" on
+     *
+     * I.e., in the case of automatic differentiation, all derivatives are with regard to
+     * the primary variables of that time index. Only "primary" DOFs can be
+     * focused on.
+     */
+    
+    void setFocusTimeIndex(unsigned timeIdx)
+    { focusTimeIdx_ = timeIdx; }
+
     /*!
      * \brief Returns the degree of freedom on which the simulator is currently "focused" on
      *
@@ -268,6 +280,14 @@ public:
     unsigned focusDofIndex() const
     { return focusDofIdx_; }
 
+    /*!
+     * \brief Returns the time index on which the simulator is currently "focused" on
+     *
+     * \copydetails setFocusDof()
+     */ 
+    unsigned focusTimeIndex() const
+    { return focusTimeIdx_; }
+
     /*!
      * \brief Return a reference to the simulator.
      */
@@ -581,7 +601,7 @@ protected:
 #endif
 
         dofVars_[dofIdx].priVars[timeIdx] = priVars;
-        dofVars_[dofIdx].intensiveQuantities[timeIdx].update(/*context=*/asImp_(), dofIdx, timeIdx);
+        dofVars_[dofIdx].intensiveQuantities[timeIdx].update(/*context=*/asImp_(), dofIdx, timeIdx, focusTimeIdx_);
     }
 
     IntensiveQuantities intensiveQuantitiesStashed_;
@@ -599,6 +619,7 @@ protected:
 
     int stashedDofIdx_;
     int focusDofIdx_;
+    int focusTimeIdx_;
     bool enableStorageCache_;
 };
 
diff --git a/ewoms/disc/common/fvbaseintensivequantities.hh b/ewoms/disc/common/fvbaseintensivequantities.hh
index 7ed89f2b6c..1928dfca83 100644
--- a/ewoms/disc/common/fvbaseintensivequantities.hh
+++ b/ewoms/disc/common/fvbaseintensivequantities.hh
@@ -67,7 +67,8 @@ public:
      */
     void update(const ElementContext& elemCtx,
                 unsigned dofIdx,
-                unsigned timeIdx)
+                unsigned timeIdx,
+                unsigned focusTimeIdx OPM_UNUSED)
     { extrusionFactor_ = elemCtx.problem().extrusionFactor(elemCtx, dofIdx, timeIdx); }
 
     /*!
diff --git a/ewoms/disc/common/fvbaselinearizer.hh b/ewoms/disc/common/fvbaselinearizer.hh
index 0bfe1eacf1..237df15c3e 100644
--- a/ewoms/disc/common/fvbaselinearizer.hh
+++ b/ewoms/disc/common/fvbaselinearizer.hh
@@ -143,6 +143,12 @@ public:
     {
         simulatorPtr_ = &simulator;
         eraseMatrix();
+        auto it = elementCtx_.begin();
+        const auto& endIt = elementCtx_.end();
+        for (; it != endIt; ++it){
+            delete *it;
+        }
+        elementCtx_.resize(0);
     }
 
     /*!
@@ -158,14 +164,15 @@ public:
     }
 
     /*!
-     * \brief Linearize the full system of non-linear equations.
+     * \brief Linearize the full system of non-linear equations about focusTimeIdx
+     * defualt is 0 i.e. current time normally time to be found.
      *
      * This means the spatial domain plus all auxiliary equations.
      */
-    void linearize()
+    void linearize(unsigned focusTimeIdx = 0)
     {
-        linearizeDomain();
-        linearizeAuxiliaryEquations();
+        linearizeDomain(focusTimeIdx);
+        linearizeAuxiliaryEquations(focusTimeIdx);
     }
 
     /*!
@@ -178,7 +185,7 @@ public:
      * The current state of affairs (esp. the previous and the current solutions) is
      * represented by the model object.
      */
-    void linearizeDomain()
+    void linearizeDomain(unsigned focustimeIndex)
     {
         // we defer the initialization of the Jacobian matrix until here because the
         // auxiliary modules usually assume the problem, model and grid to be fully
@@ -188,7 +195,7 @@ public:
 
         int succeeded;
         try {
-            linearize_();
+            linearize_(focustimeIndex);
             succeeded = 1;
         }
 #if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,5)
@@ -227,7 +234,7 @@ public:
      * \brief Linearize the part of the non-linear system of equations that is associated
      *        with the spatial domain.
      */
-    void linearizeAuxiliaryEquations()
+    void linearizeAuxiliaryEquations(unsigned focusTimeIdx OPM_UNUSED = 0)
     {
         // flush possible local caches into matrix structure
         jacobian_->commit();
@@ -237,7 +244,7 @@ public:
         for (unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx) {
             bool succeeded = true;
             try {
-                model.auxiliaryModule(auxModIdx)->linearize(*jacobian_, residual_);
+                model.auxiliaryModule(auxModIdx)->linearize(*jacobian_, residual_);//, focusTimeIdx);
             }
             catch (const std::exception& e) {
                 succeeded = false;
@@ -330,6 +337,7 @@ private:
             elementCtx_[threadId] = new ElementContext(simulator_());
     }
 
+    
     // Construct the BCRS matrix for the Jacobian of the residual function
     void createMatrix_()
     {
@@ -425,7 +433,7 @@ private:
     }
 
     // linearize the whole system
-    void linearize_()
+    void linearize_(unsigned focustimeindex)
     {
         resetSystem_();
 
@@ -472,8 +480,7 @@ private:
                     const Element& elem = *elemIt;
                     if (!linearizeNonLocalElements && elem.partitionType() != Dune::InteriorEntity)
                         continue;
-
-                    linearizeElement_(elem);
+                    linearizeElement_(elem, focustimeindex);
                 }
             }
             // If an exception occurs in the parallel block, it won't escape the
@@ -503,7 +510,7 @@ private:
     }
 
     // linearize an element in the interior of the process' grid partition
-    void linearizeElement_(const Element& elem)
+    void linearizeElement_(const Element& elem, unsigned focustimeindex)
     {
         unsigned threadId = ThreadManager::threadId();
 
@@ -511,6 +518,7 @@ private:
         auto& localLinearizer = model_().localLinearizer(threadId);
 
         // the actual work of linearization is done by the local linearizer class
+        elementCtx->setFocusTimeIndex(focustimeindex);
         localLinearizer.linearize(*elementCtx, elem);
 
         // update the right hand side and the Jacobian matrix
diff --git a/ewoms/disc/common/fvbaselocalresidual.hh b/ewoms/disc/common/fvbaselocalresidual.hh
index 1bf4f04052..607ebb1613 100644
--- a/ewoms/disc/common/fvbaselocalresidual.hh
+++ b/ewoms/disc/common/fvbaselocalresidual.hh
@@ -163,16 +163,17 @@ public:
         assert(residual.size() == elemCtx.numDof(/*timeIdx=*/0));
 
         residual = 0.0;
-
+        if(elemCtx.focusTimeIndex()==0){//NB: This is only valid if fluxeterm sis purely implicite, This was nesseary to not get false derivatives form wells in adjiont.
         // evaluate the flux terms
-        asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0);
-
+            asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0);
+        }
         // evaluate the storage and the source terms
-        asImp_().evalVolumeTerms_(residual, elemCtx);
 
+        asImp_().evalVolumeTerms_(residual, elemCtx);
+        if(elemCtx.focusTimeIndex()==0){//NB: This is only valid if fluxeterm sis purely implicite, This was nesseary to not get false derivatives form wells in adjiont.
         // evaluate the boundary conditions
-        asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0);
-
+            asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0);
+        }
         if (useVolumetricResidual) {
             // make the residual volume specific (i.e., make it incorrect mass per cubic
             // meter instead of total mass)
@@ -492,6 +493,7 @@ protected:
     {
         EvalVector tmp;
         EqVector tmp2;
+        EvalVector tmp2Der;// this is added to support calculation of storage term derivative focusTime!=0 for need for adjoint: always used if cachedStororage term si false
         RateVector sourceRate;
 
         tmp = 0.0;
@@ -572,24 +574,41 @@ protected:
             else {
                 // if the mass storage at the beginning of the time step is not cached,
                 // we re-calculate it from scratch.
-                tmp2 = 0.0;
-                asImp_().computeStorage(tmp2, elemCtx,  dofIdx, /*timeIdx=*/1);
-                Opm::Valgrind::CheckDefined(tmp2);
+                // if storage cache not enables we always use derivatives
+                tmp2Der = 0.0;
+                asImp_().computeStorage(tmp2Der, elemCtx,  dofIdx, /*timeIdx=*/1);
+                Opm::Valgrind::CheckDefined(tmp2Der);
             }
 
-            // Use the implicit Euler time discretization
-            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
-                tmp[eqIdx] -= tmp2[eqIdx];
-                tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize();
-
-                residual[dofIdx][eqIdx] += tmp[eqIdx];
+            if( !elemCtx.enableStorageCache() or elemCtx.focusTimeIndex()>0){
+                // Use the implicit Euler time discretization
+                for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+                    if (elemCtx.focusTimeIndex() == 0) {
+                        tmp2Der[eqIdx] = Toolbox::value(tmp2Der[eqIdx]);
+                    } else {
+                        tmp[eqIdx] = Toolbox::value(tmp[eqIdx]);
+                    }
+                    tmp[eqIdx] -= tmp2Der[eqIdx];
+                    tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize();
+                    residual[dofIdx][eqIdx] += tmp[eqIdx];
+                }
+            }else{
+                for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+                    // cached storage should only be used in forward mode
+                    assert(elemCtx.focusTimeIndex()==0);
+                    tmp[eqIdx] -= tmp2[eqIdx];
+                    tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize();
+                    residual[dofIdx][eqIdx] += tmp[eqIdx];
+                }
             }
 
             Opm::Valgrind::CheckDefined(residual[dofIdx]);
 
             // deal with the source term
-            asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0);
-
+            //assumes sourceterm is independent of prevois state
+            if(elemCtx.focusTimeIndex()==0){
+                asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0);
+            }
             // if the model uses extensive quantities in its storage term, and we use
             // automatic differention and current DOF is also not the one we currently
             // focus on, the storage term does not need any derivatives!
diff --git a/ewoms/disc/common/fvbaseprimaryvariables.hh b/ewoms/disc/common/fvbaseprimaryvariables.hh
index bce1404a3c..c479435a82 100644
--- a/ewoms/disc/common/fvbaseprimaryvariables.hh
+++ b/ewoms/disc/common/fvbaseprimaryvariables.hh
@@ -87,17 +87,16 @@ public:
      * it represents the a constant f = x_i. (the difference is that in the first case,
      * the derivative w.r.t. x_i is 1, while it is 0 in the second case.
      */
-    Evaluation makeEvaluation(unsigned varIdx, unsigned timeIdx) const
+    Evaluation makeEvaluation(unsigned varIdx, unsigned timeIdx, unsigned focusTimeIdx = 0) const
     {
-        if (std::is_same<Evaluation, Scalar>::value)
-            return (*this)[varIdx]; // finite differences
-        else {
-            // automatic differentiation
-            if (timeIdx == 0)
-                return Toolbox::createVariable((*this)[varIdx], varIdx);
-            else
-                return Toolbox::createConstant((*this)[varIdx]);
-        }
+    if (std::is_same<Evaluation, Scalar>::value)
+	return (*this)[varIdx]; // finite differences
+    else {
+        if (timeIdx == focusTimeIdx)
+            return Toolbox::createVariable((*this)[varIdx], varIdx);
+        else
+            return Toolbox::createConstant((*this)[varIdx]);
+    }
     }
 
     /*!
diff --git a/ewoms/disc/ecfv/ecfvdiscretization.hh b/ewoms/disc/ecfv/ecfvdiscretization.hh
index 9ae5805e5b..913f88a814 100644
--- a/ewoms/disc/ecfv/ecfvdiscretization.hh
+++ b/ewoms/disc/ecfv/ecfvdiscretization.hh
@@ -202,7 +202,7 @@ public:
     void deserialize(Restarter& res)
     {
         res.template deserializeEntities</*codim=*/0>(asImp_(), this->gridView_);
-        this->solution(/*timeIdx=*/1) = this->solution(/*timeIdx=*/0);
+        this->solution(/*timeIdx=*/1) = this->solution(/*timeIdx=*/0);// may wrong and only valid for restart
     }
 
 private:
diff --git a/ewoms/io/restart.hh b/ewoms/io/restart.hh
index fc0633c829..c98f17b718 100644
--- a/ewoms/io/restart.hh
+++ b/ewoms/io/restart.hh
@@ -31,9 +31,34 @@
 #include <fstream>
 #include <iostream>
 #include <sstream>
+# include <boost/filesystem.hpp>
 
 namespace Ewoms {
 
+
+// This small trio of functions is necessary to provide a generic
+// interface to getting grid names, when only some grids provide the
+// name() method.
+template <class GridView>
+auto getGridNameImp(const GridView& gridView, int)
+    -> decltype(gridView.grid().name())
+{
+    return gridView.grid().name();
+};
+template <class GridView>
+auto getGridNameImp(const GridView&, long long)
+    -> std::string
+{
+    return "GridWithNoNameMethod";
+};
+template <class GridView>
+auto getGridName(const GridView& gridView)
+    -> decltype(getGridNameImp(gridView, 0))
+{
+    return getGridNameImp(gridView, 0);
+};
+
+
 /*!
  * \brief Load or save a state of a problem to/from the harddisk.
  */
@@ -46,7 +71,7 @@ class Restart
     template <class GridView>
     static const std::string magicRestartCookie_(const GridView& gridView)
     {
-        static const std::string gridName = "blubb"; // gridView.grid().name();
+        static const std::string gridName = getGridName(gridView);
         static const int dim = GridView::dimension;
 
         int numVertices = gridView.size(dim);
@@ -69,22 +94,33 @@ class Restart
     /*!
      * \brief Return the restart file name.
      */
-    template <class GridView, class Scalar>
-    static const std::string restartFileName_(const GridView& gridView,
-                                              const std::string& outputDir,
-                                              const std::string& simName,
-                                              Scalar t)
+    template <class Simulator, class Scalar>
+    static const std::string restartFileName_(const Simulator& simulator,                                                                                     Scalar t)
     {
-        std::string dir = outputDir;
+        // Directory.
+        std::string dir = simulator.problem().outputDir();
         if (dir == ".")
             dir = "";
         else if (!dir.empty() && dir.back() != '/')
             dir += "/";
-
-        int rank = gridView.comm().rank();
+        namespace fs = boost::filesystem;
+        fs::path output_dir(dir);
+        fs::path subdir("ebos_restart");
+        output_dir = output_dir / subdir;
+        if(!(fs::exists(output_dir))){
+                fs::create_directory(output_dir);
+         }
+
+        // Filename.
+        int rank = simulator.gridView().comm().rank();
+        std::string simName = simulator.problem().name();
         std::ostringstream oss;
-        oss << dir << simName << "_time=" << t << "_rank=" << rank << ".ers";
-        return oss.str();
+        oss << simName << "_time=" << t << "_rank=" << rank << ".ers";
+        fs::path output_file(oss.str());
+
+        // Combine and return.
+        fs::path full_path = output_dir / output_file;
+        return full_path.string();
     }
 
 public:
@@ -98,13 +134,11 @@ public:
      * \brief Write the current state of the model to disk.
      */
     template <class Simulator>
-    void serializeBegin(Simulator& simulator)
+    void serializeBegin(Simulator& simulator, const bool atEndOfStep = true)
     {
         const std::string magicCookie = magicRestartCookie_(simulator.gridView());
-        fileName_ = restartFileName_(simulator.gridView(),
-                                     simulator.problem().outputDir(),
-                                     simulator.problem().name(),
-                                     simulator.time());
+        const double serializationTime = atEndOfStep ? simulator.time() + simulator.timeStepSize() : simulator.time();
+        fileName_ = restartFileName_(simulator, serializationTime);
 
         // open output file and write magic cookie
         outStream_.open(fileName_.c_str());
@@ -164,14 +198,16 @@ public:
     void serializeEnd()
     { outStream_.close(); }
 
-    /*!
+
+   /*!
      * \brief Start reading a restart file at a certain simulated
      *        time.
      */
+
     template <class Simulator, class Scalar>
     void deserializeBegin(Simulator& simulator, Scalar t)
     {
-        fileName_ = restartFileName_(simulator.gridView(), simulator.problem().outputDir(), simulator.problem().name(), t);
+        fileName_ = restartFileName_(simulator, t);
 
         // open input file and read magic cookie
         inStream_.open(fileName_.c_str());
diff --git a/ewoms/linear/matrixblock.hh b/ewoms/linear/matrixblock.hh
index cdb39a05e1..e6423b4efd 100644
--- a/ewoms/linear/matrixblock.hh
+++ b/ewoms/linear/matrixblock.hh
@@ -29,7 +29,7 @@
 #include <dune/istl/solvers.hh>
 #include <dune/istl/owneroverlapcopy.hh>
 #include <dune/istl/paamg/amg.hh>
-
+#include <dune/istl/matrixmarket.hh>
 #include <dune/common/fmatrix.hh>
 
 namespace Ewoms {
diff --git a/ewoms/linear/matrixmarket_ewoms.hh b/ewoms/linear/matrixmarket_ewoms.hh
new file mode 100644
index 0000000000..915eb08523
--- /dev/null
+++ b/ewoms/linear/matrixmarket_ewoms.hh
@@ -0,0 +1,198 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef MATRIXMARKET_EWOMS_HH
+#define MATRIXMARKET_EWOMS_HH
+#include <dune/istl/matrixmarket.hh>
+#include <ewoms/linear/matrixblock.hh>
+namespace Dune
+{
+
+
+  namespace MatrixMarketImpl
+  {
+
+
+    template<typename T, typename A, int i, int j>
+    struct mm_header_printer<BCRSMatrix<Ewoms::MatrixBlock<T,i,j>,A> >
+    {
+      static void print(std::ostream& os)
+      {
+        os<<"%%MatrixMarket matrix coordinate ";
+        os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
+      }
+    };
+
+
+    template<typename T, int i, int j>
+    struct mm_header_printer<Ewoms::MatrixBlock<T,i,j> >
+    {
+      static void print(std::ostream& os)
+      {
+        os<<"%%MatrixMarket matrix array ";
+        os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
+      }
+    };
+
+    /**
+     * @brief Metaprogram for writing the ISTL block
+     * structure header.
+     *
+     * Member function mm_block_structure_header::print(os, mat) writes
+     * the corresponding header to the specified ostream.
+     * @tparam The type of the matrix to generate the header for.
+     */
+
+    template<typename T, typename A, int i, int j>
+    struct mm_block_structure_header<BCRSMatrix<Ewoms::MatrixBlock<T,i,j>,A> >
+    {
+      typedef BCRSMatrix<Ewoms::MatrixBlock<T,i,j>,A> M;
+
+      static void print(std::ostream& os, const M&)
+      {
+        os<<"% ISTL_STRUCT blocked ";
+        os<<i<<" "<<j<<std::endl;
+      }
+    };
+
+
+    template<typename T, int i, int j>
+    struct mm_block_structure_header<Ewoms::MatrixBlock<T,i,j> >
+    {
+      typedef Ewoms::MatrixBlock<T,i,j> M;
+
+      static void print(std::ostream& os, const M& m)
+      {}
+    };
+
+
+
+    template<typename T, typename A, int brows, int bcols, typename D>
+    void readSparseEntries(Dune::BCRSMatrix<Ewoms::MatrixBlock<T,brows,bcols>,A>& matrix,
+                           std::istream& file, std::size_t entries,
+                           const MMHeader& mmHeader, const D&)
+    {
+      typedef Dune::BCRSMatrix<Ewoms::MatrixBlock<T,brows,bcols>,A> Matrix;
+      // First path
+      // store entries together with column index in a speparate
+      // data structure
+      std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
+
+      for(; entries>0; --entries) {
+        std::size_t row;
+        IndexData<D> data;
+        skipComments(file);
+        file>>row;
+        --row; // Index was 1 based.
+        assert(row/bcols<matrix.N());
+        file>>data;
+        assert(data.index/bcols<matrix.M());
+        rows[row].insert(data);
+      }
+
+      // TODO extend to capture the nongeneral cases.
+      if(mmHeader.structure!= general)
+        DUNE_THROW(Dune::NotImplemented, "Only general is supported right now!");
+
+      // Setup the matrix sparsity pattern
+      int nnz=0;
+      for(typename Matrix::CreateIterator iter=matrix.createbegin();
+          iter!= matrix.createend(); ++iter)
+      {
+        for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
+            brow<browend; ++brow)
+        {
+          typedef typename std::set<IndexData<D> >::const_iterator Siter;
+          for(Siter siter=rows[brow].begin(), send=rows[brow].end();
+              siter != send; ++siter, ++nnz)
+            iter.insert(siter->index/bcols);
+        }
+      }
+
+      //Set the matrix values
+      matrix=0;
+
+      MatrixValuesSetter<D,brows,bcols> Setter;
+
+      Setter(rows, matrix);
+    }
+  } // end namespace MatrixMarketImpl
+
+
+  template<typename T, typename A, int brows, int bcols>
+  void readMatrixMarket(Dune::BCRSMatrix<Ewoms::MatrixBlock<T,brows,bcols>,A>& matrix,
+                        std::istream& istr)
+  {
+    using namespace MatrixMarketImpl;
+
+    MMHeader header;
+    if(!readMatrixMarketBanner(istr, header)) {
+      std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
+                << "%%MatrixMarket matrix coordinate real general"<<std::endl;
+      // Go to the beginning of the file
+      istr.clear() ;
+      istr.seekg(0, std::ios::beg);
+    }
+    skipComments(istr);
+
+    std::size_t rows, cols, entries;
+
+    if(lineFeed(istr))
+      throw MatrixMarketFormatError();
+
+    istr >> rows;
+
+    if(lineFeed(istr))
+      throw MatrixMarketFormatError();
+    istr >> cols;
+
+    if(lineFeed(istr))
+      throw MatrixMarketFormatError();
+
+    istr >>entries;
+
+    std::size_t nnz, blockrows, blockcols;
+
+    std::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
+
+    istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
+
+
+    matrix.setSize(blockrows, blockcols);
+    matrix.setBuildMode(Dune::BCRSMatrix<Ewoms::MatrixBlock<T,brows,bcols>,A>::row_wise);
+
+    if(header.type==array_type)
+      DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
+
+    readSparseEntries(matrix, istr, entries, header, NumericWrapper<T>());
+  }
+
+
+
+  template<typename B, int i, int j, typename A>
+  struct mm_multipliers<BCRSMatrix<Ewoms::MatrixBlock<B,i,j>,A> >
+  {
+    enum {
+      rows = i,
+      cols = j
+    };
+  };
+
+  template<typename B, int i, int j>
+  void mm_print_entry(const Ewoms::MatrixBlock<B,i,j>& entry,
+                      typename Ewoms::MatrixBlock<B,i,j>::size_type rowidx,
+                      typename Ewoms::MatrixBlock<B,i,j>::size_type colidx,
+                      std::ostream& ostr)
+  {
+    typedef typename Ewoms::MatrixBlock<B,i,j>::const_iterator riterator;
+    typedef typename Ewoms::MatrixBlock<B,i,j>::ConstColIterator citerator;
+    for(riterator row=entry.begin(); row != entry.end(); ++row, ++rowidx) {
+      int coli=colidx;
+      for(citerator col = row->begin(); col != row->end(); ++col, ++coli)
+        ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
+    }
+  }
+
+  /** @} */
+}
+
+#endif // MATRIXMARKET_EWOMS_HH
diff --git a/ewoms/models/blackoil/blackoilenergymodules.hh b/ewoms/models/blackoil/blackoilenergymodules.hh
index 8dbd627e6c..0212145c09 100644
--- a/ewoms/models/blackoil/blackoilenergymodules.hh
+++ b/ewoms/models/blackoil/blackoilenergymodules.hh
@@ -352,13 +352,14 @@ public:
      */
     void updateTemperature_(const ElementContext& elemCtx,
                             unsigned dofIdx,
-                            unsigned timeIdx)
+                            unsigned timeIdx,
+                            unsigned focusTimeIdx)
     {
         auto& fs = asImp_().fluidState_;
         const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
 
         // set temperature
-        fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx));
+        fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx, focusTimeIdx));
     }
 
     /*!
@@ -420,7 +421,8 @@ class BlackOilEnergyIntensiveQuantities<TypeTag, false>
 public:
     void updateTemperature_(const ElementContext& elemCtx,
                             unsigned dofIdx,
-                            unsigned timeIdx)
+                            unsigned timeIdx,
+                            unsigned focusTimeIdx)
     {
         if (enableTemperature) {
             // even if energy is conserved, the temperature can vary over the spatial
diff --git a/ewoms/models/blackoil/blackoilintensivequantities.hh b/ewoms/models/blackoil/blackoilintensivequantities.hh
index e715715d4a..ec74bc42bd 100644
--- a/ewoms/models/blackoil/blackoilintensivequantities.hh
+++ b/ewoms/models/blackoil/blackoilintensivequantities.hh
@@ -110,14 +110,14 @@ public:
     /*!
      * \copydoc IntensiveQuantities::update
      */
-    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
+    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focusTimeIdx)
     {
-        ParentType::update(elemCtx, dofIdx, timeIdx);
+        ParentType::update(elemCtx, dofIdx, timeIdx, focusTimeIdx);
 
         const auto& problem = elemCtx.problem();
         const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
 
-        asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx);
+        asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx, focusTimeIdx);
 
         unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
         unsigned pvtRegionIdx = priVars.pvtRegionIndex();
@@ -126,21 +126,21 @@ public:
         // extract the water and the gas saturations for convenience
         Evaluation Sw = 0.0;
         if (waterEnabled)
-            Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx);
+            Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx,  focusTimeIdx);
 
         Evaluation Sg = 0.0;
         if (compositionSwitchEnabled)
         {
             if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg)
                 // -> threephase case
-                Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
+                Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx, focusTimeIdx);
             else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
                 // -> gas-water case
                 Sg = 1.0 - Sw;
 
                 // deal with solvent
                 if (enableSolvent)
-                    Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
+                    Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx, focusTimeIdx);
             }
             else
             {
@@ -157,7 +157,7 @@ public:
 
         // deal with solvent
         if (enableSolvent)
-            So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
+            So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx, focusTimeIdx);
 
         if (FluidSystem::phaseIsActive(waterPhaseIdx))
             fluidState_.setSaturation(waterPhaseIdx, Sw);
@@ -168,7 +168,7 @@ public:
         if (FluidSystem::phaseIsActive(oilPhaseIdx))
             fluidState_.setSaturation(oilPhaseIdx, So);
 
-        asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
+        asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx, focusTimeIdx);
 
         // now we compute all phase pressures
         Evaluation pC[numPhases];
@@ -177,14 +177,14 @@ public:
 
         //oil is the reference phase for pressure
         if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
-            const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
+            const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, focusTimeIdx);
             for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
                 if (FluidSystem::phaseIsActive(phaseIdx))
                     fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx]));
         }
 
         else {
-            const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
+            const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, focusTimeIdx);
             for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
                 if (FluidSystem::phaseIsActive(phaseIdx))
                     fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
@@ -196,7 +196,7 @@ public:
         Opm::Valgrind::CheckDefined(mobility_);
 
         // update the Saturation functions for the blackoil solvent module.
-        asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
+        asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx, focusTimeIdx);
 
         const Evaluation& SoMax =
             Opm::max(fluidState_.saturation(oilPhaseIdx),
@@ -236,7 +236,7 @@ public:
             Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(timeIdx, globalSpaceIdx);
 
             // oil phase, we can directly set the composition of the oil phase
-            const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
+            const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx, focusTimeIdx);
             fluidState_.setRs(Opm::min(RsMax, Rs));
 
             if (FluidSystem::enableVaporizedOil()) {
@@ -257,7 +257,7 @@ public:
         else {
             assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv);
 
-            const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
+            const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx, focusTimeIdx);
             fluidState_.setRv(Rv);
 
             if (FluidSystem::enableDissolvedGas()) {
@@ -343,12 +343,13 @@ public:
             // deal with water induced rock compaction
             porosity_ *= problem.template rockCompPoroMultiplier<Evaluation>(*this, globalSpaceIdx);
 
-        asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
-        asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
+        asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx, focusTimeIdx);
+        asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx, focusTimeIdx);
         asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache);
 
         // update the quantities which are required by the chosen
         // velocity model
+        // TODO: use the focus time here
         FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
 
 #ifndef NDEBUG
diff --git a/ewoms/models/blackoil/blackoilmodel.hh b/ewoms/models/blackoil/blackoilmodel.hh
index 6c77b4ff5e..ad70d70c17 100644
--- a/ewoms/models/blackoil/blackoilmodel.hh
+++ b/ewoms/models/blackoil/blackoilmodel.hh
@@ -512,6 +512,7 @@ public:
             }
         }
 
+        // this is restart spesific 
         this->solution(/*timeIdx=*/1) = this->solution(/*timeIdx=*/0);
     }
 
diff --git a/ewoms/models/blackoil/blackoilpolymermodules.hh b/ewoms/models/blackoil/blackoilpolymermodules.hh
index 8e08561269..438b9f358f 100644
--- a/ewoms/models/blackoil/blackoilpolymermodules.hh
+++ b/ewoms/models/blackoil/blackoilpolymermodules.hh
@@ -887,6 +887,7 @@ public:
         // Use log(v0) as initial value for u
         auto u = v0AbsLog;
         bool converged = false;
+        // TODO make this into paramters
         for (int i = 0; i < 20; ++i) {
             auto f = F(u);
             auto df = dF(u);
@@ -1048,12 +1049,13 @@ public:
      */
     void polymerPropertiesUpdate_(const ElementContext& elemCtx,
                                   unsigned dofIdx,
-                                  unsigned timeIdx)
+                                  unsigned timeIdx,
+                                  unsigned focusTimeIdx)
     {
         const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
-        polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx, timeIdx);
+        polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx, timeIdx, focusTimeIdx);
         if (enablePolymerMolarWeight) {
-            polymerMoleWeight_ = priVars.makeEvaluation(polymerMoleWeightIdx, timeIdx);
+            polymerMoleWeight_ = priVars.makeEvaluation(polymerMoleWeightIdx, timeIdx, focusTimeIdx);
         }
         const Scalar cmax = PolymerModule::plymaxMaxConcentration(elemCtx, dofIdx, timeIdx);
 
@@ -1167,7 +1169,9 @@ class BlackOilPolymerIntensiveQuantities<TypeTag, false>
 public:
     void polymerPropertiesUpdate_(const ElementContext& elemCtx OPM_UNUSED,
                                   unsigned scvIdx OPM_UNUSED,
-                                  unsigned timeIdx OPM_UNUSED)
+                                  unsigned timeIdx OPM_UNUSED,
+                                  unsigned focusTimeIdxIdx OPM_UNUSED)
+
     { }
 
     const Evaluation& polymerMoleWeight() const
diff --git a/ewoms/models/blackoil/blackoilsolventmodules.hh b/ewoms/models/blackoil/blackoilsolventmodules.hh
index b5a6dfea13..b6d89a4017 100644
--- a/ewoms/models/blackoil/blackoilsolventmodules.hh
+++ b/ewoms/models/blackoil/blackoilsolventmodules.hh
@@ -903,11 +903,12 @@ public:
      */
     void solventPreSatFuncUpdate_(const ElementContext& elemCtx,
                                   unsigned dofIdx,
-                                  unsigned timeIdx)
+                                  unsigned timeIdx,
+                                  unsigned focustimeIdx)
     {
         const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
         auto& fs = asImp_().fluidState_;
-        solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx);
+        solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx,focustimeIdx);
         hydrocarbonSaturation_ = fs.saturation(gasPhaseIdx);
 
         // apply a cut-off. Don't waste calculations if no solvent
@@ -928,7 +929,8 @@ public:
      */
     void solventPostSatFuncUpdate_(const ElementContext& elemCtx,
                                    unsigned dofIdx,
-                                   unsigned timeIdx)
+                                   unsigned timeIdx,
+                                   unsigned focustimeidx)
     {
         // revert the gas "saturation" of the fluid state back to the saturation of the
         // hydrocarbon gas.
@@ -957,9 +959,9 @@ public:
 
             //oil is the reference phase for pressure
             if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv)
-                pgMisc = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
+                pgMisc = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, focustimeidx);
             else {
-                const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
+                const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, focustimeidx);
                 pgMisc = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
             }
 
@@ -1043,7 +1045,8 @@ public:
      */
     void solventPvtUpdate_(const ElementContext& elemCtx,
                            unsigned scvIdx,
-                           unsigned timeIdx)
+                           unsigned timeIdx,
+                           unsigned focustimeidx)
     {
         const auto& iq = asImp_();
         const auto& fs = iq.fluidState();
@@ -1278,17 +1281,20 @@ class BlackOilSolventIntensiveQuantities<TypeTag, false>
 public:
     void solventPreSatFuncUpdate_(const ElementContext& elemCtx OPM_UNUSED,
                                   unsigned scvIdx OPM_UNUSED,
-                                  unsigned timeIdx OPM_UNUSED)
+                                  unsigned timeIdx OPM_UNUSED,
+                                  unsigned focustimeidx OPM_UNUSED)
     { }
 
     void solventPostSatFuncUpdate_(const ElementContext& elemCtx OPM_UNUSED,
                                    unsigned scvIdx OPM_UNUSED,
-                                   unsigned timeIdx OPM_UNUSED)
+                                   unsigned timeIdx OPM_UNUSED,
+                                   unsigned focustimeidx OPM_UNUSED )
     { }
 
     void solventPvtUpdate_(const ElementContext& elemCtx OPM_UNUSED,
                            unsigned scvIdx OPM_UNUSED,
-                           unsigned timeIdx OPM_UNUSED)
+                           unsigned timeIdx OPM_UNUSED,
+                           unsigned focustimeIdx OPM_UNUSED)
     { }
 
     const Evaluation& solventSaturation() const
diff --git a/ewoms/models/blackoil/blackoiltwophaseindices.hh b/ewoms/models/blackoil/blackoiltwophaseindices.hh
index c7f21835ba..81ac6ba62e 100644
--- a/ewoms/models/blackoil/blackoiltwophaseindices.hh
+++ b/ewoms/models/blackoil/blackoiltwophaseindices.hh
@@ -78,6 +78,7 @@ struct BlackOilTwoPhaseIndices
 
     //! Index of the oil pressure in a vector of primary variables
     static const int pressureSwitchIdx  = waterEnabled ? PVOffset + 1 : PVOffset + 0;
+    //static const int pressureSwitchIdx  = PVOffset + 0;
 
     /*!
      * \brief Index of the switching variable which determines the composition of the
diff --git a/ewoms/models/flash/flashintensivequantities.hh b/ewoms/models/flash/flashintensivequantities.hh
index cd9e632010..a2cc8a4d7d 100644
--- a/ewoms/models/flash/flashintensivequantities.hh
+++ b/ewoms/models/flash/flashintensivequantities.hh
@@ -99,9 +99,9 @@ public:
     /*!
      * \copydoc IntensiveQuantities::update
      */
-    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
+    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focusTimeIdx)
     {
-        ParentType::update(elemCtx, dofIdx, timeIdx);
+        ParentType::update(elemCtx, dofIdx, timeIdx, focusTimeIdx);
         EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);
 
         const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
diff --git a/ewoms/models/immiscible/immiscibleintensivequantities.hh b/ewoms/models/immiscible/immiscibleintensivequantities.hh
index cad51acf22..468f8fef4e 100644
--- a/ewoms/models/immiscible/immiscibleintensivequantities.hh
+++ b/ewoms/models/immiscible/immiscibleintensivequantities.hh
@@ -90,9 +90,9 @@ public:
     /*!
      * \copydoc IntensiveQuantities::update
      */
-    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
+    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focusTimeIdx)
     {
-        ParentType::update(elemCtx, dofIdx, timeIdx);
+        ParentType::update(elemCtx, dofIdx, timeIdx, focusTimeIdx);
         EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);
 
         // material law parameters
diff --git a/ewoms/models/ncp/ncpintensivequantities.hh b/ewoms/models/ncp/ncpintensivequantities.hh
index c125e9a43e..dc6ede7fc7 100644
--- a/ewoms/models/ncp/ncpintensivequantities.hh
+++ b/ewoms/models/ncp/ncpintensivequantities.hh
@@ -101,9 +101,10 @@ public:
      */
     void update(const ElementContext& elemCtx,
                 unsigned dofIdx,
-                unsigned timeIdx)
+                unsigned timeIdx,
+                unsigned focusTimeIdx)
     {
-        ParentType::update(elemCtx, dofIdx, timeIdx);
+        ParentType::update(elemCtx, dofIdx, timeIdx, focusTimeIdx);
         ParentType::checkDefined();
 
         typename FluidSystem::template ParameterCache<Evaluation> paramCache;
diff --git a/ewoms/models/pvs/pvsintensivequantities.hh b/ewoms/models/pvs/pvsintensivequantities.hh
index 3c64276a58..b62fd5959f 100644
--- a/ewoms/models/pvs/pvsintensivequantities.hh
+++ b/ewoms/models/pvs/pvsintensivequantities.hh
@@ -105,9 +105,9 @@ public:
     /*!
      * \copydoc ImmiscibleIntensiveQuantities::update
      */
-    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
+    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focusTimeIdx)
     {
-        ParentType::update(elemCtx, dofIdx, timeIdx);
+        ParentType::update(elemCtx, dofIdx, timeIdx, focusTimeIdx);
         EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);
 
         const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
diff --git a/ewoms/models/richards/richardsintensivequantities.hh b/ewoms/models/richards/richardsintensivequantities.hh
index 4c696c264a..1cb146e734 100644
--- a/ewoms/models/richards/richardsintensivequantities.hh
+++ b/ewoms/models/richards/richardsintensivequantities.hh
@@ -84,9 +84,9 @@ public:
     /*!
      * \copydoc IntensiveQuantities::update
      */
-    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
+    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focusTimeIdx)
     {
-        ParentType::update(elemCtx, dofIdx, timeIdx);
+        ParentType::update(elemCtx, dofIdx, timeIdx, focusTimeIdx);
 
         const auto& T = elemCtx.problem().temperature(elemCtx, dofIdx, timeIdx);
         fluidState_.setTemperature(T);
diff --git a/ewoms/nonlinear/newtonmethod.hh b/ewoms/nonlinear/newtonmethod.hh
index f86dcfc309..dbd621d664 100644
--- a/ewoms/nonlinear/newtonmethod.hh
+++ b/ewoms/nonlinear/newtonmethod.hh
@@ -365,10 +365,12 @@ public:
                               << std::flush;
                 }
 
+                //
+                int focusTimeIdx = 0;
                 // do the actual linearization
                 linearizeTimer_.start();
-                asImp_().linearizeDomain_();
-                asImp_().linearizeAuxiliaryEquations_();
+                asImp_().linearizeDomain_(focusTimeIdx);
+                asImp_().linearizeAuxiliaryEquations_(focusTimeIdx);
                 linearizeTimer_.stop();
 
                 solveTimer_.start();
@@ -647,14 +649,14 @@ protected:
      * \brief Linearize the global non-linear system of equations associated with the
      *        spatial domain.
      */
-    void linearizeDomain_()
+    void linearizeDomain_(int focusTimeIdx)
     {
-        model().linearizer().linearizeDomain();
+        model().linearizer().linearizeDomain(focusTimeIdx);
     }
 
-    void linearizeAuxiliaryEquations_()
+    void linearizeAuxiliaryEquations_(int focusTimeIdx)
     {
-        model().linearizer().linearizeAuxiliaryEquations();
+        model().linearizer().linearizeAuxiliaryEquations(focusTimeIdx);
         model().linearizer().finalize();
     }