diff --git a/dumux/porousmediumflow/2p1c/CMakeLists.txt b/dumux/porousmediumflow/2p1c/CMakeLists.txt
index ba8341c614f1a2c797c95f5402f602025f1087b1..0373b7ebc4e621671c07ab48e6cf5ad0a41801f4 100644
--- a/dumux/porousmediumflow/2p1c/CMakeLists.txt
+++ b/dumux/porousmediumflow/2p1c/CMakeLists.txt
@@ -1 +1,11 @@
-add_subdirectory("implicit")
+
+#install headers
+install(FILES
+darcyslaw.hh
+indices.hh
+localresidual.hh
+model.hh
+primaryvariableswitch.hh
+volumevariables.hh
+vtkoutputfields.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/2p1c/implicit)
diff --git a/dumux/porousmediumflow/2p1c/darcyslaw.hh b/dumux/porousmediumflow/2p1c/darcyslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..727dc8b9785da8219d88ac88e685672245c05ccb
--- /dev/null
+++ b/dumux/porousmediumflow/2p1c/darcyslaw.hh
@@ -0,0 +1,156 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program 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 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief This file contains the data which is required to calculate
+ *        volume and mass fluxes of fluid phases over a face of a finite volume by means
+ *        of the Darcy approximation. Specializations are provided for the different discretization methods.
+ */
+#ifndef DUMUX_2P1C_SPURIOUS_FLUX_BLOCKING_DARCYS_LAW_HH
+#define DUMUX_2P1C_SPURIOUS_FLUX_BLOCKING_DARCYS_LAW_HH
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/darcyslaw.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+     NEW_PROP_TAG(UseBlockingOfSpuriousFlow); //!< Determines whether Blocking ofspurious flow is used
+}
+
+/*!
+ * \ingroup DarcysLaw
+ * \brief Specialization of Darcy's Law for the box method.
+ */
+template <class TypeTag>
+class TwoPOneCDarcysLaw : public DarcysLaw<TypeTag>
+{
+    using ParentType = DarcysLaw<TypeTag>;
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using ElemFluxVarCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+    using FluxVarCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using IndexType = typename GridView::IndexSet::IndexType;
+    using CoordScalar = typename GridView::ctype;
+
+    enum { dim = GridView::dimension};
+    enum { dimWorld = GridView::dimensionworld};
+
+    // copy some indices for convenience
+    enum {
+        // phase indices
+        wPhaseIdx = Indices::wPhaseIdx,
+        nPhaseIdx = Indices::nPhaseIdx,
+    };
+
+    using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+public:
+
+    static Scalar flux(const Problem& problem,
+                       const Element& element,
+                       const FVElementGeometry& fvGeometry,
+                       const ElementVolumeVariables& elemVolVars,
+                       const SubControlVolumeFace& scvf,
+                       const IndexType phaseIdx,
+                       const ElemFluxVarCache& elemFluxVarCache)
+    {
+        const Scalar flux = ParentType::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarCache);
+        if((!GET_PROP_VALUE(TypeTag, UseBlockingOfSpuriousFlow)) || phaseIdx != wPhaseIdx)
+            return flux;
+
+        const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+        const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+
+        const Scalar factor = std::signbit(flux) ? factor_(outsideVolVars, insideVolVars, phaseIdx) :
+                                                   factor_(insideVolVars, outsideVolVars, phaseIdx);
+
+        return flux * factor;
+    }
+
+private:
+    /*!
+     * \brief Calculate the blocking factor which prevents spurious cold water fluxes into the steam zone (Gudbjerg et al., 2005) \cite gudbjerg2004 <BR>
+     *
+     * \param up The upstream volume variables
+     * \param dn The downstream volume variables
+     * \param phaseIdx The index of the fluid phase
+     */
+    static Scalar factor_(const VolumeVariables &up, const VolumeVariables &dn,const  int phaseIdx)
+    {
+
+
+        Scalar factor = 1.0;
+
+        Scalar tDn = dn.temperature(); //temperature of the downstream SCV (where the cold water is potentially intruding into a steam zone)
+        Scalar tUp = up.temperature(); //temperature of the upstream SCV
+
+        Scalar sgDn = dn.saturation(nPhaseIdx); //gas phase saturation of the downstream SCV
+        Scalar sgUp = up.saturation(nPhaseIdx); //gas phase saturation of the upstream SCV
+
+        bool upIsNotSteam = false;
+        bool downIsSteam = false;
+        bool spuriousFlow = false;
+
+        if(sgUp <= 1e-5)
+        {upIsNotSteam = true;}
+
+        if(sgDn > 1e-5)
+        {downIsSteam = true;}
+
+        if(upIsNotSteam && downIsSteam  && tDn > tUp && phaseIdx == wPhaseIdx)
+        {
+          spuriousFlow = true;
+        //   spuriousFlowDetected_ = true;
+        }
+
+        if(spuriousFlow)
+        {
+          Scalar deltaT = tDn - tUp;
+
+          if((deltaT) > 15 )
+          {factor = 0.0 ;}
+
+          else
+          { factor = 1-(deltaT/15); }
+
+        }
+        return factor;
+    }
+
+
+};
+
+} // end namespace Dumux
+
+#endif // DUMUX_2P1C_SPURIOUS_FLUX_BLOCKING_DARCYS_LAW_HH
diff --git a/dumux/porousmediumflow/2p1c/implicit/CMakeLists.txt b/dumux/porousmediumflow/2p1c/implicit/CMakeLists.txt
deleted file mode 100644
index d5c0899adf06164a0f194e86e8a06f13d6b90c91..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2p1c/implicit/CMakeLists.txt
+++ /dev/null
@@ -1,11 +0,0 @@
-
-#install headers
-install(FILES
-indices.hh
-localresidual.hh
-model.hh
-newtoncontroller.hh
-properties.hh
-propertydefaults.hh
-volumevariables.hh
-DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/2p1c/implicit)
diff --git a/dumux/porousmediumflow/2p1c/implicit/localresidual.hh b/dumux/porousmediumflow/2p1c/implicit/localresidual.hh
deleted file mode 100644
index 51a772df81669f7b9b2821a65cf12b2256e535fe..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2p1c/implicit/localresidual.hh
+++ /dev/null
@@ -1,272 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program 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 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Element-wise calculation of the Jacobian matrix for problems
- *        using the two-phase one-component fully implicit model.
- *
- * Important note: The 2p1c model requires the use of the non-isothermal extension found in dumux/porousmediumflow/nonisothermal/implicit/
- */
-#ifndef DUMUX_2P1C_LOCAL_RESIDUAL_HH
-#define DUMUX_2P1C_LOCAL_RESIDUAL_HH
-
-#include "properties.hh"
-
-namespace Dumux
-{
-/*!
- * \ingroup TwoPOneCModel
- * \ingroup ImplicitLocalResidual
- * \brief Element-wise calculation of the Jacobian matrix for problems
- *        using the two-phase one-component fully implicit model.
- *
- * This class depends on the non-isothermal model.
- *
- */
-template<class TypeTag>
-class TwoPOneCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
-{
-protected:
-    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-    typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
-
-
-    enum {
-        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
-        numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
-
-        conti0EqIdx = Indices::conti0EqIdx, //!< Index of the mass conservation equation for the water component
-        energyEqIdx = Indices::energyEqIdx, //!< Index of the energy conservation equation
-        wPhaseIdx = Indices::wPhaseIdx,
-        gPhaseIdx = Indices::gPhaseIdx,
-    };
-
-    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
-
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    typedef typename GridView::template Codim<0>::Entity Element;
-
-    static const bool useBlockingOfSpuriousFlow = GET_PROP_VALUE(TypeTag, UseBlockingOfSpuriousFlow);
-
-public:
-
-    /*!
-     * \brief Evaluate the amount all conservation quantities
-     *        (e.g. phase mass) within a sub-control volume.
-     *
-     * The result should be averaged over the volume (e.g. phase mass
-     * inside a sub control volume divided by the volume)
-     *
-     *  \param storage The mass of the component within the sub-control volume
-     *  \param scvIdx The SCV (sub-control-volume) index
-     *  \param usePrevSol Evaluate function with solution of current or previous time step
-     */
-    void computeStorage(PrimaryVariables &storage, const int scvIdx, bool usePrevSol) const
-    {
-        // if flag usePrevSol is set, the solution from the previous
-        // time step is used, otherwise the current solution is
-        // used. The secondary variables are used accordingly.  This
-        // is required to compute the derivative of the storage term
-        // using the implicit euler method.
-        const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
-        const VolumeVariables &volVars = elemVolVars[scvIdx];
-
-        // compute storage term of all components within all phases
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-        {
-            storage[conti0EqIdx] +=
-                volVars.porosity()
-                * volVars.saturation(phaseIdx) * volVars.density(phaseIdx);
-        }
-    }
-
-    /*!
-     * \brief Evaluates the total flux of all conservation quantities
-     *        over a face of a sub-control volume.
-     *
-     * \param flux The flux over the SCV (sub-control-volume) face for each component
-     * \param faceIdx The index of the SCV face
-     * \param onBoundary A boolean variable to specify whether the flux variables
-     *        are calculated for interior SCV faces or boundary faces, default=false
-     */
-    void computeFlux(PrimaryVariables &flux, const int faceIdx, const bool onBoundary=false) const
-    {
-        FluxVariables fluxVars;
-        fluxVars.update(this->problem_(),
-                           this->element_(),
-                           this->fvGeometry_(),
-                           faceIdx,
-                           this->curVolVars_(),
-                           onBoundary);
-
-        flux = 0;
-        computeAdvectiveFlux(flux, fluxVars);
-        asImp_()->computeDiffusiveFlux(flux, fluxVars);
-    }
-
-    /*!
-     * \brief Evaluates the advective mass flux of all components over
-     *        a face of a subcontrol volume.
-     *
-     * \param flux The advective flux over the sub-control-volume face for each component
-     * \param fluxVars The flux variables at the current SCV
-     */
-
-    void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
-    {
-        Scalar massUpwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
-
-        // loop over all phases
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-        {
-            // data attached to upstream and the downstream vertices
-            // of the current phase
-            const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
-            const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
-
-            Scalar factor = 1.0;
-
-            // factor for the blocking of spurious cold water fluxes into the steam zone (Gudbjerg, 2004)
-            if(useBlockingOfSpuriousFlow)
-            { factor = factor_(up, dn, phaseIdx); }
-
-
-            // advective mass flux
-            flux[conti0EqIdx] +=  fluxVars.volumeFlux(phaseIdx)
-                                  * (massUpwindWeight
-                                  * up.fluidState().density(phaseIdx)
-                                  + (1.0 - massUpwindWeight)
-                                  * dn.fluidState().density(phaseIdx))
-                                  * factor ;
-
-            // advective heat flux
-            flux[energyEqIdx] += fluxVars.volumeFlux(phaseIdx)
-                                 *  (massUpwindWeight
-                                 *  up.density(phaseIdx)
-                                 *  up.enthalpy(phaseIdx)
-                                 +  (1.0 - massUpwindWeight)
-                                 *  dn.density(phaseIdx)
-                                 *  dn.enthalpy(phaseIdx))
-                                 *  factor ;
-        }
-    }
-
-    /*!
-     * \brief Adds the diffusive mass flux of all components over
-     *        a face of a subcontrol volume.
-     *
-     * \param flux The diffusive flux over the sub-control-volume face for each component
-     * \param fluxVars The flux variables at the current SCV
-     */
-    void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
-    {}
-
-    /*!
-     * \brief Returns true if a spurious flow has been detected
-     *
-     */
-    bool spuriousFlowDetected() const
-    {
-        return spuriousFlowDetected_;
-    }
-
-    /*!
-     * \brief Used to reset the respective flag
-     *
-     */
-    void resetSpuriousFlowDetected()
-    {
-        spuriousFlowDetected_ = false;
-    }
-
-protected:
-
-    /*!
-     * \brief Calculate the blocking factor which prevents spurious cold water fluxes into the steam zone (Gudbjerg et al., 2005) \cite gudbjerg2004 <BR>
-     *
-     * \param up The upstream volume variables
-     * \param dn The downstream volume variables
-     * \param phaseIdx The index of the fluid phase
-     */
-    Scalar factor_(const VolumeVariables &up, const VolumeVariables &dn,const  int phaseIdx) const {
-
-        Scalar factor = 1.0;
-
-        Scalar tDn = dn.temperature(); //temperature of the downstream SCV (where the cold water is potentially intruding into a steam zone)
-        Scalar tUp = up.temperature(); //temperature of the upstream SCV
-
-        Scalar sgDn = dn.saturation(gPhaseIdx); //gas phase saturation of the downstream SCV
-        Scalar sgUp = up.saturation(gPhaseIdx); //gas phase saturation of the upstream SCV
-
-        bool upIsNotSteam = false;
-        bool downIsSteam = false;
-        bool spuriousFlow = false;
-
-        if(sgUp <= 1e-5)
-        {upIsNotSteam = true;}
-
-        if(sgDn > 1e-5)
-        {downIsSteam = true;}
-
-        if(upIsNotSteam && downIsSteam  && tDn > tUp && phaseIdx == wPhaseIdx)
-        {
-          spuriousFlow = true;
-          spuriousFlowDetected_ = true;
-        }
-
-        if(spuriousFlow)
-        {
-          Scalar deltaT = tDn - tUp;
-
-          if((deltaT) > 15 )
-          {factor = 0.0 ;}
-
-          else
-          { factor = 1-(deltaT/15); }
-
-        }
-        return factor;
-    }
-
-    Implementation *asImp_()
-    {
-        return static_cast<Implementation *> (this);
-    }
-
-    const Implementation *asImp_() const
-    {
-        return static_cast<const Implementation *> (this);
-    }
-
-private:
-    mutable bool spuriousFlowDetected_;// true if spurious cold-water flow into a steam zone is happening
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/porousmediumflow/2p1c/implicit/model.hh b/dumux/porousmediumflow/2p1c/implicit/model.hh
deleted file mode 100644
index 0ec78078c25bc6ad6d20375d57715ea7fcff2146..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2p1c/implicit/model.hh
+++ /dev/null
@@ -1,624 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program 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 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Adaption of the fully implicit scheme to the two-phase one-component flow model.
- *
- */
-#ifndef DUMUX_2P1C_MODEL_HH
-#define DUMUX_2P1C_MODEL_HH
-
-#include <dumux/porousmediumflow/implicit/velocityoutput.hh>
-#include "properties.hh"
-
-namespace Dumux
-{
-/*!
- * \ingroup TwoPOneCModel
- * \brief Adaption of the fully implicit scheme to the two-phase one-component flow model.
- *
- * \note The 2p1c model requires the use of the non-isothermal extension found in dumux/implicit/nonisothermal.
- *
- * This model is designed for simulating two fluid phases with water as the only component.
- * It is particularly suitable for the simulation of steam injection in saturated conditions.
- *
- * The model implements the flow of two phases and one component, i.e. a pure liquid (e.g. water)
- * and its vapor (e.g. steam),
- * \f$\alpha \in \{ w, n \}\f$ using a standard multiphase Darcy
- * approach as the equation for the conservation of momentum, i.e.
- \f[
- v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K}
- \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right)
- \f]
- *
- * By inserting this into the equation for the conservation of the
- * phase mass, one gets
- \f[
-\phi \frac{\partial\ \sum_\alpha (\rho_\alpha S_\alpha)}{\partial t} \\-\sum \limits_ \alpha \text{div} \left \{\rho_\alpha \frac{k_{r\alpha}}{\mu_\alpha}
-\mathbf{K} (\mathbf{grad}p_\alpha - \rho_\alpha \mathbf{g}) \right \} -q^w =0
- \f]
- *
- * All equations are discretized using a vertex-centered finite volume (box)
- * or cell-centered finite volume scheme as spatial
- * and the implicit Euler method as time discretization.
- *
- * By using constitutive relations for the capillary pressure \f$p_c =
- * p_n - p_w\f$ and relative permeability \f$k_{r\alpha}\f$ and taking
- * advantage of the fact that \f$S_w + S_n = 1\f$, the number of
- * unknowns can be reduced to two. The model features a primary variable switch.
- * If only one phase is present, \f$p_g\f$ and \f$T\f$ are the primary variables.
- * In the presence of two phases, \f$p_g\f$ and \f$S_w\f$ become the new primary variables.
- */
-template<class TypeTag>
-class TwoPOneCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
-{
-    typedef typename GET_PROP_TYPE(TypeTag, BaseModel) ParentType;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-
-    enum {
-        dim = GridView::dimension,
-        dimWorld = GridView::dimensionworld,
-
-        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
-        numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
-
-        pressureIdx = Indices::pressureIdx,
-        switch1Idx = Indices::switch1Idx,
-
-        wPhaseIdx = Indices::wPhaseIdx,
-        gPhaseIdx = Indices::gPhaseIdx,
-
-        twoPhases = Indices::twoPhases,
-        wPhaseOnly  = Indices::wPhaseOnly,
-        gPhaseOnly  = Indices::gPhaseOnly,
-    };
-
-    typedef typename GridView::template Codim<dim>::Entity Vertex;
-    typedef typename GridView::template Codim<0>::Entity Element;
-    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
-    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
-    typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix;
-
-    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
-
-
-    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
-    enum { dofCodim = isBox ? dim : 0 };
-
-public:
-    /*!
-     * \brief Initialize the static data with the initial solution.
-     *
-     * \param problem The problem to be solved
-     */
-    void init(Problem &problem)
-    {
-        ParentType::init(problem);
-
-        staticDat_.resize(this->numDofs());
-
-        setSwitched_(false);
-
-        if (isBox)
-        {
-            for(const auto& vertex: vertices(this->gridView_()))
-            {
-                int globalIdx = this->dofMapper().index(vertex);
-                const GlobalPosition &globalPos = vertex.geometry().corner(0);
-
-                // initialize phase presence
-                staticDat_[globalIdx].phasePresence
-                    = this->problem_().initialPhasePresence(vertex, globalIdx,
-                                                        globalPos);
-                staticDat_[globalIdx].wasSwitched = false;
-
-                staticDat_[globalIdx].oldPhasePresence
-                    = staticDat_[globalIdx].phasePresence;
-            }
-        }
-        else
-        {
-            for (const auto& element : elements(this->gridView_()))
-            {
-                int globalIdx = this->dofMapper().index(element);
-                const GlobalPosition &globalPos = element.geometry().center();
-
-                // initialize phase presence
-                staticDat_[globalIdx].phasePresence
-                    = this->problem_().initialPhasePresence(*this->gridView_().template begin<dim> (),
-                                                            globalIdx, globalPos);
-                staticDat_[globalIdx].wasSwitched = false;
-
-                staticDat_[globalIdx].oldPhasePresence
-                    = staticDat_[globalIdx].phasePresence;
-            }
-        }
-    }
-
-    /*!
-     * \brief Called by the update() method if applying the newton
-     *         method was unsuccessful.
-     */
-    void updateFailed()
-    {
-        ParentType::updateFailed();
-
-        setSwitched_(false);
-        resetPhasePresence_();
-    };
-
-    /*!
-     * \brief Called by the problem if a time integration was
-     *        successful, post processing of the solution is done and the
-     *        result has been written to disk.
-     *
-     * This should prepare the model for the next time integration.
-     */
-    void advanceTimeLevel()
-    {
-        ParentType::advanceTimeLevel();
-
-        // update the phase state
-        updateOldPhasePresence_();
-        setSwitched_(false);
-    }
-
-    /*!
-     * \brief Return true if the primary variables were switched for
-     *        at least one vertex after the last timestep.
-     */
-    bool switched() const
-    {
-        return switchFlag_;
-    }
-
-    /*!
-     * \brief Returns the phase presence of the current or the old solution of a degree of freedom.
-     *
-     * \param globalIdx The global index of the degree of freedom
-     * \param oldSol Evaluate function with solution of current or previous time step
-     */
-    int phasePresence(int globalIdx, bool oldSol) const
-    {
-        return
-            oldSol
-            ? staticDat_[globalIdx].oldPhasePresence
-            : staticDat_[globalIdx].phasePresence;
-    }
-
-    /*!
-     * \brief Append all quantities of interest which can be derived
-     *        from the solution of the current time step to the VTK
-     *        writer.
-     *
-     * \param sol The solution vector
-     * \param writer The writer for multi-file VTK datasets
-     */
-    template<class MultiWriter>
-    void addOutputVtkFields(const SolutionVector &sol,
-                            MultiWriter &writer)
-    {
-        typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
-        typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorField;
-
-        // get the number of degrees of freedom
-        unsigned numDofs = this->numDofs();
-
-        // create the required scalar fields
-        ScalarField *saturation[numPhases];
-        ScalarField *pressure[numPhases];
-        ScalarField *density[numPhases];
-
-        ScalarField *mobility[numPhases];
-        ScalarField *viscosity[numPhases];
-
-
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-            saturation[phaseIdx] = writer.allocateManagedBuffer(numDofs);
-            pressure[phaseIdx] = writer.allocateManagedBuffer(numDofs);
-            density[phaseIdx] = writer.allocateManagedBuffer(numDofs);
-
-            mobility[phaseIdx] = writer.allocateManagedBuffer(numDofs);
-            viscosity[phaseIdx] = writer.allocateManagedBuffer(numDofs);
-        }
-
-        ScalarField *phasePresence = writer.allocateManagedBuffer (numDofs);
-        ScalarField *poro = writer.allocateManagedBuffer(numDofs);
-        ScalarField *permXX = writer.allocateManagedBuffer(numDofs);
-        ScalarField *permYY = writer.allocateManagedBuffer(numDofs);
-        ScalarField *permZZ = writer.allocateManagedBuffer(numDofs);
-        VectorField *velocityW = writer.template allocateManagedBuffer<double, dim>(numDofs);
-        VectorField *velocityG = writer.template allocateManagedBuffer<double, dim>(numDofs);
-        ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
-
-        if (velocityOutput.enableOutput()) // check if velocity output is demanded
-        {
-            // initialize velocity fields
-            for (unsigned int i = 0; i < numDofs; ++i)
-            {
-                (*velocityW)[i] = Scalar(0);
-                (*velocityG)[i] = Scalar(0);
-            }
-        }
-
-        unsigned numElements = this->gridView_().size(0);
-        ScalarField *rank = writer.allocateManagedBuffer (numElements);
-
-        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
-        {
-            int eIdx = this->elementMapper().index(element);
-            (*rank)[eIdx] = this->gridView_().comm().rank();
-
-            FVElementGeometry fvGeometry;
-            fvGeometry.update(this->gridView_(), element);
-
-
-            ElementVolumeVariables elemVolVars;
-            elemVolVars.update(this->problem_(),
-                            element,
-                            fvGeometry,
-                            false /* oldSol? */);
-
-            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-            {
-                int globalIdx = this->dofMapper().subIndex(element, scvIdx, dofCodim);
-
-                for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-                    (*saturation[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().saturation(phaseIdx);
-                    (*pressure[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().pressure(phaseIdx);
-                    (*density[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().density(phaseIdx);
-
-                    (*mobility[phaseIdx])[globalIdx] = elemVolVars[scvIdx].mobility(phaseIdx);
-                    (*viscosity[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().viscosity(phaseIdx);
-                }
-
-                (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
-                (*phasePresence)[globalIdx] = staticDat_[globalIdx].phasePresence;
-
-                FieldMatrix K = this->problem_().spatialParams().intrinsicPermeability(element,
-                                                                                    fvGeometry,
-                                                                                    scvIdx);
-                (*permXX)[globalIdx] = K[0][0];
-                if(dimWorld > 1)
-                    (*permYY)[globalIdx] = K[1][1];
-                if(dimWorld > 2)
-                    (*permZZ)[globalIdx] = K[2][2];
-            }
-
-            // velocity output
-            velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
-            velocityOutput.calculateVelocity(*velocityG, elemVolVars, fvGeometry, element, gPhaseIdx);
-        }
-
-        writer.attachDofData(*saturation[wPhaseIdx], "sw", isBox);
-        writer.attachDofData(*saturation[gPhaseIdx], "sg", isBox);
-        writer.attachDofData(*pressure[wPhaseIdx], "pw", isBox);
-        writer.attachDofData(*pressure[gPhaseIdx], "pg", isBox);
-        writer.attachDofData(*density[wPhaseIdx], "rhow", isBox);
-        writer.attachDofData(*density[gPhaseIdx], "rhog", isBox);
-
-        writer.attachDofData(*mobility[wPhaseIdx], "MobW", isBox);
-        writer.attachDofData(*mobility[gPhaseIdx], "MobG", isBox);
-
-        writer.attachDofData(*viscosity[wPhaseIdx], "ViscosW", isBox);
-        writer.attachDofData(*viscosity[gPhaseIdx], "ViscosG", isBox);
-
-        writer.attachDofData(*poro, "porosity", isBox);
-        writer.attachDofData(*permXX, "permeabilityXX", isBox);
-        if(dimWorld > 1)
-            writer.attachDofData(*permYY, "permeabilityYY", isBox);
-        if(dimWorld > 2)
-            writer.attachDofData(*permZZ, "permeabilityZZ", isBox);
-        writer.attachDofData(*phasePresence, "phase presence", isBox);
-
-        if (velocityOutput.enableOutput()) // check if velocity output is demanded
-        {
-            writer.attachDofData(*velocityW,  "velocityW", isBox, dim);
-            writer.attachDofData(*velocityG,  "velocityG", isBox, dim);
-        }
-
-        writer.attachCellData(*rank, "process rank");
-    }
-
-    /*!
-     * \brief Write the current solution to a restart file.
-     *
-     * \param outStream The output stream of one vertex for the restart file
-     * \param entity The entity, either a vertex or an element
-     */
-    template<class Entity>
-    void serializeEntity(std::ostream &outStream, const Entity &entity)
-    {
-        // write primary variables
-        ParentType::serializeEntity(outStream, entity);
-
-        int dofIdxGlobal = this->dofMapper().index(entity);
-
-        if (!outStream.good())
-            DUNE_THROW(Dune::IOError, "Could not serialize entity " << dofIdxGlobal);
-
-        outStream << staticDat_[dofIdxGlobal].phasePresence << " ";
-    }
-
-    /*!
-     * \brief Reads the current solution from a restart file.
-     *
-     * \param inStream The input stream of one vertex from the restart file
-     * \param entity The entity, either a vertex or an element
-     */
-    template<class Entity>
-    void deserializeEntity(std::istream &inStream, const Entity &entity)
-    {
-        // read primary variables
-        ParentType::deserializeEntity(inStream, entity);
-
-        // read phase presence
-        int dofIdxGlobal = this->dofMapper().index(entity);
-
-        if (!inStream.good())
-            DUNE_THROW(Dune::IOError,
-                       "Could not deserialize entity " << dofIdxGlobal);
-
-        inStream >> staticDat_[dofIdxGlobal].phasePresence;
-        staticDat_[dofIdxGlobal].oldPhasePresence
-            = staticDat_[dofIdxGlobal].phasePresence;
-
-    }
-
-    /*!
-     * \brief Update the static data of all vertices in the grid.
-     *
-     * \param curGlobalSol The current global solution
-     * \param oldGlobalSol The previous global solution
-     */
-    void updateStaticData(SolutionVector &curGlobalSol,
-                          const SolutionVector &oldGlobalSol)
-    {
-        bool wasSwitched = false;
-        int succeeded;
-        try {
-        for (unsigned i = 0; i < staticDat_.size(); ++i)
-            staticDat_[i].visited = false;
-
-        FVElementGeometry fvGeometry;
-        static VolumeVariables volVars;
-        for (const auto& element : elements(this->gridView_()))
-        {
-            fvGeometry.update(this->gridView_(), element);
-            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-            {
-                int globalIdx = this->dofMapper().subIndex(element, scvIdx, dofCodim);
-
-                if (staticDat_[globalIdx].visited)
-                    continue;
-
-                staticDat_[globalIdx].visited = true;
-                volVars.update(curGlobalSol[globalIdx],
-                               this->problem_(),
-                               element,
-                               fvGeometry,
-                               scvIdx,
-                               false);
-                const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
-
-                    if (primaryVarSwitch_(curGlobalSol,
-                                          volVars,
-                                          globalIdx,
-                                          globalPos))
-                    {
-                        this->jacobianAssembler().markDofRed(globalIdx);
-                        wasSwitched = true;
-                    }
-            }
-        }
-         succeeded = 1;
-        }
-
-        catch (Dumux::NumericalProblem &e)
-        {
-            std::cout << "\n"
-                      << "Rank " << this->problem_().gridView().comm().rank()
-                      << " caught an exception while updating the static data." << e.what()
-                      << "\n";
-            succeeded = 0;
-        }
-
-        //make sure that all processes succeeded. If not throw a NumericalProblem to decrease the time step size.
-        if (this->gridView_().comm().size() > 1)
-            succeeded = this->gridView_().comm().min(succeeded);
-
-        if (!succeeded) {
-            updateFailed();
-
-            DUNE_THROW(NumericalProblem,
-                       "A process did not succeed in updating the static data.");
-
-            return;
-        }
-
-        // make sure that if there was a variable switch in an
-        // other partition we will also set the switch flag
-        // for our partition.
-        if (this->gridView_().comm().size() > 1)
-            wasSwitched = this->gridView_().comm().max(wasSwitched);
-
-        setSwitched_(wasSwitched);
-    }
-
-protected:
-    /*!
-     * \brief Data which is attached to each vertex and is not only
-     *        stored locally.
-     */
-    struct StaticVars
-    {
-        int phasePresence;
-        bool wasSwitched;
-
-        int oldPhasePresence;
-        bool visited;
-    };
-
-    /*!
-     * \brief Reset the current phase presence of all vertices to the old one.
-     *
-     * This is done after an update failed.
-     */
-    void resetPhasePresence_()
-    {
-        for (unsigned int i = 0; i < this->numDofs(); ++i)
-        {
-            staticDat_[i].phasePresence
-                = staticDat_[i].oldPhasePresence;
-            staticDat_[i].wasSwitched = false;
-        }
-    }
-
-    /*!
-     * \brief Set the old phase of all verts state to the current one.
-     */
-    void updateOldPhasePresence_()
-    {
-        for (unsigned int i = 0; i < this->numDofs(); ++i)
-        {
-            staticDat_[i].oldPhasePresence
-                = staticDat_[i].phasePresence;
-            staticDat_[i].wasSwitched = false;
-        }
-    }
-
-    /*!
-     * \brief Set whether there was a primary variable switch after in
-     *        the last timestep.
-     */
-    void setSwitched_(bool yesno)
-    {
-        switchFlag_ = yesno;
-    }
-
-    //  perform variable switch at a vertex; Returns true if a
-    //  variable switch was performed.
-    bool primaryVarSwitch_(SolutionVector &globalSol,
-                           const VolumeVariables &volVars,
-                           int globalIdx,
-                           const GlobalPosition &globalPos)
-    {
-        // evaluate primary variable switch
-        bool wouldSwitch = false;
-        int phasePresence = staticDat_[globalIdx].phasePresence;
-        int newPhasePresence = phasePresence;
-
-        globalSol[globalIdx][pressureIdx] = volVars.fluidState().pressure(gPhaseIdx);
-
-        // check if a primary var switch is necessary
-        if (phasePresence == twoPhases)
-        {
-            Scalar Smin = 0;
-            if (staticDat_[globalIdx].wasSwitched)
-                Smin = -0.01;
-
-            if (volVars.saturation(gPhaseIdx) <= Smin)
-            {
-                wouldSwitch = true;
-                // gas phase disappears
-                std::cout << "Gas phase disappears at vertex " << globalIdx
-                          << ", coordinates: " << globalPos << ", sg: "
-                          << volVars.saturation(gPhaseIdx) << std::endl;
-                newPhasePresence = wPhaseOnly;
-
-                globalSol[globalIdx][switch1Idx] = volVars.fluidState().temperature();
-            }
-            else if (volVars.saturation(wPhaseIdx) <= Smin)
-            {
-                wouldSwitch = true;
-                // water phase disappears
-                std::cout << "Wetting phase disappears at vertex " << globalIdx
-                          << ", coordinates: " << globalPos << ", sw: "
-                          << volVars.saturation(wPhaseIdx) << std::endl;
-                newPhasePresence = gPhaseOnly;
-
-                globalSol[globalIdx][switch1Idx] = volVars.fluidState().temperature();
-            }
-
-        }
-        else if (phasePresence == wPhaseOnly)
-        {
-            Scalar temp = volVars.fluidState().temperature();
-            Scalar tempVap = volVars.vaporTemperature();
-
-            // if the the temperature would be larger than
-            // the vapor temperature at the given pressure, gas phase appears
-            if (temp >= tempVap)
-            {
-                wouldSwitch = true;
-                // gas phase appears
-                std::cout << "gas phase appears at vertex " << globalIdx
-                          << ", coordinates: " << globalPos  << std::endl;
-
-               newPhasePresence = twoPhases;
-               globalSol[globalIdx][switch1Idx] = 0.9999; //wetting phase saturation
-            }
-        }
-
-
-        else if (phasePresence == gPhaseOnly)
-        {
-
-            Scalar temp = volVars.fluidState().temperature();
-            Scalar tempVap = volVars.vaporTemperature();
-
-            if (temp < tempVap)
-            {
-                wouldSwitch = true;
-                // wetting phase appears
-                std::cout << "wetting phase appears at vertex " << globalIdx
-                          << ", coordinates: " << globalPos  << std::endl;
-
-
-               newPhasePresence = twoPhases;
-               globalSol[globalIdx][switch1Idx] = 0.0001; //arbitrary small value
-            }
-    }
-        staticDat_[globalIdx].phasePresence = newPhasePresence;
-        staticDat_[globalIdx].wasSwitched = wouldSwitch;
-        return phasePresence != newPhasePresence;
-    }
-
-    // parameters given in constructor
-    std::vector<StaticVars> staticDat_;
-    bool switchFlag_;
-};
-
-}
-
-#include "propertydefaults.hh"
-
-#endif
diff --git a/dumux/porousmediumflow/2p1c/implicit/newtoncontroller.hh b/dumux/porousmediumflow/2p1c/implicit/newtoncontroller.hh
deleted file mode 100644
index bc5c7e32f493186f9dbdb23cdc9e0125f8132d20..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2p1c/implicit/newtoncontroller.hh
+++ /dev/null
@@ -1,86 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program 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 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- * \brief A 2p1c specific controller for the newton solver.
- *
- * This controller 'knows' what a 'physically meaningful' solution is
- * which allows the newton method to abort quicker if the solution is
- * way out of bounds.
- */
-#ifndef DUMUX_2P1CNI_NEWTON_CONTROLLER_HH
-#define DUMUX_2P1CNI_NEWTON_CONTROLLER_HH
-
-#include "properties.hh"
-
-#include <dumux/porousmediumflow/2p2c/implicit/newtoncontroller.hh>
-#include <dumux/nonlinear/newtoncontroller.hh>
-
-namespace Dumux {
-/*!
- * \ingroup TwoPOneCModel
- * \brief A 2p1cni specific controller for the newton solver.
- *
- * This controller 'knows' what a 'physically meaningful' solution is
- * which allows the newton method to abort quicker if the solution is
- * way out of bounds.
- */
-template <class TypeTag>
-class TwoPOneCNINewtonController : public TwoPTwoCNewtonController<TypeTag>
-{
-    typedef TwoPTwoCNewtonController<TypeTag> ParentType;
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
-
-public:
-    TwoPOneCNINewtonController(const Problem &problem)
-        : ParentType(problem)
-    {}
-
-    void newtonBeginStep()
-    {
-        // call the method of the base class
-        ParentType::newtonBeginStep();
-
-        //reset the indicator of spurious cold water fluxes into the steam zone
-        this->model_().localJacobian().localResidual().resetSpuriousFlowDetected();
-    }
-
-
-    /*!
-     * \brief Called after each Newton update
-     *
-     * \param uCurrentIter The current global solution vector
-     * \param uLastIter The previous global solution vector
-     */
-    void newtonEndStep(SolutionVector &uCurrentIter,
-                       const SolutionVector &uLastIter)
-    {
-        // call the method of the base class
-        ParentType::newtonEndStep(uCurrentIter, uLastIter);
-
-        //Output message after each Newton Step if spurious fluxes have been blocked
-        if (this->model_().localJacobian().localResidual().spuriousFlowDetected())
-            std::cout <<"Spurious fluxes blocked!" <<std::endl;
-    }
-
-};
-}
-
-#endif
diff --git a/dumux/porousmediumflow/2p1c/implicit/properties.hh b/dumux/porousmediumflow/2p1c/implicit/properties.hh
deleted file mode 100644
index 99ad0da677c20a8e53e1bec5066b36ff38162a8f..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2p1c/implicit/properties.hh
+++ /dev/null
@@ -1,69 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program 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 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \ingroup TwoPOneCModel
- */
-/*!
- * \file
- *
- * \brief Defines the properties required for the 2p1cn fully implicit model.
- *
- *Important note: The 2p1c model requires the use of the non-isothermal extension found in dumux/implicit/nonisothermal
- */
-#ifndef DUMUX_2P1C_PROPERTIES_HH
-#define DUMUX_2P1C_PROPERTIES_HH
-
-#include <dumux/implicit/box/properties.hh>
-#include <dumux/implicit/cellcentered/properties.hh>
-#include <dumux/porousmediumflow/nonisothermal/implicit/model.hh>
-
-namespace Dumux
-{
-
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Type tags
-//////////////////////////////////////////////////////////////////
-NEW_TYPE_TAG(TwoPOneCNI, INHERITS_FROM(NonIsothermal));
-NEW_TYPE_TAG(BoxTwoPOneCNI, INHERITS_FROM(BoxModel, TwoPOneCNI));
-NEW_TYPE_TAG(CCTwoPOneCNI, INHERITS_FROM(CCModel, TwoPOneCNI));
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-
-NEW_PROP_TAG(NumPhases);   //!< Number of fluid phases in the system
-NEW_PROP_TAG(NumComponents); //!< Number of fluid components in the system
-NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters
-NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
-
-NEW_PROP_TAG(MaterialLaw);   //!< The material law which ought to be used (extracted from the spatial parameters)
-
-NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
-
-NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the upwind parameter for the mobility
-NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation
-NEW_PROP_TAG(UseBlockingOfSpuriousFlow); //!< Determines whether Blocking ofspurious flow is used
-NEW_PROP_TAG(BaseFluxVariables); //! The base flux variables
-NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
-}
-}
-
-#endif
diff --git a/dumux/porousmediumflow/2p1c/implicit/propertydefaults.hh b/dumux/porousmediumflow/2p1c/implicit/propertydefaults.hh
deleted file mode 100644
index 19b589d5ae175bd5d915d437eb15ac8a064782ca..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2p1c/implicit/propertydefaults.hh
+++ /dev/null
@@ -1,171 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program 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 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \ingroup TwoPOneCModel
- */
-/*!
- * \file
- *
- * \brief Defines default values for most properties required by the
- *        2p1cni fully implicit model.
- *
- *Important note: The 2p1c model requires the use of the non-isothermal extension found in dumux/implicit/nonisothermal
- */
-#ifndef DUMUX_2P1C_PROPERTY_DEFAULTS_HH
-#define DUMUX_2P1C_PROPERTY_DEFAULTS_HH
-
-#include "model.hh"
-#include "indices.hh"
-#include "volumevariables.hh"
-#include "properties.hh"
-#include "newtoncontroller.hh"
-#include "localresidual.hh"
-
-#include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
-#include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
-#include <dumux/material/fluidsystems/2pliquidvapor.hh>
-#include <dumux/material/components/constant.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
-#include <dumux/material/spatialparams/implicit.hh>
-
-namespace Dumux
-{
-
-namespace Properties {
-//////////////////////////////////////////////////////////////////
-// Property values
-//////////////////////////////////////////////////////////////////
-
-/*!
- * \brief Set the property for the number of components.
- *
- * We just forward the number from the fluid system and use an static
- * assert to make sure it is 1.
- */
-SET_PROP(TwoPOneCNI, NumComponents)
-{
- private:
-    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-
- public:
-    static const int value = FluidSystem::numComponents;
-
-    static_assert(value == 1,
-                  "Only fluid systems with 1 component are supported by the 2p1cni model!");
-};
-
-/*!
- * \brief Set the property for the number of fluid phases.
- *
- * We just forward the number from the fluid system and use an static
- * assert to make sure it is 2.
- */
-SET_PROP(TwoPOneCNI, NumPhases)
-{
- private:
-    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-
- public:
-    static const int value = FluidSystem::numPhases;
-    static_assert(value == 2,
-                  "Only fluid systems with 2 phases are supported by the 2p1cni model!");
-};
-
-//! Use the 2p2c specific newton controller for dealing with phase switches
-SET_TYPE_PROP(TwoPOneCNI, NewtonController, TwoPOneCNINewtonController<TypeTag>);
-
-//! the FluxVariables property
-// SET_TYPE_PROP(TwoPOneCNI, FluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
-
-//! the upwind factor for the mobility.
-SET_SCALAR_PROP(TwoPOneCNI, ImplicitMassUpwindWeight, 1.0);
-
-//! set default mobility upwind weight to 1.0, i.e. fully upwind
-SET_SCALAR_PROP(TwoPOneCNI, ImplicitMobilityUpwindWeight, 1.0);
-
-//! The fluid system to use by default
-SET_PROP(TwoPOneCNI, FluidSystem)
-{ private:
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-public:
-    typedef FluidSystems::TwoPLiquidVaporFluidsystem<Scalar, Components::Constant<1, Scalar> > type;
-};
-
-//! Determines whether Blocking ofspurious flow is used
-SET_BOOL_PROP(TwoPOneCNI, UseBlockingOfSpuriousFlow, false);
-
-//! The spatial parameters to be employed.
-//! Use ImplicitSpatialParams by default.
-SET_TYPE_PROP(TwoPOneCNI, SpatialParams, ImplicitSpatialParams<TypeTag>);
-
-// disable velocity output by default
-
-// enable gravity by default
-SET_BOOL_PROP(TwoPOneCNI, ProblemEnableGravity, true);
-
-//! default value for the forchheimer coefficient
-// Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90.
-//        Actually the Forchheimer coefficient is also a function of the dimensions of the
-//        porous medium. Taking it as a constant is only a first approximation
-//        (Nield, Bejan, Convection in porous media, 2006, p. 10)
-SET_SCALAR_PROP(BoxModel, SpatialParamsForchCoeff, 0.55);
-
-
-//! Somerton is used as default model to compute the effective thermal heat conductivity
-SET_PROP(NonIsothermal, ThermalConductivityModel)
-{
-private:
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-public:
-    typedef ThermalConductivitySomerton<Scalar> type;
-};
-
-//////////////////////////////////////////////////////////////////
-// Property values for isothermal model required for the general non-isothermal model
-//////////////////////////////////////////////////////////////////
-
-// set isothermal Model
-SET_TYPE_PROP(TwoPOneCNI, IsothermalModel, TwoPOneCModel<TypeTag>);
-
-// set isothermal FluxVariables
-SET_TYPE_PROP(TwoPOneCNI, IsothermalFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
-
-//set isothermal VolumeVariables
-SET_TYPE_PROP(TwoPOneCNI, IsothermalVolumeVariables, TwoPOneCVolumeVariables<TypeTag>);
-
-//set isothermal LocalResidual
-SET_TYPE_PROP(TwoPOneCNI, IsothermalLocalResidual, TwoPOneCLocalResidual<TypeTag>);
-
-//set isothermal Indices
-SET_PROP(TwoPOneCNI, IsothermalIndices)
-{
-
-public:
-    typedef TwoPOneCIndices<TypeTag, 0> type;
-};
-
-//set isothermal NumEq
-SET_INT_PROP(TwoPOneCNI, IsothermalNumEq, 1);
-
-}
-
-}
-
-#endif
diff --git a/dumux/porousmediumflow/2p1c/implicit/indices.hh b/dumux/porousmediumflow/2p1c/indices.hh
similarity index 89%
rename from dumux/porousmediumflow/2p1c/implicit/indices.hh
rename to dumux/porousmediumflow/2p1c/indices.hh
index 49a07a06e8fea6035c1e9e8d98799e83212b734b..6bd3e53b970df39299c7cdc9e88174b6d1474dd3 100644
--- a/dumux/porousmediumflow/2p1c/implicit/indices.hh
+++ b/dumux/porousmediumflow/2p1c/indices.hh
@@ -25,11 +25,22 @@
 #ifndef DUMUX_2P1C_INDICES_HH
 #define DUMUX_2P1C_INDICES_HH
 
-#include "properties.hh"
-
 namespace Dumux
 {
 
+/*!
+ * \ingroup TwoPNCModel
+ * \ingroup ImplicitIndices
+ * \brief Enumerates the formulations which the two-phase n-component model accepts.
+ */
+struct TwoPOneCFormulation
+{
+    enum {
+        pnsw,
+        pwsn
+        };
+};
+
 /*!
  * \ingroup TwoPOneCModel
  * \ingroup ImplicitIndices
@@ -46,12 +57,12 @@ class TwoPOneCIndices
 public:
     // Phase indices
     static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< index of the wetting liquid phase
-    static const int gPhaseIdx = FluidSystem::nPhaseIdx; //!< index of the gas phase
+    static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< index of the gas phase
 
     // present phases (-> 'pseudo' primary variable)
     static const int twoPhases = 1; //!< All three phases are present
     static const int wPhaseOnly = 2; //!< Only the water phase is present
-    static const int gPhaseOnly = 3; //!< Only gas phase is present
+    static const int nPhaseOnly = 3; //!< Only gas phase is present
 
     // Primary variable indices
     static const int pressureIdx = PVOffset + 0; //!< Index for phase pressure in a solution vector
diff --git a/dumux/porousmediumflow/2p1c/localresidual.hh b/dumux/porousmediumflow/2p1c/localresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..de98aa2c269ff54df7939b3941d3abdf4af1ac0a
--- /dev/null
+++ b/dumux/porousmediumflow/2p1c/localresidual.hh
@@ -0,0 +1,132 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program 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 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Element-wise calculation of the residual for problems
+ *        using the two-phase one-component fully implicit models.
+ */
+#ifndef DUMUX_2P1C_LOCAL_RESIDUAL_HH
+#define DUMUX_2P1C_LOCAL_RESIDUAL_HH
+
+#include <dumux/porousmediumflow/immiscible/localresidual.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPOneC
+ * \ingroup ImplicitLocalResidual
+ * \brief Element-wise calculation of the residual for problems
+ *        using the two-phase one-component fully implicit models.
+ */
+template<class TypeTag>
+class TwoPOneCLocalResidual : public ImmiscibleLocalResidual<TypeTag>
+{
+    using ParentType = ImmiscibleLocalResidual<TypeTag>;
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+    using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+    using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    // first index for the mass balance
+    enum { conti0EqIdx = Indices::conti0EqIdx };
+
+    static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+
+public:
+    using ParentType::ParentType;
+
+    /*!
+     * \brief Evaluate the rate of change of all conservation
+     *        quantites (e.g. phase mass) within a sub-control
+     *        volume of a finite volume element for the immiscible models.
+     * \param scv The sub control volume
+     * \param volVars The current or previous volVars
+     * \note This function should not include the source and sink terms.
+     * \note The volVars can be different to allow computing
+     *       the implicit euler time derivative here
+     */
+    ResidualVector computeStorage(const Problem& problem,
+                                  const SubControlVolume& scv,
+                                  const VolumeVariables& volVars) const
+    {
+        ResidualVector storage(0.0);
+        // compute storage term of all components within all phases
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+        {
+            storage[conti0EqIdx] +=
+                volVars.porosity()
+                * volVars.saturation(phaseIdx) * volVars.density(phaseIdx);
+
+            EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
+        }
+
+        //! The energy storage in the solid matrix
+        EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
+
+        return storage;
+    }
+
+    /*!
+     * \brief Evaluate the mass flux over a face of a sub control volume
+     * \param scvf The sub control volume face to compute the flux on
+     */
+    ResidualVector computeFlux(const Problem& problem,
+                               const Element& element,
+                               const FVElementGeometry& fvGeometry,
+                               const ElementVolumeVariables& elemVolVars,
+                               const SubControlVolumeFace& scvf,
+                               const ElementFluxVariablesCache& elemFluxVarsCache) const
+    {
+        FluxVariables fluxVars;
+        fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+
+        ResidualVector flux;
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+        {
+            // the physical quantities for which we perform upwinding
+            auto upwindTerm = [phaseIdx](const auto& volVars)
+                              { return volVars.density(phaseIdx)*volVars.mobility(phaseIdx); };
+
+            flux[conti0EqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
+
+            //! Add advective phase energy fluxes. For isothermal model the contribution is zero.
+            EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
+        }
+
+        //! Add diffusive energy fluxes. For isothermal model the contribution is zero.
+        EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
+
+        return flux;
+    }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/porousmediumflow/2p1c/model.hh b/dumux/porousmediumflow/2p1c/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7e0178d9c3ac186858d4d810d94e0f5b1514e8ce
--- /dev/null
+++ b/dumux/porousmediumflow/2p1c/model.hh
@@ -0,0 +1,201 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program 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 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Adaption of the fully implicit scheme to the two-phase one-component flow model.
+ *
+ */
+#ifndef DUMUX_2P1C_MODEL_HH
+#define DUMUX_2P1C_MODEL_HH
+
+#include <dumux/common/properties.hh>
+
+#include <dumux/material/spatialparams/implicit.hh>
+#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
+#include <dumux/material/fluidstates/immiscible.hh>
+#include <dumux/material/fluidstates/compositional.hh>
+
+
+#include <dumux/porousmediumflow/properties.hh>
+#include <dumux/porousmediumflow/immiscible/localresidual.hh>
+#include <dumux/porousmediumflow/compositional/switchableprimaryvariables.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/model.hh>
+
+#include <dumux/porousmediumflow/2p/implicit/vtkoutputfields.hh>
+
+#include "darcyslaw.hh"
+#include "vtkoutputfields.hh"
+#include "localresidual.hh"
+#include "indices.hh"
+#include "volumevariables.hh"
+#include "primaryvariableswitch.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPOneCModel
+ * \brief Adaption of the fully implicit scheme to the two-phase one-component flow model.
+ *
+ * \note The 2p1c model requires the use of the non-isothermal extension found in dumux/implicit/nonisothermal.
+ *
+ * This model is designed for simulating two fluid phases with water as the only component.
+ * It is particularly suitable for the simulation of steam injection in saturated conditions.
+ *
+ * The model implements the flow of two phases and one component, i.e. a pure liquid (e.g. water)
+ * and its vapor (e.g. steam),
+ * \f$\alpha \in \{ w, n \}\f$ using a standard multiphase Darcy
+ * approach as the equation for the conservation of momentum, i.e.
+ \f[
+ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K}
+ \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right)
+ \f]
+ *
+ * By inserting this into the equation for the conservation of the
+ * phase mass, one gets
+ \f[
+\phi \frac{\partial\ \sum_\alpha (\rho_\alpha S_\alpha)}{\partial t} \\-\sum \limits_ \alpha \text{div} \left \{\rho_\alpha \frac{k_{r\alpha}}{\mu_\alpha}
+\mathbf{K} (\mathbf{grad}p_\alpha - \rho_\alpha \mathbf{g}) \right \} -q^w =0
+ \f]
+ *
+ * All equations are discretized using a vertex-centered finite volume (box)
+ * or cell-centered finite volume scheme as spatial
+ * and the implicit Euler method as time discretization.
+ *
+ * By using constitutive relations for the capillary pressure \f$p_c =
+ * p_n - p_w\f$ and relative permeability \f$k_{r\alpha}\f$ and taking
+ * advantage of the fact that \f$S_w + S_n = 1\f$, the number of
+ * unknowns can be reduced to two. The model features a primary variable switch.
+ * If only one phase is present, \f$p_g\f$ and \f$T\f$ are the primary variables.
+ * In the presence of two phases, \f$p_g\f$ and \f$S_w\f$ become the new primary variables.
+ */
+
+namespace Properties
+{
+ //////////////////////////////////////////////////////////////////
+ // Type tags
+ //////////////////////////////////////////////////////////////////
+ NEW_TYPE_TAG(TwoPOneCNI, INHERITS_FROM(PorousMediumFlow, NonIsothermal));
+
+ //////////////////////////////////////////////////////////////////
+ // Property tags
+ //////////////////////////////////////////////////////////////////
+
+ NEW_PROP_TAG(UseBlockingOfSpuriousFlow); //!< Determines whether Blocking ofspurious flow is used
+
+ /*!
+  * \brief Set the property for the number of components.
+  *
+  * We just forward the number from the fluid system and use an static
+  * assert to make sure it is 1.
+  */
+ SET_PROP(TwoPOneCNI, NumComponents)
+ {
+  private:
+     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+  public:
+     static constexpr auto value = FluidSystem::numComponents;
+
+     static_assert(value == 1,
+                   "Only fluid systems with 1 component are supported by the 2p1cni model!");
+ };
+
+ /*!
+  * \brief Set the property for the number of fluid phases.
+  *
+  * We just forward the number from the fluid system and use an static
+  * assert to make sure it is 2.
+  */
+ SET_PROP(TwoPOneCNI, NumPhases)
+ {
+  private:
+     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+  public:
+     static constexpr auto value = FluidSystem::numPhases;
+     static_assert(value == 2,
+                   "Only fluid systems with 2 phases are supported by the 2p1cni model!");
+ };
+
+SET_TYPE_PROP(TwoPOneCNI, LocalResidual, TwoPOneCLocalResidual<TypeTag>); //! The local residual function
+
+SET_BOOL_PROP(TwoPOneCNI, EnableAdvection, true);                           //! The one-phase model considers advection
+SET_BOOL_PROP(TwoPOneCNI, EnableMolecularDiffusion, false);                 //! The one-phase model has no molecular diffusion
+
+ /*!
+  * \brief The fluid state which is used by the volume variables to
+  *        store the thermodynamic state. This should be chosen
+  *        appropriately for the model ((non-)isothermal, equilibrium, ...).
+  *        This can be done in the problem.
+  */
+ SET_PROP(TwoPOneCNI, FluidState)
+ {
+ private:
+     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ public:
+     using type = CompositionalFluidState<Scalar, FluidSystem>;
+    //  using type = ImmiscibleFluidState<Scalar, FluidSystem>;
+ };
+
+ //! Determines whether Blocking ofspurious flow is used
+ SET_BOOL_PROP(TwoPOneCNI, UseBlockingOfSpuriousFlow, false);
+
+ SET_TYPE_PROP(TwoPOneCNI, AdvectionType, TwoPOneCDarcysLaw<TypeTag>);
+
+  SET_TYPE_PROP(TwoPOneCNI, VolumeVariables, TwoPOneCVolumeVariables<TypeTag>);
+
+ //! The primary variable switch for the 2p1c model
+ SET_TYPE_PROP(TwoPOneCNI, PrimaryVariableSwitch, TwoPOneCPrimaryVariableSwitch<TypeTag>);
+
+ //! The primary variables vector for the 2p1c model
+ SET_TYPE_PROP(TwoPOneCNI, PrimaryVariables, SwitchablePrimaryVariables<TypeTag, int>);
+
+ //! Somerton is used as default model to compute the effective thermal heat conductivity
+ SET_PROP(TwoPOneCNI, ThermalConductivityModel)
+ {
+ private:
+     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ public:
+     using type = ThermalConductivitySomerton<Scalar>;
+ };
+
+ //////////////////////////////////////////////////////////////////
+ // Property values for isothermal model required for the general non-isothermal model
+ //////////////////////////////////////////////////////////////////
+ //set isothermal VolumeVariables
+ SET_TYPE_PROP(TwoPOneCNI, IsothermalVolumeVariables, TwoPOneCVolumeVariables<TypeTag>);
+
+ //set isothermal LocalResidual
+ SET_TYPE_PROP(TwoPOneCNI, IsothermalLocalResidual, ImmiscibleLocalResidual<TypeTag>);
+
+ //set isothermal Indices
+ SET_TYPE_PROP(TwoPOneCNI, IsothermalIndices, TwoPOneCIndices<TypeTag, 0>);
+
+//! the isothermal vtk output fields
+ SET_TYPE_PROP(TwoPOneCNI, IsothermalVtkOutputFields, TwoPOneCVtkOutputFields<TypeTag>);
+
+ //set isothermal NumEq
+ SET_INT_PROP(TwoPOneCNI, IsothermalNumEq, 1);
+
+} // end namespace Properties
+} // end namespace Dumux
+
+#endif // DUMUX_2P1C_MODEL_HH
diff --git a/dumux/porousmediumflow/2p1c/primaryvariableswitch.hh b/dumux/porousmediumflow/2p1c/primaryvariableswitch.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1cb87b50f522c25b497199a3501876f37760aeda
--- /dev/null
+++ b/dumux/porousmediumflow/2p1c/primaryvariableswitch.hh
@@ -0,0 +1,156 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program 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 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief The primary variable switch for the 2p2c model
+ */
+#ifndef DUMUX_2P1C_PRIMARY_VARIABLE_SWITCH_HH
+#define DUMUX_2P1C_PRIMARY_VARIABLE_SWITCH_HH
+
+#include <dumux/porousmediumflow/compositional/primaryvariableswitch.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPTwoCModel
+ * \brief The primary variable switch controlling the phase presence state variable
+ */
+template<class TypeTag>
+class TwoPOneCPrimaryVariableSwitch : public PrimaryVariableSwitch<TypeTag>
+{
+    friend typename Dumux::PrimaryVariableSwitch<TypeTag>;
+    using ParentType = Dumux::PrimaryVariableSwitch<TypeTag>;
+
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using IndexType = typename GridView::IndexSet::IndexType;
+    using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
+
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+    enum {
+        switchIdx = Indices::switch1Idx,
+
+        pressureIdx = Indices::pressureIdx,
+
+        wPhaseIdx = Indices::wPhaseIdx,
+        nPhaseIdx = Indices::nPhaseIdx,
+
+        twoPhases = Indices::twoPhases,
+        wPhaseOnly  = Indices::wPhaseOnly,
+        nPhaseOnly  = Indices::nPhaseOnly
+    };
+
+public:
+    using ParentType::ParentType;
+
+protected:
+
+    // perform variable switch at a degree of freedom location
+    bool update_(PrimaryVariables& priVars,
+                 const VolumeVariables& volVars,
+                 IndexType dofIdxGlobal,
+                 const GlobalPosition& globalPos)
+    {
+        // evaluate primary variable switch
+        bool wouldSwitch = false;
+        int phasePresence =  priVars.state();
+        int newPhasePresence = phasePresence;
+
+        // check if a primary var switch is necessary
+        if (phasePresence == twoPhases)
+        {
+            Scalar Smin = 0;
+            if (this->wasSwitched_[dofIdxGlobal])
+                Smin = -0.01;
+
+            if (volVars.saturation(nPhaseIdx) <= Smin)
+            {
+                wouldSwitch = true;
+                // gas phase disappears
+                std::cout << "Gas phase disappears at vertex " << dofIdxGlobal
+                          << ", coordinates: " << globalPos << ", sg: "
+                          << volVars.saturation(nPhaseIdx) << std::endl;
+                newPhasePresence = wPhaseOnly;
+
+                priVars[switchIdx] = volVars.fluidState().temperature();
+            }
+            else if (volVars.saturation(wPhaseIdx) <= Smin)
+            {
+                wouldSwitch = true;
+                // water phase disappears
+                std::cout << "Wetting phase disappears at vertex " << dofIdxGlobal
+                          << ", coordinates: " << globalPos << ", sw: "
+                          << volVars.saturation(wPhaseIdx) << std::endl;
+                newPhasePresence = nPhaseOnly;
+
+                priVars[switchIdx] = volVars.fluidState().temperature();
+            }
+
+        }
+        else if (phasePresence == wPhaseOnly)
+        {
+            const Scalar temp = volVars.fluidState().temperature();
+            const Scalar tempVap = volVars.vaporTemperature();
+
+            // if the the temperature would be larger than
+            // the vapor temperature at the given pressure, gas phase appears
+            if (temp >= tempVap)
+            {
+                wouldSwitch = true;
+                // gas phase appears
+                std::cout << "gas phase appears at vertex " << dofIdxGlobal
+                          << ", coordinates: " << globalPos  << std::endl;
+
+               newPhasePresence = twoPhases;
+               priVars[switchIdx] = 0.9999; //wetting phase saturation
+            }
+        }
+
+
+        else if (phasePresence == nPhaseOnly)
+        {
+
+            const Scalar temp = volVars.fluidState().temperature();
+            const Scalar tempVap = volVars.vaporTemperature();
+
+            if (temp < tempVap)
+            {
+                wouldSwitch = true;
+                // wetting phase appears
+                std::cout << "wetting phase appears at vertex " << dofIdxGlobal
+                          << ", coordinates: " << globalPos  << std::endl;
+
+
+               newPhasePresence = twoPhases;
+               priVars[switchIdx] = 0.0001; //arbitrary small value
+            }
+    }
+        priVars.setState(newPhasePresence);
+        this->wasSwitched_[dofIdxGlobal] = wouldSwitch;
+        return phasePresence != newPhasePresence;
+    }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/porousmediumflow/2p1c/implicit/volumevariables.hh b/dumux/porousmediumflow/2p1c/volumevariables.hh
similarity index 59%
rename from dumux/porousmediumflow/2p1c/implicit/volumevariables.hh
rename to dumux/porousmediumflow/2p1c/volumevariables.hh
index cf406f1ff303196525a8db14fbfe9f77f95760e2..24b9949ac61bd91b162c747fef0da723b7b96160 100644
--- a/dumux/porousmediumflow/2p1c/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/2p1c/volumevariables.hh
@@ -27,35 +27,37 @@
 #ifndef DUMUX_2P1C_VOLUME_VARIABLES_HH
 #define DUMUX_2P1C_VOLUME_VARIABLES_HH
 
-#include <dumux/implicit/volumevariables.hh>
+#include <dumux/porousmediumflow/volumevariables.hh>
 #include <dumux/material/fluidstates/compositional.hh>
-#include "properties.hh"
+#include "indices.hh"
 
 namespace Dumux
 {
 
 /*!
  * \ingroup TwoPOneCModel
- * \ingroup ImplicitVolumeVariables
+ * \ingroup PorousmediumflowVolumeVariables
  * \brief Contains the quantities which are are constant within a
  *        finite volume in the two-phase, two-component model.
  */
 template <class TypeTag>
-class TwoPOneCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
+class TwoPOneCVolumeVariables : public PorousMediumFlowVolumeVariables<TypeTag>
 {
-    typedef ImplicitVolumeVariables<TypeTag> ParentType;
-    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) Implementation;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
-    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params MaterialLawParams;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    using ParentType = PorousMediumFlowVolumeVariables<TypeTag>;
+    using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
+    using PermeabilityType = typename SpatialParams::PermeabilityType;
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+    using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params;
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
     enum {
         dim = GridView::dimension,
 
@@ -63,7 +65,7 @@ class TwoPOneCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
         numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
 
         wPhaseIdx = Indices::wPhaseIdx,
-        gPhaseIdx = Indices::gPhaseIdx,
+        nPhaseIdx = Indices::nPhaseIdx,
 
         switch1Idx = Indices::switch1Idx,
         pressureIdx = Indices::pressureIdx
@@ -73,42 +75,70 @@ class TwoPOneCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
     enum {
         twoPhases = Indices::twoPhases,
         wPhaseOnly  = Indices::wPhaseOnly,
-        gPhaseOnly  = Indices::gPhaseOnly,
+        nPhaseOnly  = Indices::nPhaseOnly,
     };
 
-    typedef typename GridView::template Codim<0>::Entity Element;
-    typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix;
-
-    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
-    enum { dofCodim = isBox ? dim : 0 };
+    using Element = typename GridView::template Codim<0>::Entity;
 
 public:
     //! The type of the object returned by the fluidState() method
-    typedef Dumux::CompositionalFluidState<Scalar, FluidSystem> FluidState;
+    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
 
     /*!
      * \copydoc ImplicitVolumeVariables::update
      */
-    void update(const PrimaryVariables &priVars,
+    void update(const ElementSolutionVector &elemSol,
                 const Problem &problem,
                 const Element &element,
-                const FVElementGeometry &fvGeometry,
-                const int scvIdx,
-                bool isOldSol)
+                const SubControlVolume& scv)
+    {
+        ParentType::update(elemSol, problem, element, scv);
+
+        completeFluidState(elemSol, problem, element, scv, fluidState_);
+
+        /////////////
+        // calculate the remaining quantities
+        /////////////
+        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
+
+        // Second instance of a parameter cache.
+        // Could be avoided if diffusion coefficients also
+        // became part of the fluid state.
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updateAll(fluidState_);
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+        {
+            // relative permeabilities
+            Scalar kr;
+            if (phaseIdx == wPhaseIdx)
+                kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx));
+            else // ATTENTION: krn requires the wetting phase saturation
+                // as parameter!
+                kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx));
+            relativePermeability_[phaseIdx] = kr;
+            Valgrind::CheckDefined(relativePermeability_[phaseIdx]);
+        }
+
+        // porosity & permeability
+        porosity_ = problem.spatialParams().porosity(element, scv, elemSol);
+        permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
+    }
+
+    /*!
+     * \copydoc ImplicitModel::completeFluidState
+     */
+    static void completeFluidState(const ElementSolutionVector& elemSol,
+                                   const Problem& problem,
+                                   const Element& element,
+                                   const SubControlVolume& scv,
+                                   FluidState& fluidState)
     {
-        ParentType::update(priVars,
-                           problem,
-                           element,
-                           fvGeometry,
-                           scvIdx,
-                           isOldSol);
 
-        // capillary pressure parameters
-        const MaterialLawParams &materialParams =
-            problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
+        // // capillary pressure parameters
+        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
 
-        int globalIdx = problem.model().dofMapper().subIndex(element, scvIdx, dofCodim);
-        int phasePresence = problem.model().phasePresence(globalIdx, isOldSol);
+        const auto& priVars = ParentType::extractDofPriVars(elemSol, scv);
+        const auto phasePresence = priVars.state();
 
         // get saturations
         Scalar sw(0.0);
@@ -123,7 +153,7 @@ public:
             sw = 1.0;
             sg = 0.0;
         }
-        else if (phasePresence == gPhaseOnly)
+        else if (phasePresence == nPhaseOnly)
         {
             sw = 0.0;
             sg = 1.0;
@@ -131,8 +161,8 @@ public:
         else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
         Valgrind::CheckDefined(sg);
 
-        fluidState_.setSaturation(wPhaseIdx, sw);
-        fluidState_.setSaturation(gPhaseIdx, sg);
+        fluidState.setSaturation(wPhaseIdx, sw);
+        fluidState.setSaturation(nPhaseIdx, sg);
 
         // get gas phase pressure
         const Scalar pg = priVars[pressureIdx];
@@ -144,64 +174,45 @@ public:
         const Scalar pw = pg - pc;
 
         //set pressures
-        fluidState_.setPressure(wPhaseIdx, pw);
-        fluidState_.setPressure(gPhaseIdx, pg);
+        fluidState.setPressure(wPhaseIdx, pw);
+        fluidState.setPressure(nPhaseIdx, pg);
 
         // get temperature
         Scalar temperature;
-        if (phasePresence == wPhaseOnly || phasePresence == gPhaseOnly)
+        if (phasePresence == wPhaseOnly || phasePresence == nPhaseOnly)
             temperature = priVars[switch1Idx];
         else if (phasePresence == twoPhases)
-            temperature = FluidSystem::vaporTemperature(fluidState_, wPhaseIdx);
+            temperature = FluidSystem::vaporTemperature(fluidState, wPhaseIdx);
         else
             DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
 
         Valgrind::CheckDefined(temperature);
 
-        fluidState_.setTemperature(temperature);
+        fluidState.setTemperature(temperature);
 
         // set the densities
-        const Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
-        const Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
+        const Scalar rhoW = FluidSystem::density(fluidState, wPhaseIdx);
+        const Scalar rhoG = FluidSystem::density(fluidState, nPhaseIdx);
 
-        fluidState_.setDensity(wPhaseIdx, rhoW);
-        fluidState_.setDensity(gPhaseIdx, rhoG);
+        fluidState.setDensity(wPhaseIdx, rhoW);
+        fluidState.setDensity(nPhaseIdx, rhoG);
 
         //get the viscosity and mobility
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
         {
             // Mobilities
             const Scalar mu =
-                FluidSystem::viscosity(fluidState_,
+                FluidSystem::viscosity(fluidState,
                                        phaseIdx);
-            fluidState_.setViscosity(phaseIdx,mu);
-
-            Scalar kr;
-            if (phaseIdx == wPhaseIdx)
-                kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx));
-            else // ATTENTION: krn requires the wetting phase saturation
-                // as parameter!
-                kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx));
-
-            mobility_[phaseIdx] = kr / mu;
-            Valgrind::CheckDefined(mobility_[phaseIdx]);
+            fluidState.setViscosity(phaseIdx,mu);
         }
 
-        // porosity
-        porosity_ = problem.spatialParams().porosity(element,
-                                                         fvGeometry,
-                                                         scvIdx);
-        Valgrind::CheckDefined(porosity_);
-
         // the enthalpies (internal energies are directly calculated in the fluidstate
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
         {
-            Scalar h = FluidSystem::enthalpy(fluidState_, phaseIdx);
-            fluidState_.setEnthalpy(phaseIdx, h);
+            const Scalar h = FluidSystem::enthalpy(fluidState, phaseIdx);
+            fluidState.setEnthalpy(phaseIdx, h);
         }
-
-        // energy related quantities not contained in the fluid state
-        asImp_().updateEnergy_(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
     }
 
     /*!
@@ -253,8 +264,8 @@ public:
      * temperatures of the rock matrix and of all fluid phases are
      * identical.
      */
-    Scalar temperature() const
-    { return fluidState_.temperature(/*phaseIdx=*/0); }
+    Scalar temperature(const int phaseIdx = 0) const
+    { return fluidState_.temperature(phaseIdx); }
 
     /*!
      * \brief Returns the effective mobility of a given phase within
@@ -264,14 +275,15 @@ public:
      */
     Scalar mobility(const int phaseIdx) const
     {
-        return mobility_[phaseIdx];
+        return relativePermeability_[phaseIdx]/fluidState_.viscosity(phaseIdx);
     }
 
     /*!
-     * \brief Returns the effective capillary pressure within the control volume.
+     * \brief Returns the effective capillary pressure within the control volume
+     *        in \f$[kg/(m*s^2)=N/m^2=Pa]\f$.
      */
     Scalar capillaryPressure() const
-    { return fluidState_.capillaryPressure(); }
+    { return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx); }
 
     /*!
      * \brief Returns the average porosity within the control volume.
@@ -279,6 +291,12 @@ public:
     Scalar porosity() const
     { return porosity_; }
 
+    /*!
+     * \brief Returns the average permeability within the control volume in \f$[m^2]\f$.
+     */
+    const PermeabilityType& permeability() const
+    { return permeability_; }
+
     /*!
      * \brief Returns the vapor temperature (T_{vap}(p_g) of the fluid within the control volume.
      */
@@ -289,9 +307,9 @@ protected:
     FluidState fluidState_;
 
 private:
-
-    Scalar porosity_;
-    Scalar mobility_[numPhases];
+    Scalar porosity_;               //!< Effective porosity within the control volume
+    PermeabilityType permeability_; //!> Effective permeability within the control volume
+    std::array<Scalar, numPhases> relativePermeability_;
 
     Implementation &asImp_()
     { return *static_cast<Implementation*>(this); }
diff --git a/dumux/porousmediumflow/2p1c/vtkoutputfields.hh b/dumux/porousmediumflow/2p1c/vtkoutputfields.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ee2702256c58fa2d6aaa30df9c93993d5710cb02
--- /dev/null
+++ b/dumux/porousmediumflow/2p1c/vtkoutputfields.hh
@@ -0,0 +1,53 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program 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 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Adds vtk output fields specific to the twop-onec model
+ */
+#ifndef DUMUX_TWOP_OneC_VTK_OUTPUT_FIELDS_HH
+#define DUMUX_TWOP_OneC_VTK_OUTPUT_FIELDS_HH
+
+#include <dumux/porousmediumflow/2p/implicit/vtkoutputfields.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPOneC, InputOutput
+ * \brief Adds vtk output fields specific to the TwoPOneC model
+ */
+template<class TypeTag>
+class TwoPOneCVtkOutputFields
+{
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+public:
+    template <class VtkOutputModule>
+    static void init(VtkOutputModule& vtk)
+    {
+        // use default fields from the 2p model
+        TwoPVtkOutputFields<TypeTag>::init(vtk);
+
+        //output additional to TwoP output:
+        vtk.addVolumeVariable([](const VolumeVariables& v){ return v.priVars().state(); }, "phasePresence");
+    }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/2p1c/implicit/CMakeLists.txt b/test/porousmediumflow/2p1c/implicit/CMakeLists.txt
index 436aaaac55c891152be6a21e14741cf0a592031a..4f936a7af39acd4f6c9c9c73586cb50bbabc2038 100644
--- a/test/porousmediumflow/2p1c/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/2p1c/implicit/CMakeLists.txt
@@ -1,16 +1,20 @@
 #add links to input files
 add_input_file_links()
 
-add_dumux_test(test_boxsteaminjection test_boxsteaminjection test_boxsteaminjection.cc
-               python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
-                 --script fuzzy
-                 --files ${CMAKE_SOURCE_DIR}/test/references/steaminjectionbox-reference.vtu
-                         ${CMAKE_CURRENT_BINARY_DIR}/test_boxsteaminjection-00008.vtu
-                 --command "${CMAKE_CURRENT_BINARY_DIR}/test_boxsteaminjection")
+dune_add_test(NAME test_2p1cni_box
+              SOURCES test_2p1c_fv.cc
+              COMPILE_DEFINITIONS TYPETAG=TwoPOneCNIBoxTypeTag
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS  --script fuzzy
+                        --files ${CMAKE_SOURCE_DIR}/test/references/steaminjectionbox-reference.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_boxsteaminjection-00007.vtu
+                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p1cni_box test_boxsteaminjection.input")
 
-add_dumux_test(test_ccsteaminjection test_ccsteaminjection test_ccsteaminjection.cc
-               python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
-                 --script fuzzy
-                 --files ${CMAKE_SOURCE_DIR}/test/references/steaminjectioncc-reference.vtu
-                         ${CMAKE_CURRENT_BINARY_DIR}/test_ccsteaminjection-00009.vtu
-                 --command "${CMAKE_CURRENT_BINARY_DIR}/test_ccsteaminjection")
+dune_add_test(NAME test_2p1cni_tpfa
+              SOURCES test_2p1c_fv.cc
+              COMPILE_DEFINITIONS TYPETAG=TwoPOneCNICCTpfaTypeTag
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS  --script fuzzy
+                        --files ${CMAKE_SOURCE_DIR}/test/references/steaminjectioncc-reference.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_ccsteaminjection-00009.vtu
+                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p1cni_tpfa test_ccsteaminjection.input")
diff --git a/test/porousmediumflow/2p1c/implicit/steaminjectionproblem.hh b/test/porousmediumflow/2p1c/implicit/steaminjectionproblem.hh
index 97a3347bd0eba5a200f2ec79e9a2983240d388ad..9af5bc1b872fb80531789df908082d89095c1584 100644
--- a/test/porousmediumflow/2p1c/implicit/steaminjectionproblem.hh
+++ b/test/porousmediumflow/2p1c/implicit/steaminjectionproblem.hh
@@ -25,13 +25,18 @@
 #ifndef DUMUX_STEAM_INJECTIONPROBLEM_HH
 #define DUMUX_STEAM_INJECTIONPROBLEM_HH
 
-#include <dumux/porousmediumflow/2p1c/implicit/model.hh>
-#include <dumux/porousmediumflow/implicit/problem.hh>
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/discretization/box/properties.hh>
+#include <dumux/porousmediumflow/2p1c/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include <dumux/material/fluidsystems/2pliquidvapor.hh>
 
-#include "steaminjectionspatialparams.hh"
 #include <dumux/material/components/tabulatedcomponent.hh>
 #include <dumux/material/components/h2o.hh>
 
+#include "steaminjectionspatialparams.hh"
+
 namespace Dumux
 {
 template <class TypeTag>
@@ -39,55 +44,59 @@ class InjectionProblem;
 
 namespace Properties
 {
-NEW_TYPE_TAG(InjectionProblem, INHERITS_FROM(TwoPOneCNI, InjectionProblemSpatialParams));
-NEW_TYPE_TAG(InjectionBoxProblem, INHERITS_FROM(BoxModel, InjectionProblem));
-NEW_TYPE_TAG(InjectionCCProblem, INHERITS_FROM(CCModel, InjectionProblem));
+NEW_TYPE_TAG(InjectionProblemTypeTag, INHERITS_FROM(TwoPOneCNI, InjectionProblemSpatialParams));
+NEW_TYPE_TAG(TwoPOneCNIBoxTypeTag, INHERITS_FROM(BoxModel, InjectionProblemTypeTag));
+NEW_TYPE_TAG(TwoPOneCNICCTpfaTypeTag, INHERITS_FROM(CCTpfaModel, InjectionProblemTypeTag));
 
-SET_TYPE_PROP(InjectionProblem, Grid, Dune::YaspGrid<2>);
+SET_TYPE_PROP(InjectionProblemTypeTag, Grid, Dune::YaspGrid<2>);
 
 // Set the problem property
-SET_PROP(InjectionProblem, Problem)
-{
-    typedef Dumux::InjectionProblem<TypeTag> type;
-};
+SET_TYPE_PROP(InjectionProblemTypeTag, Problem, InjectionProblem<TypeTag>);
 
-// Set the fluid system
-SET_PROP(InjectionProblem, FluidSystem)
+
+// Set fluid configuration
+SET_PROP(InjectionProblemTypeTag, FluidSystem)
 {
 private:
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar> >  H2OType ;
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using H2OType = Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar> >;
 public:
-    typedef Dumux::FluidSystems::TwoPLiquidVaporFluidsystem<Scalar, H2OType > type;
+    using type = Dumux::FluidSystems::TwoPLiquidVaporFluidsystem<Scalar, H2OType >;
 };
 
-// Enable gravity
-SET_BOOL_PROP(InjectionProblem, ProblemEnableGravity, true);
 
 //Define whether spurious cold-water flow into the steam is blocked
-SET_BOOL_PROP(InjectionProblem, UseBlockingOfSpuriousFlow, true);
+SET_BOOL_PROP(InjectionProblemTypeTag, UseBlockingOfSpuriousFlow, true);
 }
 
-//TODO: Names
 /*!
- * \ingroup ThreePTwoCNIBoxModel
- * \ingroup ImplicitTestProblems
+ * \ingroup TwoPOneC
  * \brief Non-isothermal 2D problem where steam is injected on the lower left side of the domain.
  *
- * This problem uses the \ref ThreePTwoCNIModel.
+ * This problem uses the \ref TwoPOneC model.
  *
  *  */
 template <class TypeTag>
-class InjectionProblem : public ImplicitPorousMediaProblem<TypeTag>
+class InjectionProblem : public PorousMediumFlowProblem<TypeTag>
 {
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    typedef typename GridView::Grid Grid;
-
-    typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
+    using ParentType = PorousMediumFlowProblem<TypeTag>;
+
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     // copy some indices for convenience
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
     enum {
         pressureIdx = Indices::pressureIdx,
         switch1Idx = Indices::switch1Idx,
@@ -96,51 +105,27 @@ class InjectionProblem : public ImplicitPorousMediaProblem<TypeTag>
 
         // phase and component indices
         wPhaseIdx = Indices::wPhaseIdx,
-        gPhaseIdx = Indices::gPhaseIdx,
+        nPhaseIdx = Indices::nPhaseIdx,
 
         // Phase State
         wPhaseOnly = Indices::wPhaseOnly,
-        gPhaseOnly = Indices::gPhaseOnly,
+        nPhaseOnly = Indices::nPhaseOnly,
         twoPhases = Indices::twoPhases,
-
-        // Grid and world dimension
-        dim = GridView::dimension,
-        dimWorld = GridView::dimensionworld
     };
 
-    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
-    typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
-
-    typedef typename GridView::template Codim<0>::Entity Element;
-    typedef typename GridView::template Codim<dim>::Entity Vertex;
-
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-
-    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
-
-    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+    static constexpr int dim = GridView::dimension;
+    static constexpr int dimWorld = GridView::dimensionworld;
+    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
 public:
 
     /*!
      * \brief The constructor
      *
-     * \param timeManager The time manager
-     * \param gridView The grid view
      */
-
-    InjectionProblem(TimeManager &timeManager, const GridView &gridView)
-        : ParentType(timeManager, gridView)
+    InjectionProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry)
     {
-
-        name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,
-                                             std::string,
-                                             Problem,
-                                             Name);
-
         FluidSystem::init();
     }
 
@@ -149,20 +134,15 @@ public:
      */
     // \{
 
-    /*!
-     * \brief The problem name.
-     *
-     * This is used as a prefix for files generated by the simulation.
-     */
-    const std::string name() const
-    { return name_; }
-
 
-    void sourceAtPos(PrimaryVariables &values,
-                     const GlobalPosition &globalPos) const
-     {
-          values = 0.0;
-     }
+    //! \copydoc Dumux::FVProblem::source()
+    Sources source(const Element &element,
+                   const FVElementGeometry& fvGeometry,
+                   const ElementVolumeVariables& elemVolVars,
+                   const SubControlVolume &scv) const
+    {
+          return Sources(0.0);
+    }
 
     /*!
      * \name Boundary conditions
@@ -171,58 +151,71 @@ public:
 
     /*!
      * \brief Specifies which kind of boundary condition should be
-     *        used for which equation on a given boundary segment.
+     *        used for which equation on a given boundary segment
      *
-     * \param bcTypes The boundary types for the conservation equations
-     * \param globalPos The position for which the bc type should be evaluated
+     * \param globalPos The global position
      */
-   void boundaryTypesAtPos(BoundaryTypes &bcTypes,
-            const GlobalPosition &globalPos) const
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
     {
-        if(globalPos[1] > this->bBoxMax()[1] - eps_ || globalPos[0] > this->bBoxMax()[0] - eps_)
+        BoundaryTypes bcTypes;
+
+        if(globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_ || globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_)
            bcTypes.setAllDirichlet();
         else
            bcTypes.setAllNeumann();
+
+         return bcTypes;
     }
 
     /*!
-     * \brief Evaluate the boundary conditions for a dirichlet
-     *        boundary segment.
-     *
-     * \param values The dirichlet values for the primary variables
-     * \param globalPos The position for which the bc type should be evaluated
+     * \brief Evaluates the boundary conditions for a Dirichlet
+     *        boundary segment
      *
-     * For this method, the \a values parameter stores primary variables.
+     * \param globalPos The global position
      */
-    void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+    PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
     {
-       initial_(values, globalPos);
+       return initialAtPos(globalPos);
     }
 
     /*!
      * \brief Evaluate the boundary conditions for a neumann
      *        boundary segment.
      *
-     * \param values The dirichlet values for the primary variables
-     * \param globalPos The position for which the bc type should be evaluated
+     * This is the method for the case where the Neumann condition is
+     * potentially solution dependent and requires some quantities that
+     * are specific to the fully-implicit method.
      *
-     * For this method, the \a values parameter stores the mass flux
+     * \param values The neumann values for the conservation equations in units of
+     *                 \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$
+     * \param element The finite element
+     * \param fvGeometry The finite-volume geometry
+     * \param elemVolVars All volume variables for the element
+     * \param scvf The sub control volume face
+     *
+     * For this method, the \a values parameter stores the flux
      * in normal direction of each phase. Negative values mean influx.
+     * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
      */
-      void neumannAtPos(PrimaryVariables &values,
-                      const GlobalPosition &globalPos) const
+    ResidualVector neumann(const Element& element,
+                           const FVElementGeometry& fvGeometry,
+                           const ElementVolumeVariables& elemVolVars,
+                           const SubControlVolumeFace& scvf) const
     {
-        values = 0.0;
+        ResidualVector values(0.0);
+
+        const auto& ipGlobal = scvf.ipGlobal();
 
-        if (globalPos[0] < eps_)
+        if (ipGlobal[0] < eps_)
         {
-            if(globalPos[1] > 2.0 - eps_ && globalPos[1] < 3.0 + eps_)
+            if(ipGlobal[1] > 2.0 - eps_ && ipGlobal[1] < 3.0 + eps_)
             {
-            Scalar massRate = 1e-1;
-            values[Indices::conti0EqIdx] = -massRate;
-            values[Indices::energyEqIdx] = -massRate * 2690e3;
+                const Scalar massRate = 1e-1;
+                values[Indices::conti0EqIdx] = -massRate;
+                values[Indices::energyEqIdx] = -massRate * 2690e3;
             }
         }
+        return values;
     }
 
     // \}
@@ -233,48 +226,26 @@ public:
     // \{
 
     /*!
-     * \brief Evaluate the initial value for a control volume.
-     *
-     * \param values The initial values for the primary variables
-     * \param globalPos The position for which the initial condition should be evaluated
-     *
-     * For this method, the \a values parameter stores primary
-     * variables.
-     */
-    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
-    {
-        initial_(values, globalPos);
-    }
-
-    /*!
-     * \brief Return the initial phase state inside a control volume.
+     * \brief Evaluates the initial values for a control volume
      *
-     * \param vert The vertex
-     * \param globalIdx The index of the global vertex
      * \param globalPos The global position
      */
-    int initialPhasePresence(const Vertex &vert,
-                             int &globalIdx,
-                             const GlobalPosition &globalPos) const
+    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
     {
-        return wPhaseOnly;
-    }
+        PrimaryVariables values(0.0);
 
+        const Scalar densityW = 1000.0;
+        values[pressureIdx] = 101300.0 + (this->fvGridGeometry().bBoxMax()[1] - globalPos[1])*densityW*9.81; // hydrostatic pressure
+        values[switch1Idx] = 283.13;
 
+        values.setState(wPhaseOnly);
 
-private:
-    // internal method for the initial condition (reused for the
-    // dirichlet conditions!)
-    void initial_(PrimaryVariables &values,
-                  const GlobalPosition &globalPos) const
-    {
-        Scalar densityW = 1000.0;
-        values[pressureIdx] = 101300.0 + (this->bBoxMax()[1] - globalPos[1])*densityW*9.81; // hydrostatic pressure
-        values[switch1Idx] = 283.13;
+        return values;
     }
 
+private:
+
     static constexpr Scalar eps_ = 1e-6;
-    std::string name_;
 };
 } //end namespace
 
diff --git a/test/porousmediumflow/2p1c/implicit/steaminjectionspatialparams.hh b/test/porousmediumflow/2p1c/implicit/steaminjectionspatialparams.hh
index 43f02fc690c822b56ef00a58b59f5808f7a0a756..48a00ca88b29986d2095d11a6d8cc7ddcc25478c 100644
--- a/test/porousmediumflow/2p1c/implicit/steaminjectionspatialparams.hh
+++ b/test/porousmediumflow/2p1c/implicit/steaminjectionspatialparams.hh
@@ -64,34 +64,36 @@ SET_PROP(InjectionProblemSpatialParams, MaterialLaw)
 template<class TypeTag>
 class InjectionProblemSpatialParams : public ImplicitSpatialParams<TypeTag>
 {
-    typedef ImplicitSpatialParams<TypeTag> ParentType;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename Grid::ctype CoordScalar;
-    enum {
-        dim=GridView::dimension,
-        dimWorld=GridView::dimensionworld
-    };
-
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-    typedef typename GridView::template Codim<0>::Entity Element;
-
-    typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix;
-
-
+    using ParentType = ImplicitSpatialParams<TypeTag>;
+
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+
+    using MaterialLawParams = typename MaterialLaw::Params;
+    using CoordScalar = typename GridView::ctype;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+
+    static constexpr int dim = GridView::dimension;
+    static constexpr int dimWorld = GridView::dimensionworld;
+    static constexpr int wPhaseIdx = FluidSystem::wPhaseIdx;
+
+    using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
+    using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
 public:
-    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
-    typedef typename MaterialLaw::Params MaterialLawParams;
+    using PermeabilityType = DimWorldMatrix;
 
     /*!
      * \brief The constructor
      *
      * \param gridView The grid view
      */
-    InjectionProblemSpatialParams(const GridView &gridView)
-        : ParentType(gridView)
+    InjectionProblemSpatialParams(const Problem& problem)
+    : ParentType(problem)
     {
         // set Van Genuchten Parameters
         materialParams_.setSwr(0.1);
@@ -101,18 +103,13 @@ public:
     }
 
     /*!
-     * \brief Apply the intrinsic permeability tensor to a pressure
-     *        potential gradient.
+     * \brief Returns the hydraulic conductivity \f$[m^2]\f$
      *
-     * \param element The current finite element
-     * \param fvElemGeom The current finite volume geometry of the element
-     * \param scvIdx The index of the sub-control volume
+     * \param globalPos The global position
      */
-    const FieldMatrix intrinsicPermeability(const Element &element,
-                                       const FVElementGeometry &fvElemGeom,
-                                       int scvIdx) const
+    DimWorldMatrix permeabilityAtPos(const GlobalPosition& globalPos) const
     {
-        FieldMatrix permMatrix(0.0);
+        DimWorldMatrix permMatrix(0.0);
 
         // intrinsic permeability
         permMatrix[0][0] = 1e-9;
@@ -124,78 +121,62 @@ public:
     /*!
      * \brief Define the porosity \f$[-]\f$ of the spatial parameters
      *
-     * \param element The finite element
-     * \param fvGeometry The finite volume geometry
-     * \param scvIdx The local index of the sub-control volume where
-     *                    the porosity needs to be defined
+     * \param globalPos The global position
      */
-    double porosity(const Element &element,
-                    const FVElementGeometry &fvGeometry,
-                    const int scvIdx) const
+    Scalar porosityAtPos(const GlobalPosition& globalPos) const
     {
         return 0.4;
     }
 
-
     /*!
-     * \brief return the parameter object for the Brooks-Corey material law which depends on the position
+     * \brief Returns the parameter object for the capillary-pressure/
+     *        saturation material law
      *
-     * \param element The current finite element
-     * \param fvGeometry The current finite volume geometry of the element
-     * \param scvIdx The index of the sub-control volume
+     * \param globalPos The global position
      */
-    const MaterialLawParams& materialLawParams(const Element &element,
-                                               const FVElementGeometry &fvGeometry,
-                                               const int scvIdx) const
+    const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
     {
         return materialParams_;
     }
-
     /*!
      * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
      *
      * This is only required for non-isothermal models.
      *
-     * \param element The finite element
-     * \param fvGeometry The finite volume geometry
-     * \param scvIdx The local index of the sub-control volume
+     * \param element The element
+     * \param scv The sub control volume
+     * \param elemSol The element solution vector
      */
-     Scalar solidHeatCapacity(const Element &element,
-                             const FVElementGeometry &fvGeometry,
-                             const int scvIdx) const
-    {
-        return 850.0; // specific heat capacity of granite [J / (kg K)]
-    }
+    Scalar solidHeatCapacity(const Element &element,
+                             const SubControlVolume& scv,
+                             const ElementSolutionVector& elemSol) const
+    { return 850.0; /*specific heat capacity of granite [J / (kg K)]*/ }
 
     /*!
      * \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix.
      *
      * This is only required for non-isothermal models.
      *
-     * \param element The finite element
-     * \param fvGeometry The finite volume geometry
-     * \param scvIdx The local index of the sub-control volume
+     * \param element The element
+     * \param scv The sub control volume
+     * \param elemSol The element solution vector
      */
     Scalar solidDensity(const Element &element,
-                        const FVElementGeometry &fvGeometry,
-                        const int scvIdx) const
-    {
-        return 2650.0; // density of granite [kg/m^3]
-    }
+                        const SubControlVolume& scv,
+                        const ElementSolutionVector& elemSol) const
+    { return 2650; /*density of granite [kg/m^3]*/ }
 
     /*!
-     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
      *
-     * \param element The finite element
-     * \param fvGeometry The finite volume geometry
-     * \param scvIdx The local index of the sub-control volume
+     * \param element The element
+     * \param scv The sub control volume
+     * \param elemSol The element solution vector
      */
     Scalar solidThermalConductivity(const Element &element,
-                                    const FVElementGeometry &fvGeometry,
-                                    const int scvIdx) const
-    {
-        return 2.8;
-    }
+                                    const SubControlVolume& scv,
+                                    const ElementSolutionVector& elemSol) const
+    { return 2.8; }
 
 private:
     MaterialLawParams materialParams_;
diff --git a/test/porousmediumflow/2p1c/implicit/test_2p1c_fv.cc b/test/porousmediumflow/2p1c/implicit/test_2p1c_fv.cc
new file mode 100644
index 0000000000000000000000000000000000000000..94aafde415fa6cbe7771fdff0bf8dd5ad5607499
--- /dev/null
+++ b/test/porousmediumflow/2p1c/implicit/test_2p1c_fv.cc
@@ -0,0 +1,218 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program 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 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief test for the 1pnc model
+ */
+ #include <config.h>
+
+ #include "steaminjectionproblem.hh"
+
+ #include <ctime>
+ #include <iostream>
+
+ #include <dune/common/parallel/mpihelper.hh>
+ #include <dune/common/timer.hh>
+ #include <dune/grid/io/file/dgfparser/dgfexception.hh>
+ #include <dune/grid/io/file/vtk.hh>
+ #include <dune/istl/io.hh>
+
+ #include <dumux/discretization/methods.hh>
+
+ #include <dumux/common/properties.hh>
+ #include <dumux/common/parameters.hh>
+ #include <dumux/common/valgrind.hh>
+ #include <dumux/common/dumuxmessage.hh>
+ #include <dumux/common/defaultusagemessage.hh>
+ #include <dumux/common/parameterparser.hh>
+
+ #include <dumux/porousmediumflow/compositional/privarswitchnewtoncontroller.hh>
+ #include <dumux/linear/seqsolverbackend.hh>
+ #include <dumux/nonlinear/newtonmethod.hh>
+
+ #include <dumux/assembly/fvassembler.hh>
+
+ #include <dumux/io/vtkoutputmodule.hh>
+
+ int main(int argc, char** argv) try
+ {
+     using namespace Dumux;
+
+     // define the type tag for this problem
+     using TypeTag = TTAG(TYPETAG);
+
+     ////////////////////////////////////////////////////////////
+     ////////////////////////////////////////////////////////////
+
+     // initialize MPI, finalize is done automatically on exit
+     const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+     // print dumux start message
+     if (mpiHelper.rank() == 0)
+         DumuxMessage::print(/*firstCall=*/true);
+
+     // initialize parameter tree
+     Parameters::init(argc, argv);
+
+     //////////////////////////////////////////////////////////////////////
+     // try to create a grid (from the given grid file or the input file)
+     /////////////////////////////////////////////////////////////////////
+
+     using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+     GridCreator::makeGrid();
+     GridCreator::loadBalance();
+
+     ////////////////////////////////////////////////////////////
+     // run instationary non-linear problem on this grid
+     ////////////////////////////////////////////////////////////
+
+     // we compute on the leaf grid view
+     const auto& leafGridView = GridCreator::grid().leafGridView();
+
+     // create the finite volume grid geometry
+     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+     auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+     fvGridGeometry->update();
+
+     // the problem (initial and boundary conditions)
+     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+     auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+     // the solution vector
+     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+     SolutionVector x(fvGridGeometry->numDofs());
+     problem->applyInitialSolution(x);
+     auto xOld = x;
+
+     // the grid variables
+     using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+     auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+     gridVariables->init(x, xOld);
+
+     // get some time loop parameters
+     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+     auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+     auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+     auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+     auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+
+     // intialize the vtk output module
+     VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+     using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+     VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+     vtkWriter.write(0.0);
+
+     // instantiate time loop
+     auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
+     timeLoop->setMaxTimeStepSize(maxDt);
+
+     // the assembler with time loop for instationary problem
+     using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+     auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+     // the linear solver
+     using LinearSolver = ILU0BiCGSTABBackend<TypeTag>;
+     auto linearSolver = std::make_shared<LinearSolver>();
+
+     // the non-linear solver
+     using NewtonController = Dumux::PriVarSwitchNewtonController<TypeTag>;
+     auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+     NewtonMethod<NewtonController, Assembler, LinearSolver> nonLinearSolver(newtonController, assembler, linearSolver);
+
+     // time loop
+     timeLoop->start(); do
+     {
+         // set previous solution for storage evaluations
+         assembler->setPreviousSolution(xOld);
+
+         // try solving the non-linear system
+         for (int i = 0; i < maxDivisions; ++i)
+         {
+             // linearize & solve
+             auto converged = nonLinearSolver.solve(x);
+
+             if (converged)
+                 break;
+
+             if (!converged && i == maxDivisions-1)
+                 DUNE_THROW(Dune::MathError,
+                            "Newton solver didn't converge after "
+                            << maxDivisions
+                            << " time-step divisions. dt="
+                            << timeLoop->timeStepSize()
+                            << ".\nThe solutions of the current and the previous time steps "
+                            << "have been saved to restart files.");
+         }
+
+         // make the new solution the old solution
+         xOld = x;
+         gridVariables->advanceTimeStep();
+
+         // advance to the time loop to the next step
+         timeLoop->advanceTimeStep();
+
+         // write vtk output
+         vtkWriter.write(timeLoop->time());
+
+         // report statistics of this time step
+         timeLoop->reportTimeStep();
+
+         // set new dt as suggested by newton controller
+         timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+
+     } while (!timeLoop->finished());
+
+     timeLoop->finalize(leafGridView.comm());
+
+     ////////////////////////////////////////////////////////////
+     // finalize, print dumux message to say goodbye
+     ////////////////////////////////////////////////////////////
+
+     // print dumux end message
+     if (mpiHelper.rank() == 0)
+         DumuxMessage::print(/*firstCall=*/false);
+
+     return 0;
+
+ }
+ catch (Dumux::ParameterException &e)
+ {
+     std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+     return 1;
+ }
+ catch (Dune::DGFException & e)
+ {
+     std::cerr << "DGF exception thrown (" << e <<
+                  "). Most likely, the DGF file name is wrong "
+                  "or the DGF file is corrupted, "
+                  "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+                  << " ---> Abort!" << std::endl;
+     return 2;
+ }
+ catch (Dune::Exception &e)
+ {
+     std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+     return 3;
+ }
+ catch (...)
+ {
+     std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+     return 4;
+ }
diff --git a/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.cc b/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.cc
deleted file mode 100644
index bd83890846bb39f76a3e9233e83fa289156e7bfe..0000000000000000000000000000000000000000
--- a/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.cc
+++ /dev/null
@@ -1,62 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program 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 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Test for the 2p1cni box model
- */
-#include "config.h"
-#include "steaminjectionproblem.hh"
-#include <dumux/common/start.hh>
-
-/*!
- * \brief Provides an interface for customizing error messages associated with
- *        reading in parameters.
- *
- * \param progName  The name of the program, that was tried to be started.
- * \param errorMsg  The error message that was issued by the start function.
- *                  Comprises the thing that went wrong and a general help message.
- */
-void usage(const char *progName, const std::string &errorMsg)
-{
-    if (errorMsg.size() > 0) {
-        std::string errorMessageOut = "\nUsage: ";
-        errorMessageOut += progName;
-        errorMessageOut += " [options]\n";
-        errorMessageOut += errorMsg;
-        errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
-                           "\t-TimeManager.TEnd              End of the simulation [s] \n"
-                           "\t-TimeManager.DtInitial         Initial timestep size [s] \n"
-                           "\t-Grid.UpperRight               coordinates of the upper right corner of the grid [m] \n"
-                           "\t-Grid.Cells                    Number of cells in respective coordinate directions\n"
-                           "\t-Problem.Name                  String for naming of the output files \n"
-                           "\n";
-
-        std::cout << errorMessageOut << std::endl;
-    }
-}
-
-////////////////////////
-// the main function
-////////////////////////
-int main(int argc, char** argv)
-{
-    typedef TTAG(InjectionBoxProblem) TypeTag;
-    return Dumux::start<TypeTag>(argc, argv, usage);
-}
diff --git a/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.input b/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.input
index 4696345154e120977b7f4c2d8377e28ec4cdfc07..59aa81a23af877876c9ffba67dd3c141602d507b 100644
--- a/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.input
+++ b/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.input
@@ -1,4 +1,4 @@
-[TimeManager]
+[TimeLoop]
 DtInitial = 50 # [s]
 TEnd = 300 # [s]
 
diff --git a/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.cc b/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.cc
deleted file mode 100644
index 1273bd90ede4491753202a2f9cfef6b72e669adc..0000000000000000000000000000000000000000
--- a/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.cc
+++ /dev/null
@@ -1,62 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program 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 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Test for the 2p1cni cell-centered model
- */
-#include "config.h"
-#include "steaminjectionproblem.hh"
-#include <dumux/common/start.hh>
-
-/*!
- * \brief Provides an interface for customizing error messages associated with
- *        reading in parameters.
- *
- * \param progName  The name of the program, that was tried to be started.
- * \param errorMsg  The error message that was issued by the start function.
- *                  Comprises the thing that went wrong and a general help message.
- */
-void usage(const char *progName, const std::string &errorMsg)
-{
-    if (errorMsg.size() > 0) {
-        std::string errorMessageOut = "\nUsage: ";
-        errorMessageOut += progName;
-        errorMessageOut += " [options]\n";
-        errorMessageOut += errorMsg;
-        errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
-                           "\t-TimeManager.TEnd              End of the simulation [s] \n"
-                           "\t-TimeManager.DtInitial         Initial timestep size [s] \n"
-                           "\t-Grid.UpperRight               coordinates of the upper right corner of the grid [m] \n"
-                           "\t-Grid.Cells                    Number of cells in respective coordinate directions\n"
-                           "\t-Problem.Name                  String for naming of the output files \n"
-                           "\n";
-
-        std::cout << errorMessageOut << std::endl;
-    }
-}
-
-////////////////////////
-// the main function
-////////////////////////
-int main(int argc, char** argv)
-{
-    typedef TTAG(InjectionCCProblem) TypeTag;
-    return Dumux::start<TypeTag>(argc, argv, usage);
-}
diff --git a/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.input b/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.input
index c5b8d8556b1fd4b9c149baa9212b1684d3db858d..fabf21bed2ff303018d5c1d26fc55ec3f9f00e87 100644
--- a/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.input
+++ b/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.input
@@ -1,4 +1,4 @@
-[TimeManager]
+[TimeLoop]
 DtInitial = 50 # [s]
 TEnd = 300 # [s]
 
diff --git a/test/references/steaminjectionbox-reference.vtu b/test/references/steaminjectionbox-reference.vtu
index 138f79f6d4849ce5604fd82f9179b3315bcf566b..feca2fc2a7f2e369f1266dd2a48269993a8dfc4f 100644
--- a/test/references/steaminjectionbox-reference.vtu
+++ b/test/references/steaminjectionbox-reference.vtu
@@ -2,15 +2,15 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="200" NumberOfPoints="231">
-      <PointData Scalars="sw">
-        <DataArray type="Float32" Name="sw" NumberOfComponents="1" format="ascii">
+      <PointData Scalars="Sw">
+        <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
           1 1 1 1 1 1 1 1 1 1 1 1
           1 1 1 1 1 1 1 1 1 1 1 1
           1 1 1 1 1 1 1 1 1 1 1 1
           1 1 1 1 1 1 1 1 1 1 1 1
           1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 0.822808 1 1 1 1 1
-          1 1 1 1 1 0.74307 1 1 1 1 1 1
+          1 1 1 1 1 1 0.79171 1 1 1 1 1
+          1 1 1 1 1 0.71017 1 1 1 1 1 1
           1 1 1 1 1 1 1 1 1 1 1 1
           1 1 1 1 1 1 1 1 1 1 1 1
           1 1 1 1 1 1 1 1 1 1 1 1
@@ -25,14 +25,14 @@
           1 1 1 1 1 1 1 1 1 1 1 1
           1 1 1
         </DataArray>
-        <DataArray type="Float32" Name="sg" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0.177192 0 0 0 0 0
-          0 0 0 0 0 0.25693 0 0 0 0 0 0
+          0 0 0 0 0 0 0.20829 0 0 0 0 0
+          0 0 0 0 0 0.28983 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
@@ -48,123 +48,145 @@
           0 0 0
         </DataArray>
         <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
-          180140 180133 176223 176217 180115 176198 180087 176168 180052 176131 180011 176090
-          179967 176045 179921 175998 179874 175951 179826 175903 179780 175856 172320 172312
-          172290 172257 172216 172171 172124 172076 172028 171980 171932 168434 168423 168394
-          168352 168304 168254 168204 168155 168105 168057 168008 164568 164554 164511 164452
-          164394 164338 164284 164233 164182 164133 164084 160692 160744 160632 160550 160480
-          160419 160362 160309 160258 160209 160160 157303 156850 156733 156637 156560 156495
-          156437 156384 156333 156284 156236 153356 152873 152801 152705 152628 152564 152507
-          152455 152406 152359 152312 148899 148944 148833 148753 148686 148626 148573 148524
-          148477 148432 148388 144946 144901 144849 144790 144734 144683 144635 144590 144547
-          144505 144464 140913 140897 140863 140822 140779 140736 140695 140655 140616 140578
-          140540 136914 136906 136885 136856 136823 136788 136753 136718 136684 136650 136616
-          132933 132928 132914 132893 132868 132840 132811 132781 132751 132722 132692 128964
-          128960 128950 128934 128915 128893 128869 128844 128819 128793 128768 125001 124999
-          124991 124980 124964 124947 124928 124908 124887 124865 124844 121044 121043 121037
-          121028 121017 121003 120988 120972 120955 120937 120920 117091 117090 117086 117079
-          117071 117060 117049 117036 117023 117010 116996 113141 113140 113137 113133 113127
-          113119 113111 113102 113092 113082 113072 109193 109192 109191 109188 109184 109179
-          109173 109167 109161 109155 109148 105246 105246 105245 105244 105242 105239 105237
-          105234 105231 105227 105224 101300 101300 101300 101300 101300 101300 101300 101300
+          179934 179932 176014 176011 179924 176003 179913 175992 179898 175976 179880 175958
+          179861 175939 179841 175919 179821 175898 179800 175877 179780 175856 172098 172096
+          172087 172074 172057 172038 172017 171996 171975 171953 171932 168189 168186 168176
+          168160 168140 168118 168096 168074 168052 168030 168008 164284 164285 164271 164248
+          164224 164200 164176 164152 164129 164107 164084 160395 160410 160371 160338 160308
+          160280 160254 160230 160206 160183 160160 156677 156518 156468 156426 156390 156359
+          156332 156306 156282 156259 156236 152955 152588 152555 152506 152467 152435 152407
+          152382 152358 152335 152312 148678 148703 148622 148575 148538 148507 148480 148455
+          148432 148410 148388 144739 144708 144672 144635 144604 144576 144551 144528 144506
+          144485 144464 140747 140737 140715 140690 140665 140642 140620 140599 140579 140559
+          140540 136780 136775 136762 136744 136725 136706 136687 136669 136651 136633 136616
+          132824 132821 132812 132800 132785 132770 132754 132738 132723 132707 132692 128874
+          128872 128866 128857 128846 128834 128821 128808 128795 128781 128768 124929 124928
+          124923 124917 124908 124899 124888 124877 124866 124855 124844 120987 120986 120983
+          120978 120972 120964 120956 120947 120938 120929 120920 117047 117047 117044 117041
+          117036 117030 117024 117017 117010 117003 116996 113109 113109 113107 113105 113101
+          113097 113093 113088 113083 113077 113072 109172 109172 109171 109169 109167 109164
+          109162 109158 109155 109152 109148 105236 105236 105235 105234 105233 105232 105231
+          105229 105227 105226 105224 101300 101300 101300 101300 101300 101300 101300 101300
           101300 101300 101300
         </DataArray>
-        <DataArray type="Float32" Name="pg" NumberOfComponents="1" format="ascii">
-          180140 180133 176223 176217 180115 176198 180087 176168 180052 176131 180011 176090
-          179967 176045 179921 175998 179874 175951 179826 175903 179780 175856 172320 172312
-          172290 172257 172216 172171 172124 172076 172028 171980 171932 168434 168423 168394
-          168352 168304 168254 168204 168155 168105 168057 168008 164568 164554 164511 164452
-          164394 164338 164284 164233 164182 164133 164084 160692 160744 160632 160550 160480
-          160419 160362 160309 160258 160209 160160 157568 156850 156733 156637 156560 156495
-          156437 156384 156333 156284 156236 153705 152873 152801 152705 152628 152564 152507
-          152455 152406 152359 152312 148899 148944 148833 148753 148686 148626 148573 148524
-          148477 148432 148388 144946 144901 144849 144790 144734 144683 144635 144590 144547
-          144505 144464 140913 140897 140863 140822 140779 140736 140695 140655 140616 140578
-          140540 136914 136906 136885 136856 136823 136788 136753 136718 136684 136650 136616
-          132933 132928 132914 132893 132868 132840 132811 132781 132751 132722 132692 128964
-          128960 128950 128934 128915 128893 128869 128844 128819 128793 128768 125001 124999
-          124991 124980 124964 124947 124928 124908 124887 124865 124844 121044 121043 121037
-          121028 121017 121003 120988 120972 120955 120937 120920 117091 117090 117086 117079
-          117071 117060 117049 117036 117023 117010 116996 113141 113140 113137 113133 113127
-          113119 113111 113102 113092 113082 113072 109193 109192 109191 109188 109184 109179
-          109173 109167 109161 109155 109148 105246 105246 105245 105244 105242 105239 105237
-          105234 105231 105227 105224 101300 101300 101300 101300 101300 101300 101300 101300
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+          179934 179932 176014 176011 179924 176003 179913 175992 179898 175976 179880 175958
+          179861 175939 179841 175919 179821 175898 179800 175877 179780 175856 172098 172096
+          172087 172074 172057 172038 172017 171996 171975 171953 171932 168189 168186 168176
+          168160 168140 168118 168096 168074 168052 168030 168008 164284 164285 164271 164248
+          164224 164200 164176 164152 164129 164107 164084 160395 160410 160371 160338 160308
+          160280 160254 160230 160206 160183 160160 156974 156518 156468 156426 156390 156359
+          156332 156306 156282 156259 156236 153342 152588 152555 152506 152467 152435 152407
+          152382 152358 152335 152312 148678 148703 148622 148575 148538 148507 148480 148455
+          148432 148410 148388 144739 144708 144672 144635 144604 144576 144551 144528 144506
+          144485 144464 140747 140737 140715 140690 140665 140642 140620 140599 140579 140559
+          140540 136780 136775 136762 136744 136725 136706 136687 136669 136651 136633 136616
+          132824 132821 132812 132800 132785 132770 132754 132738 132723 132707 132692 128874
+          128872 128866 128857 128846 128834 128821 128808 128795 128781 128768 124929 124928
+          124923 124917 124908 124899 124888 124877 124866 124855 124844 120987 120986 120983
+          120978 120972 120964 120956 120947 120938 120929 120920 117047 117047 117044 117041
+          117036 117030 117024 117017 117010 117003 116996 113109 113109 113107 113105 113101
+          113097 113093 113088 113083 113077 113072 109172 109172 109171 109169 109167 109164
+          109162 109158 109155 109152 109148 105236 105236 105235 105234 105233 105232 105231
+          105229 105227 105226 105224 101300 101300 101300 101300 101300 101300 101300 101300
           101300 101300 101300
         </DataArray>
-        <DataArray type="Float32" Name="rhow" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 297.293 0 0 0 0 0
+          0 0 0 0 0 387.236 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
           999.729 999.729 999.727 999.727 999.729 999.727 999.729 999.727 999.729 999.727 999.729 999.727
           999.729 999.727 999.729 999.727 999.729 999.727 999.729 999.727 999.729 999.727 999.725 999.725
-          999.725 999.725 999.725 999.725 999.725 999.725 999.725 999.725 999.725 999.723 999.723 999.723
-          999.723 999.723 999.723 999.723 999.723 999.723 999.723 999.723 999.657 999.72 999.721 999.721
-          999.721 999.721 999.721 999.721 999.721 999.721 999.721 978.951 999.651 999.719 999.719 999.719
-          999.719 999.719 999.719 999.719 999.719 999.719 948.769 998.836 999.707 999.717 999.717 999.717
-          999.717 999.717 999.717 999.717 999.717 949.345 997.01 999.666 999.714 999.716 999.716 999.716
-          999.716 999.716 999.715 999.715 988.278 999.483 999.708 999.714 999.714 999.714 999.714 999.714
-          999.714 999.714 999.714 999.437 999.704 999.712 999.712 999.712 999.712 999.712 999.712 999.712
-          999.712 999.712 999.704 999.71 999.71 999.71 999.71 999.71 999.71 999.71 999.71 999.71
+          999.725 999.725 999.725 999.725 999.725 999.725 999.725 999.725 999.725 999.722 999.723 999.723
+          999.723 999.723 999.723 999.723 999.723 999.723 999.723 999.723 999.616 999.719 999.721 999.721
+          999.721 999.721 999.721 999.721 999.721 999.721 999.721 977.363 999.605 999.718 999.719 999.719
+          999.719 999.719 999.719 999.719 999.719 999.719 948.857 998.746 999.705 999.717 999.717 999.717
+          999.717 999.717 999.717 999.717 999.717 949.4 997.06 999.668 999.714 999.716 999.716 999.716
+          999.716 999.715 999.715 999.715 988.284 999.513 999.709 999.714 999.714 999.714 999.714 999.714
+          999.714 999.714 999.714 999.487 999.706 999.712 999.712 999.712 999.712 999.712 999.712 999.712
+          999.712 999.712 999.705 999.71 999.71 999.71 999.71 999.71 999.71 999.71 999.71 999.71
           999.71 999.708 999.708 999.708 999.708 999.708 999.708 999.708 999.708 999.708 999.708 999.708
           999.706 999.706 999.706 999.706 999.706 999.706 999.706 999.706 999.706 999.706 999.706 999.704
           999.704 999.704 999.704 999.704 999.704 999.704 999.704 999.704 999.704 999.704 999.702 999.702
           999.702 999.702 999.702 999.702 999.702 999.702 999.702 999.702 999.702 999.701 999.701 999.701
-          999.701 999.701 999.701 999.701 999.701 999.701 999.7 999.7 999.699 999.699 999.699 999.699
+          999.701 999.701 999.701 999.701 999.701 999.7 999.7 999.7 999.699 999.699 999.699 999.699
           999.699 999.699 999.699 999.699 999.699 999.699 999.699 999.697 999.697 999.697 999.697 999.697
           999.697 999.697 999.697 999.697 999.697 999.697 999.695 999.695 999.695 999.695 999.695 999.695
           999.695 999.695 999.695 999.695 999.695 999.693 999.693 999.693 999.693 999.693 999.693 999.693
           999.693 999.693 999.693 999.693 999.691 999.691 999.691 999.691 999.691 999.691 999.691 999.691
           999.691 999.691 999.691
         </DataArray>
-        <DataArray type="Float32" Name="rhog" NumberOfComponents="1" format="ascii">
-          1.38113 1.38108 1.3511 1.35105 1.38094 1.35091 1.38073 1.35068 1.38046 1.3504 1.38014 1.35008
-          1.3798 1.34973 1.37945 1.34938 1.37909 1.34901 1.37873 1.34865 1.37837 1.34829 1.32118 1.32112
-          1.32095 1.32069 1.32038 1.32003 1.31967 1.31931 1.31894 1.31857 1.3182 1.29136 1.2913 1.29108
-          1.29075 1.29039 1.29 1.28962 1.28924 1.28886 1.28849 1.28812 1.25832 1.26155 1.2613 1.26085
-          1.2604 1.25998 1.25957 1.25917 1.25878 1.25841 1.25803 1.03357 1.22892 1.23152 1.23093 1.2304
-          1.22993 1.2295 1.22909 1.2287 1.22832 1.22795 0.90318 1.17527 1.20107 1.20092 1.20034 1.19984
-          1.1994 1.19899 1.1986 1.19823 1.19786 0.882446 1.11391 1.16881 1.17072 1.1702 1.16971 1.16927
-          1.16887 1.16849 1.16813 1.16777 1.00739 1.1334 1.14078 1.14048 1.13997 1.13952 1.13911 1.13873
-          1.13837 1.13803 1.13769 1.10158 1.11054 1.11054 1.1101 1.10968 1.10928 1.10892 1.10857 1.10824
-          1.10792 1.1076 1.08005 1.08024 1.08 1.07968 1.07935 1.07902 1.0787 1.0784 1.0781 1.07781
-          1.07752 1.04971 1.04965 1.04949 1.04927 1.04902 1.04875 1.04848 1.04822 1.04795 1.04769 1.04743
-          1.0192 1.01916 1.01905 1.01889 1.01869 1.01848 1.01826 1.01803 1.0178 1.01757 1.01735 0.988761
-          0.988734 0.988655 0.988536 0.988386 0.988216 0.988034 0.987844 0.987651 0.987457 0.987261 0.958382 0.958362
-          0.958305 0.958214 0.958099 0.957965 0.957818 0.957662 0.957502 0.95734 0.957176 0.928044 0.928029 0.927987
-          0.927919 0.92783 0.927726 0.92761 0.927485 0.927356 0.927224 0.92709 0.897736 0.897725 0.897694 0.897644
-          0.897577 0.897498 0.897409 0.897313 0.897213 0.897109 0.897005 0.86745 0.867443 0.86742 0.867385 0.867337
-          0.867281 0.867216 0.867146 0.867072 0.866996 0.866919 0.837179 0.837175 0.83716 0.837138 0.837107 0.83707
-          0.837028 0.836983 0.836934 0.836884 0.836834 0.806919 0.806916 0.806909 0.806898 0.806883 0.806865 0.806845
-          0.806822 0.806798 0.806773 0.806748 0.776663 0.776663 0.776663 0.776663 0.776663 0.776663 0.776663 0.776663
+        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
+          1.37955 1.37953 1.3495 1.34948 1.37948 1.34942 1.37939 1.34933 1.37928 1.34921 1.37914 1.34907
+          1.379 1.34892 1.37884 1.34877 1.37868 1.34861 1.37853 1.34845 1.37837 1.34829 1.31948 1.31946
+          1.31939 1.31929 1.31916 1.31901 1.31886 1.31869 1.31853 1.31837 1.3182 1.28945 1.28948 1.2894
+          1.28928 1.28913 1.28896 1.28879 1.28862 1.28845 1.28828 1.28812 1.25465 1.25945 1.25946 1.25929
+          1.2591 1.25892 1.25873 1.25855 1.25838 1.2582 1.25803 1.0241 1.22476 1.22949 1.22931 1.22908
+          1.22887 1.22867 1.22848 1.2283 1.22812 1.22795 0.899997 1.17059 1.19893 1.1993 1.19904 1.1988
+          1.19859 1.1984 1.19821 1.19804 1.19786 0.880498 1.1125 1.16702 1.16919 1.16896 1.16872 1.1685
+          1.16831 1.16812 1.16795 1.16777 1.00593 1.13255 1.13922 1.13912 1.13884 1.1386 1.13839 1.1382
+          1.13803 1.13786 1.13769 1.10158 1.10915 1.10919 1.10892 1.10867 1.10846 1.10827 1.10809 1.10792
+          1.10776 1.1076 1.07885 1.07902 1.07886 1.07867 1.07848 1.0783 1.07813 1.07797 1.07781 1.07766
+          1.07752 1.04868 1.04865 1.04855 1.04842 1.04827 1.04812 1.04798 1.04784 1.0477 1.04757 1.04743
+          1.01836 1.01833 1.01827 1.01817 1.01806 1.01794 1.01782 1.0177 1.01758 1.01746 1.01735 0.988075
+          0.988059 0.988014 0.987945 0.987861 0.987767 0.987668 0.987567 0.987465 0.987363 0.987261 0.957828 0.957817
+          0.957784 0.957733 0.957669 0.957595 0.957515 0.957432 0.957347 0.957262 0.957176 0.927605 0.927597 0.927573
+          0.927535 0.927486 0.927429 0.927366 0.927299 0.927231 0.927161 0.92709 0.897399 0.897393 0.897375 0.897348
+          0.897311 0.897268 0.89722 0.897169 0.897115 0.89706 0.897005 0.867205 0.867201 0.867189 0.867169 0.867143
+          0.867112 0.867077 0.86704 0.867 0.86696 0.866919 0.83702 0.837017 0.837009 0.836997 0.83698 0.83696
+          0.836937 0.836913 0.836887 0.836861 0.836834 0.80684 0.806839 0.806835 0.806829 0.806821 0.806811 0.8068
+          0.806787 0.806775 0.806762 0.806748 0.776663 0.776663 0.776663 0.776663 0.776663 0.776663 0.776663 0.776663
           0.776663 0.776663 0.776663
         </DataArray>
-        <DataArray type="Float32" Name="MobW" NumberOfComponents="1" format="ascii">
-          764.381 764.381 764.379 764.379 764.381 764.379 764.381 764.379 764.381 764.379 764.381 764.379
-          764.381 764.379 764.381 764.379 764.381 764.379 764.381 764.378 764.381 764.378 764.377 764.376
-          764.376 764.376 764.376 764.376 764.376 764.376 764.376 764.376 764.376 764.517 764.379 764.374
-          764.374 764.374 764.374 764.374 764.374 764.374 764.374 764.374 782.422 764.828 764.378 764.372
-          764.372 764.372 764.372 764.372 764.372 764.372 764.372 2408.16 783.213 764.612 764.373 764.37
-          764.37 764.37 764.37 764.37 764.37 764.37 590.728 918.384 767.734 764.417 764.369 764.368
-          764.368 764.368 764.368 764.368 764.368 305.539 1128.42 780.017 764.747 764.372 764.366 764.366
-          764.366 764.366 764.366 764.366 1815.35 812.654 766.195 764.404 764.365 764.364 764.364 764.364
-          764.364 764.364 764.364 821.035 766.832 764.435 764.364 764.362 764.362 764.362 764.362 764.362
-          764.362 764.362 766.371 764.432 764.362 764.36 764.36 764.36 764.36 764.36 764.36 764.36
-          764.36 764.397 764.36 764.359 764.358 764.358 764.358 764.358 764.358 764.358 764.358 764.358
+        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+          764.381 764.381 764.378 764.378 764.381 764.378 764.381 764.378 764.381 764.378 764.381 764.378
+          764.381 764.378 764.381 764.378 764.381 764.378 764.381 764.378 764.381 764.378 764.378 764.376
+          764.376 764.376 764.376 764.376 764.376 764.376 764.376 764.376 764.376 764.642 764.382 764.374
+          764.374 764.374 764.374 764.374 764.374 764.374 764.374 764.374 789.74 765.006 764.381 764.372
+          764.372 764.372 764.372 764.372 764.372 764.372 764.372 2504.04 791.26 764.776 764.375 764.37
+          764.37 764.37 764.37 764.37 764.37 764.37 458.525 932.15 768.303 764.432 764.369 764.368
+          764.368 764.368 764.368 764.368 764.368 231.31 1123.86 779.537 764.749 764.372 764.366 764.366
+          764.366 764.366 764.366 764.366 1814.98 807.009 765.892 764.397 764.365 764.364 764.364 764.364
+          764.364 764.364 764.364 811.577 766.334 764.42 764.363 764.362 764.362 764.362 764.362 764.362
+          764.362 764.362 765.918 764.415 764.362 764.36 764.36 764.36 764.36 764.36 764.36 764.36
+          764.36 764.388 764.36 764.358 764.358 764.358 764.358 764.358 764.358 764.358 764.358 764.358
           764.357 764.356 764.356 764.356 764.356 764.356 764.356 764.356 764.356 764.356 764.355 764.354
           764.354 764.354 764.354 764.354 764.354 764.354 764.354 764.354 764.353 764.353 764.352 764.352
           764.352 764.352 764.352 764.352 764.352 764.352 764.351 764.351 764.351 764.35 764.35 764.35
-          764.35 764.35 764.35 764.349 764.349 764.349 764.349 764.349 764.348 764.348 764.348 764.347
+          764.35 764.35 764.349 764.349 764.349 764.349 764.349 764.349 764.348 764.348 764.348 764.347
           764.347 764.347 764.347 764.347 764.347 764.347 764.347 764.345 764.345 764.345 764.345 764.345
           764.345 764.345 764.345 764.345 764.345 764.345 764.343 764.343 764.343 764.343 764.343 764.343
           764.343 764.343 764.343 764.343 764.343 764.341 764.341 764.341 764.341 764.341 764.341 764.341
           764.341 764.341 764.341 764.341 764.339 764.339 764.339 764.339 764.339 764.339 764.339 764.339
           764.339 764.339 764.339
         </DataArray>
-        <DataArray type="Float32" Name="MobG" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 16289.6 0 0 0 0 0
-          0 0 0 0 0 25473.2 0 0 0 0 0 0
+          0 0 0 0 0 0 19827.9 0 0 0 0 0
+          0 0 0 0 0 29280.5 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
@@ -179,50 +201,6 @@
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0
         </DataArray>
-        <DataArray type="Float32" Name="ViscosW" NumberOfComponents="1" format="ascii">
-          0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825
-          0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130826
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130802 0.00130825 0.00130826
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00127808 0.00130748 0.00130825 0.00130826
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.000415254 0.00127679 0.00130785 0.00130826 0.00130827
-          0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.000247832 0.00108887 0.00130253 0.00130819 0.00130827 0.00130827
-          0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.000249633 0.000886195 0.00128202 0.00130762 0.00130826 0.00130827 0.00130827
-          0.00130827 0.00130827 0.00130827 0.00130827 0.000550859 0.00123054 0.00130515 0.00130821 0.00130828 0.00130828 0.00130828 0.00130828
-          0.00130828 0.00130828 0.00130828 0.00121797 0.00130407 0.00130816 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828
-          0.00130828 0.00130828 0.00130485 0.00130816 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828
-          0.00130828 0.00130822 0.00130828 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829
-          0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829
-          0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083
-          0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083
-          0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.00130831 0.00130831 0.00130831 0.00130831
-          0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831
-          0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831
-          0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832
-          0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832
-          0.00130832 0.00130832 0.00130832
-        </DataArray>
-        <DataArray type="Float32" Name="ViscosG" NumberOfComponents="1" format="ascii">
-          8.67141e-06 8.67143e-06 8.68382e-06 8.68384e-06 8.67149e-06 8.6839e-06 8.67158e-06 8.68399e-06 8.67169e-06 8.68411e-06 8.67182e-06 8.68424e-06
-          8.67196e-06 8.68438e-06 8.67211e-06 8.68453e-06 8.67226e-06 8.68468e-06 8.67241e-06 8.68483e-06 8.67255e-06 8.68498e-06 8.69618e-06 8.6962e-06
-          8.69627e-06 8.69638e-06 8.69651e-06 8.69665e-06 8.6968e-06 8.69695e-06 8.6971e-06 8.69725e-06 8.69741e-06 8.70871e-06 8.70852e-06 8.70861e-06
-          8.70874e-06 8.70889e-06 8.70905e-06 8.70921e-06 8.70937e-06 8.70952e-06 8.70968e-06 8.70983e-06 8.75024e-06 8.72151e-06 8.72091e-06 8.72109e-06
-          8.72128e-06 8.72145e-06 8.72162e-06 8.72179e-06 8.72195e-06 8.7221e-06 8.72226e-06 1.09277e-05 8.76365e-06 8.73358e-06 8.73345e-06 8.73367e-06
-          8.73386e-06 8.73404e-06 8.73421e-06 8.73437e-06 8.73453e-06 8.73468e-06 1.26779e-05 8.99726e-06 8.75092e-06 8.74592e-06 8.74608e-06 8.74629e-06
-          8.74647e-06 8.74664e-06 8.7468e-06 8.74696e-06 8.74711e-06 1.26519e-05 9.32498e-06 8.78276e-06 8.7589e-06 8.75854e-06 8.75874e-06 8.75892e-06
-          8.75908e-06 8.75924e-06 8.75939e-06 8.75953e-06 1.02456e-05 8.85123e-06 8.77346e-06 8.77087e-06 8.77102e-06 8.7712e-06 8.77137e-06 8.77153e-06
-          8.77168e-06 8.77182e-06 8.77196e-06 8.87719e-06 8.7869e-06 8.78328e-06 8.78336e-06 8.78353e-06 8.78369e-06 8.78384e-06 8.78399e-06 8.78412e-06
-          8.78425e-06 8.78439e-06 8.79879e-06 8.79579e-06 8.79579e-06 8.79592e-06 8.79605e-06 8.79619e-06 8.79632e-06 8.79645e-06 8.79657e-06 8.79669e-06
-          8.79681e-06 8.80836e-06 8.80832e-06 8.80839e-06 8.80848e-06 8.80858e-06 8.80869e-06 8.8088e-06 8.80891e-06 8.80902e-06 8.80913e-06 8.80924e-06
-          8.8209e-06 8.82092e-06 8.82096e-06 8.82103e-06 8.82111e-06 8.8212e-06 8.82129e-06 8.82138e-06 8.82147e-06 8.82157e-06 8.82166e-06 8.83347e-06
-          8.83348e-06 8.83351e-06 8.83356e-06 8.83362e-06 8.83369e-06 8.83377e-06 8.83385e-06 8.83393e-06 8.83401e-06 8.83409e-06 8.84602e-06 8.84602e-06
-          8.84605e-06 8.84609e-06 8.84613e-06 8.84619e-06 8.84625e-06 8.84631e-06 8.84638e-06 8.84645e-06 8.84651e-06 8.85855e-06 8.85855e-06 8.85857e-06
-          8.8586e-06 8.85863e-06 8.85868e-06 8.85873e-06 8.85878e-06 8.85883e-06 8.85888e-06 8.85894e-06 8.87106e-06 8.87107e-06 8.87108e-06 8.8711e-06
-          8.87113e-06 8.87116e-06 8.8712e-06 8.87124e-06 8.87128e-06 8.87132e-06 8.87137e-06 8.88357e-06 8.88357e-06 8.88358e-06 8.8836e-06 8.88362e-06
-          8.88364e-06 8.88367e-06 8.8837e-06 8.88373e-06 8.88376e-06 8.88379e-06 8.89607e-06 8.89608e-06 8.89608e-06 8.89609e-06 8.8961e-06 8.89612e-06
-          8.89614e-06 8.89616e-06 8.89618e-06 8.8962e-06 8.89622e-06 8.90857e-06 8.90857e-06 8.90858e-06 8.90858e-06 8.90859e-06 8.90859e-06 8.9086e-06
-          8.90861e-06 8.90862e-06 8.90863e-06 8.90864e-06 8.92107e-06 8.92107e-06 8.92107e-06 8.92107e-06 8.92107e-06 8.92107e-06 8.92107e-06 8.92107e-06
-          8.92107e-06 8.92107e-06 8.92107e-06
-        </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
@@ -245,51 +223,7 @@
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
           0.4 0.4 0.4
         </DataArray>
-        <DataArray type="Float32" Name="permeabilityXX" NumberOfComponents="1" format="ascii">
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09
-        </DataArray>
-        <DataArray type="Float32" Name="permeabilityYY" NumberOfComponents="1" format="ascii">
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09
-        </DataArray>
-        <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="phasePresence" NumberOfComponents="1" format="ascii">
           2 2 2 2 2 2 2 2 2 2 2 2
           2 2 2 2 2 2 2 2 2 2 2 2
           2 2 2 2 2 2 2 2 2 2 2 2
@@ -314,15 +248,15 @@
         <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
           283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13
           283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13
-          283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.136 283.13 283.13
-          283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.913 283.15 283.13 283.13
-          283.13 283.13 283.13 283.13 283.13 283.13 283.13 341.114 283.951 283.14 283.13 283.13
-          283.13 283.13 283.13 283.13 283.13 283.13 385.984 289.877 283.274 283.132 283.13 283.13
-          283.13 283.13 283.13 283.13 283.13 385.234 298.364 283.796 283.146 283.13 283.13 283.13
-          283.13 283.13 283.13 283.13 322.683 285.32 283.209 283.132 283.13 283.13 283.13 283.13
-          283.13 283.13 283.13 285.691 283.236 283.133 283.13 283.13 283.13 283.13 283.13 283.13
-          283.13 283.13 283.216 283.133 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13
-          283.13 283.132 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13
+          283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.142 283.13 283.13
+          283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 284.263 283.157 283.13 283.13
+          283.13 283.13 283.13 283.13 283.13 283.13 283.13 343.922 284.335 283.147 283.13 283.13
+          283.13 283.13 283.13 283.13 283.13 283.13 385.87 290.431 283.299 283.133 283.13 283.13
+          283.13 283.13 283.13 283.13 283.13 385.163 298.177 283.773 283.146 283.13 283.13 283.13
+          283.13 283.13 283.13 283.13 322.671 285.065 283.196 283.131 283.13 283.13 283.13 283.13
+          283.13 283.13 283.13 285.271 283.215 283.133 283.13 283.13 283.13 283.13 283.13 283.13
+          283.13 283.13 283.197 283.132 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13
+          283.13 283.131 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13
           283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13
           283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13
           283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13 283.13
diff --git a/test/references/steaminjectioncc-reference.vtu b/test/references/steaminjectioncc-reference.vtu
index 7c1956e3deaf2a7a2370bc2a88e89373c1d5a5d4..70ac918116b56e90a4386bafe509082f4d480acd 100644
--- a/test/references/steaminjectioncc-reference.vtu
+++ b/test/references/steaminjectioncc-reference.vtu
@@ -2,8 +2,8 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="800" NumberOfPoints="861">
-      <CellData Scalars="sw">
-        <DataArray type="Float32" Name="sw" NumberOfComponents="1" format="ascii">
+      <CellData Scalars="Sw">
+        <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
           1 1 1 1 1 1 1 1 1 1 1 1
           1 1 1 1 1 1 1 1 1 1 1 1
           1 1 1 1 1 1 1 1 1 1 1 1
@@ -72,7 +72,7 @@
           1 1 1 1 1 1 1 1 1 1 1 1
           1 1 1 1 1 1 1 1
         </DataArray>
-        <DataArray type="Float32" Name="sg" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
@@ -210,7 +210,7 @@
           102283 102283 102283 102283 102283 102283 102283 102283 102283 102283 102283 102282
           102282 102282 102282 102282 102282 102281 102281 102281
         </DataArray>
-        <DataArray type="Float32" Name="pg" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
           178915 178914 178912 178909 178906 178902 178897 178891 178885 178879 178872 178865
           178857 178849 178842 178834 178826 178818 178810 178803 176954 176953 176951 176949
           176945 176941 176936 176930 176924 176918 176911 176903 176896 176888 176880 176872
@@ -279,7 +279,76 @@
           102283 102283 102283 102283 102283 102283 102283 102283 102283 102283 102283 102282
           102282 102282 102282 102282 102282 102281 102281 102281
         </DataArray>
-        <DataArray type="Float32" Name="rhow" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 211.319 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          287.019 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 338.655 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 350.336 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
           999.728 999.728 999.728 999.728 999.728 999.728 999.728 999.728 999.728 999.728 999.728 999.728
           999.728 999.728 999.728 999.728 999.728 999.728 999.728 999.728 999.727 999.727 999.727 999.727
           999.727 999.727 999.727 999.727 999.727 999.727 999.727 999.727 999.727 999.727 999.727 999.727
@@ -348,7 +417,7 @@
           999.692 999.692 999.692 999.692 999.692 999.692 999.692 999.692 999.692 999.692 999.692 999.692
           999.692 999.692 999.692 999.692 999.692 999.692 999.692 999.692
         </DataArray>
-        <DataArray type="Float32" Name="rhog" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
           1.37174 1.37173 1.37172 1.3717 1.37167 1.37164 1.3716 1.37156 1.37151 1.37146 1.37141 1.37135
           1.3713 1.37124 1.37118 1.37112 1.37106 1.371 1.37094 1.37088 1.35671 1.3567 1.35668 1.35666
           1.35664 1.3566 1.35657 1.35652 1.35648 1.35643 1.35637 1.35632 1.35626 1.3562 1.35614 1.35608
@@ -417,7 +486,7 @@
           0.784203 0.784203 0.784203 0.784202 0.784202 0.784201 0.7842 0.784199 0.784199 0.784198 0.784196 0.784195
           0.784194 0.784193 0.784191 0.78419 0.784189 0.784187 0.784186 0.784185
         </DataArray>
-        <DataArray type="Float32" Name="MobW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
           764.38 764.38 764.38 764.38 764.38 764.38 764.38 764.38 764.38 764.38 764.38 764.38
           764.38 764.38 764.38 764.38 764.38 764.38 764.38 764.38 764.379 764.379 764.379 764.379
           764.379 764.379 764.379 764.379 764.379 764.379 764.379 764.379 764.379 764.379 764.379 764.379
@@ -486,7 +555,7 @@
           764.34 764.34 764.34 764.34 764.34 764.34 764.34 764.34 764.339 764.339 764.339 764.339
           764.339 764.339 764.339 764.339 764.339 764.339 764.339 764.339
         </DataArray>
-        <DataArray type="Float32" Name="MobG" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
@@ -555,144 +624,6 @@
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0
         </DataArray>
-        <DataArray type="Float32" Name="ViscosW" NumberOfComponents="1" format="ascii">
-          0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825
-          0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825
-          0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825
-          0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825
-          0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825 0.00130825
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130823 0.00130826 0.00130826 0.00130826
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130621 0.00130816 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826
-          0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826 0.00130826
-          0.00122413 0.00130452 0.00130815 0.00130826 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827
-          0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.000258928 0.00124998 0.00130626 0.00130821
-          0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827
-          0.00130827 0.00130827 0.00130827 0.00130827 0.000247772 0.000998563 0.00129267 0.00130782 0.00130826 0.00130827 0.00130827 0.00130827
-          0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827
-          0.000248618 0.000859732 0.00127198 0.00130655 0.00130821 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827
-          0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.0002495 0.000700015 0.0012223 0.00130166
-          0.00130794 0.00130826 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827
-          0.00130827 0.00130827 0.00130827 0.00130827 0.000250431 0.000563979 0.00111777 0.00128735 0.00130703 0.00130822 0.00130827 0.00130827
-          0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827 0.00130827
-          0.000329334 0.00100797 0.00126767 0.00130496 0.00130809 0.00130827 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828
-          0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00094886 0.00125577 0.00130321 0.00130795
-          0.00130826 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828
-          0.00130828 0.00130828 0.00130828 0.00130828 0.0012611 0.00130303 0.00130787 0.00130825 0.00130828 0.00130828 0.00130828 0.00130828
-          0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828
-          0.0013047 0.00130794 0.00130826 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828
-          0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.0013081 0.00130826 0.00130828 0.00130828
-          0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828
-          0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828 0.00130828
-          0.00130828 0.00130828 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829
-          0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829
-          0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829
-          0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829
-          0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829
-          0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829
-          0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829
-          0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829
-          0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829 0.00130829
-          0.00130829 0.00130829 0.00130829 0.00130829 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083
-          0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083
-          0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083
-          0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083
-          0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083
-          0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083
-          0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083
-          0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083
-          0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083
-          0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083 0.0013083
-          0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831
-          0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831
-          0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831
-          0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831
-          0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831
-          0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831
-          0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831
-          0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831
-          0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130831 0.00130832 0.00130832 0.00130832 0.00130832
-          0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832
-          0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832
-          0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832
-          0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832
-          0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832 0.00130832
-        </DataArray>
-        <DataArray type="Float32" Name="ViscosG" NumberOfComponents="1" format="ascii">
-          8.67529e-06 8.6753e-06 8.6753e-06 8.67531e-06 8.67532e-06 8.67534e-06 8.67535e-06 8.67537e-06 8.67539e-06 8.67541e-06 8.67543e-06 8.67545e-06
-          8.67548e-06 8.6755e-06 8.67553e-06 8.67555e-06 8.67557e-06 8.6756e-06 8.67562e-06 8.67565e-06 8.6815e-06 8.68151e-06 8.68151e-06 8.68152e-06
-          8.68153e-06 8.68154e-06 8.68156e-06 8.68158e-06 8.6816e-06 8.68162e-06 8.68164e-06 8.68166e-06 8.68169e-06 8.68171e-06 8.68174e-06 8.68176e-06
-          8.68179e-06 8.68181e-06 8.68184e-06 8.68186e-06 8.68771e-06 8.68771e-06 8.68772e-06 8.68773e-06 8.68774e-06 8.68775e-06 8.68777e-06 8.68779e-06
-          8.68781e-06 8.68783e-06 8.68785e-06 8.68787e-06 8.6879e-06 8.68792e-06 8.68795e-06 8.68797e-06 8.688e-06 8.68802e-06 8.68805e-06 8.68807e-06
-          8.69391e-06 8.69391e-06 8.69392e-06 8.69393e-06 8.69394e-06 8.69396e-06 8.69397e-06 8.69399e-06 8.69401e-06 8.69403e-06 8.69406e-06 8.69408e-06
-          8.69411e-06 8.69413e-06 8.69416e-06 8.69418e-06 8.69421e-06 8.69423e-06 8.69426e-06 8.69429e-06 8.70011e-06 8.70011e-06 8.70012e-06 8.70013e-06
-          8.70014e-06 8.70016e-06 8.70018e-06 8.7002e-06 8.70022e-06 8.70024e-06 8.70026e-06 8.70029e-06 8.70032e-06 8.70034e-06 8.70037e-06 8.70039e-06
-          8.70042e-06 8.70045e-06 8.70047e-06 8.7005e-06 8.70631e-06 8.70631e-06 8.70632e-06 8.70633e-06 8.70634e-06 8.70636e-06 8.70638e-06 8.7064e-06
-          8.70642e-06 8.70645e-06 8.70647e-06 8.7065e-06 8.70652e-06 8.70655e-06 8.70658e-06 8.7066e-06 8.70663e-06 8.70666e-06 8.70668e-06 8.70671e-06
-          8.7125e-06 8.7125e-06 8.71251e-06 8.71252e-06 8.71253e-06 8.71255e-06 8.71258e-06 8.7126e-06 8.71262e-06 8.71265e-06 8.71268e-06 8.7127e-06
-          8.71273e-06 8.71276e-06 8.71279e-06 8.71281e-06 8.71284e-06 8.71287e-06 8.7129e-06 8.71292e-06 8.71873e-06 8.71869e-06 8.7187e-06 8.71871e-06
-          8.71873e-06 8.71875e-06 8.71877e-06 8.7188e-06 8.71883e-06 8.71885e-06 8.71888e-06 8.71891e-06 8.71894e-06 8.71897e-06 8.719e-06 8.71903e-06
-          8.71905e-06 8.71908e-06 8.71911e-06 8.71914e-06 8.72684e-06 8.72497e-06 8.72488e-06 8.72489e-06 8.72491e-06 8.72494e-06 8.72497e-06 8.725e-06
-          8.72503e-06 8.72506e-06 8.72509e-06 8.72512e-06 8.72515e-06 8.72518e-06 8.72521e-06 8.72524e-06 8.72526e-06 8.72529e-06 8.72532e-06 8.72535e-06
-          8.82044e-06 8.73458e-06 8.73115e-06 8.73107e-06 8.73109e-06 8.73113e-06 8.73116e-06 8.7312e-06 8.73123e-06 8.73126e-06 8.73129e-06 8.73133e-06
-          8.73136e-06 8.73139e-06 8.73142e-06 8.73145e-06 8.73148e-06 8.7315e-06 8.73153e-06 8.73156e-06 1.24977e-05 8.79766e-06 8.73909e-06 8.73728e-06
-          8.73727e-06 8.73732e-06 8.73736e-06 8.73739e-06 8.73743e-06 8.73747e-06 8.7375e-06 8.73753e-06 8.73757e-06 8.7376e-06 8.73763e-06 8.73766e-06
-          8.73769e-06 8.73772e-06 8.73775e-06 8.73777e-06 1.26788e-05 9.12458e-06 8.758e-06 8.74382e-06 8.74346e-06 8.7435e-06 8.74355e-06 8.74359e-06
-          8.74364e-06 8.74367e-06 8.74371e-06 8.74374e-06 8.74378e-06 8.74381e-06 8.74384e-06 8.74387e-06 8.7439e-06 8.74393e-06 8.74396e-06 8.74399e-06
-          1.26666e-05 9.36982e-06 8.78542e-06 8.75119e-06 8.7497e-06 8.7497e-06 8.74975e-06 8.7498e-06 8.74984e-06 8.74988e-06 8.74992e-06 8.74995e-06
-          8.74999e-06 8.75002e-06 8.75005e-06 8.75008e-06 8.75011e-06 8.75014e-06 8.75017e-06 8.7502e-06 1.26539e-05 9.74208e-06 8.84621e-06 8.76194e-06
-          8.75614e-06 8.75591e-06 8.75596e-06 8.75601e-06 8.75605e-06 8.75609e-06 8.75613e-06 8.75617e-06 8.7562e-06 8.75623e-06 8.75626e-06 8.7563e-06
-          8.75633e-06 8.75635e-06 8.75638e-06 8.75641e-06 1.26404e-05 1.01898e-05 8.97542e-06 8.78149e-06 8.7632e-06 8.76216e-06 8.76217e-06 8.76222e-06
-          8.76226e-06 8.76231e-06 8.76234e-06 8.76238e-06 8.76241e-06 8.76245e-06 8.76248e-06 8.76251e-06 8.76254e-06 8.76257e-06 8.7626e-06 8.76263e-06
-          1.16371e-05 9.13208e-06 8.8085e-06 8.77127e-06 8.76843e-06 8.76834e-06 8.76839e-06 8.76844e-06 8.76848e-06 8.76852e-06 8.76856e-06 8.7686e-06
-          8.76863e-06 8.76866e-06 8.76869e-06 8.76872e-06 8.76875e-06 8.76878e-06 8.76881e-06 8.76884e-06 9.22985e-06 8.82757e-06 8.77905e-06 8.77473e-06
-          8.77451e-06 8.77456e-06 8.77461e-06 8.77466e-06 8.7747e-06 8.77474e-06 8.77478e-06 8.77481e-06 8.77485e-06 8.77488e-06 8.77491e-06 8.77494e-06
-          8.77497e-06 8.775e-06 8.77502e-06 8.77505e-06 8.82787e-06 8.78541e-06 8.78099e-06 8.7807e-06 8.78074e-06 8.78079e-06 8.78084e-06 8.78089e-06
-          8.78093e-06 8.78096e-06 8.781e-06 8.78103e-06 8.78107e-06 8.7811e-06 8.78113e-06 8.78116e-06 8.78118e-06 8.78121e-06 8.78124e-06 8.78127e-06
-          8.79013e-06 8.78717e-06 8.78692e-06 8.78694e-06 8.78699e-06 8.78703e-06 8.78708e-06 8.78712e-06 8.78715e-06 8.78719e-06 8.78722e-06 8.78726e-06
-          8.78729e-06 8.78732e-06 8.78734e-06 8.78737e-06 8.7874e-06 8.78743e-06 8.78745e-06 8.78748e-06 8.79329e-06 8.79315e-06 8.79317e-06 8.7932e-06
-          8.79324e-06 8.79328e-06 8.79331e-06 8.79335e-06 8.79338e-06 8.79342e-06 8.79345e-06 8.79348e-06 8.79351e-06 8.79354e-06 8.79356e-06 8.79359e-06
-          8.79362e-06 8.79364e-06 8.79367e-06 8.79369e-06 8.79941e-06 8.79941e-06 8.79943e-06 8.79946e-06 8.79949e-06 8.79952e-06 8.79955e-06 8.79958e-06
-          8.79962e-06 8.79965e-06 8.79967e-06 8.7997e-06 8.79973e-06 8.79976e-06 8.79978e-06 8.79981e-06 8.79983e-06 8.79986e-06 8.79988e-06 8.79991e-06
-          8.80567e-06 8.80568e-06 8.8057e-06 8.80572e-06 8.80574e-06 8.80577e-06 8.80579e-06 8.80582e-06 8.80585e-06 8.80587e-06 8.8059e-06 8.80593e-06
-          8.80595e-06 8.80598e-06 8.806e-06 8.80603e-06 8.80605e-06 8.80607e-06 8.8061e-06 8.80612e-06 8.81193e-06 8.81194e-06 8.81195e-06 8.81197e-06
-          8.81199e-06 8.81201e-06 8.81203e-06 8.81206e-06 8.81208e-06 8.8121e-06 8.81213e-06 8.81215e-06 8.81218e-06 8.8122e-06 8.81222e-06 8.81224e-06
-          8.81227e-06 8.81229e-06 8.81231e-06 8.81233e-06 8.81819e-06 8.81819e-06 8.8182e-06 8.81822e-06 8.81823e-06 8.81825e-06 8.81827e-06 8.81829e-06
-          8.81831e-06 8.81833e-06 8.81836e-06 8.81838e-06 8.8184e-06 8.81842e-06 8.81844e-06 8.81846e-06 8.81848e-06 8.8185e-06 8.81853e-06 8.81855e-06
-          8.82444e-06 8.82444e-06 8.82445e-06 8.82446e-06 8.82448e-06 8.82449e-06 8.82451e-06 8.82453e-06 8.82454e-06 8.82456e-06 8.82458e-06 8.8246e-06
-          8.82462e-06 8.82464e-06 8.82466e-06 8.82468e-06 8.8247e-06 8.82472e-06 8.82474e-06 8.82476e-06 8.83069e-06 8.83069e-06 8.8307e-06 8.8307e-06
-          8.83072e-06 8.83073e-06 8.83074e-06 8.83076e-06 8.83078e-06 8.83079e-06 8.83081e-06 8.83083e-06 8.83085e-06 8.83086e-06 8.83088e-06 8.8309e-06
-          8.83092e-06 8.83094e-06 8.83095e-06 8.83097e-06 8.83693e-06 8.83693e-06 8.83694e-06 8.83694e-06 8.83695e-06 8.83697e-06 8.83698e-06 8.83699e-06
-          8.83701e-06 8.83702e-06 8.83704e-06 8.83705e-06 8.83707e-06 8.83709e-06 8.8371e-06 8.83712e-06 8.83714e-06 8.83715e-06 8.83717e-06 8.83719e-06
-          8.84317e-06 8.84317e-06 8.84318e-06 8.84318e-06 8.84319e-06 8.8432e-06 8.84321e-06 8.84322e-06 8.84324e-06 8.84325e-06 8.84326e-06 8.84328e-06
-          8.84329e-06 8.84331e-06 8.84332e-06 8.84334e-06 8.84335e-06 8.84337e-06 8.84338e-06 8.8434e-06 8.84941e-06 8.84941e-06 8.84941e-06 8.84942e-06
-          8.84943e-06 8.84943e-06 8.84944e-06 8.84945e-06 8.84947e-06 8.84948e-06 8.84949e-06 8.8495e-06 8.84952e-06 8.84953e-06 8.84954e-06 8.84956e-06
-          8.84957e-06 8.84959e-06 8.8496e-06 8.84961e-06 8.85565e-06 8.85565e-06 8.85565e-06 8.85565e-06 8.85566e-06 8.85567e-06 8.85568e-06 8.85568e-06
-          8.85569e-06 8.85571e-06 8.85572e-06 8.85573e-06 8.85574e-06 8.85575e-06 8.85576e-06 8.85578e-06 8.85579e-06 8.8558e-06 8.85581e-06 8.85583e-06
-          8.86188e-06 8.86188e-06 8.86188e-06 8.86189e-06 8.86189e-06 8.8619e-06 8.86191e-06 8.86191e-06 8.86192e-06 8.86193e-06 8.86194e-06 8.86195e-06
-          8.86196e-06 8.86197e-06 8.86198e-06 8.862e-06 8.86201e-06 8.86202e-06 8.86203e-06 8.86204e-06 8.86811e-06 8.86811e-06 8.86812e-06 8.86812e-06
-          8.86812e-06 8.86813e-06 8.86814e-06 8.86814e-06 8.86815e-06 8.86816e-06 8.86817e-06 8.86818e-06 8.86819e-06 8.86819e-06 8.8682e-06 8.86821e-06
-          8.86822e-06 8.86823e-06 8.86824e-06 8.86825e-06 8.87435e-06 8.87435e-06 8.87435e-06 8.87435e-06 8.87436e-06 8.87436e-06 8.87437e-06 8.87437e-06
-          8.87438e-06 8.87438e-06 8.87439e-06 8.8744e-06 8.87441e-06 8.87442e-06 8.87442e-06 8.87443e-06 8.87444e-06 8.87445e-06 8.87446e-06 8.87447e-06
-          8.88058e-06 8.88058e-06 8.88058e-06 8.88058e-06 8.88059e-06 8.88059e-06 8.88059e-06 8.8806e-06 8.8806e-06 8.88061e-06 8.88062e-06 8.88062e-06
-          8.88063e-06 8.88064e-06 8.88064e-06 8.88065e-06 8.88066e-06 8.88067e-06 8.88067e-06 8.88068e-06 8.88681e-06 8.88681e-06 8.88681e-06 8.88681e-06
-          8.88681e-06 8.88682e-06 8.88682e-06 8.88683e-06 8.88683e-06 8.88683e-06 8.88684e-06 8.88685e-06 8.88685e-06 8.88686e-06 8.88686e-06 8.88687e-06
-          8.88688e-06 8.88688e-06 8.88689e-06 8.88689e-06 8.89304e-06 8.89304e-06 8.89304e-06 8.89304e-06 8.89304e-06 8.89305e-06 8.89305e-06 8.89305e-06
-          8.89306e-06 8.89306e-06 8.89306e-06 8.89307e-06 8.89307e-06 8.89308e-06 8.89308e-06 8.89309e-06 8.89309e-06 8.8931e-06 8.8931e-06 8.89311e-06
-          8.89927e-06 8.89927e-06 8.89927e-06 8.89927e-06 8.89927e-06 8.89927e-06 8.89928e-06 8.89928e-06 8.89928e-06 8.89928e-06 8.89929e-06 8.89929e-06
-          8.89929e-06 8.8993e-06 8.8993e-06 8.89931e-06 8.89931e-06 8.89931e-06 8.89932e-06 8.89932e-06 8.9055e-06 8.9055e-06 8.9055e-06 8.9055e-06
-          8.9055e-06 8.9055e-06 8.9055e-06 8.9055e-06 8.90551e-06 8.90551e-06 8.90551e-06 8.90551e-06 8.90552e-06 8.90552e-06 8.90552e-06 8.90552e-06
-          8.90553e-06 8.90553e-06 8.90553e-06 8.90553e-06 8.91173e-06 8.91173e-06 8.91173e-06 8.91173e-06 8.91173e-06 8.91173e-06 8.91173e-06 8.91173e-06
-          8.91173e-06 8.91173e-06 8.91173e-06 8.91174e-06 8.91174e-06 8.91174e-06 8.91174e-06 8.91174e-06 8.91174e-06 8.91174e-06 8.91175e-06 8.91175e-06
-          8.91795e-06 8.91795e-06 8.91795e-06 8.91795e-06 8.91795e-06 8.91795e-06 8.91796e-06 8.91796e-06 8.91796e-06 8.91796e-06 8.91796e-06 8.91796e-06
-          8.91796e-06 8.91796e-06 8.91796e-06 8.91796e-06 8.91796e-06 8.91796e-06 8.91796e-06 8.91796e-06
-        </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
@@ -762,144 +693,6 @@
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
         </DataArray>
-        <DataArray type="Float32" Name="permeabilityXX" NumberOfComponents="1" format="ascii">
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-        </DataArray>
-        <DataArray type="Float32" Name="permeabilityYY" NumberOfComponents="1" format="ascii">
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-          1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09 1e-09
-        </DataArray>
         <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
           2 2 2 2 2 2 2 2 2 2 2 2
           2 2 2 2 2 2 2 2 2 2 2 2