diff --git a/dumux/porousmediumflow/1pnc/CMakeLists.txt b/dumux/porousmediumflow/1pnc/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ba8341c614f1a2c797c95f5402f602025f1087b1
--- /dev/null
+++ b/dumux/porousmediumflow/1pnc/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory("implicit")
diff --git a/dumux/porousmediumflow/1pnc/implicit/CMakeLists.txt b/dumux/porousmediumflow/1pnc/implicit/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f37794a19ed79beaac326cd9ec72b7f68d7c7075
--- /dev/null
+++ b/dumux/porousmediumflow/1pnc/implicit/CMakeLists.txt
@@ -0,0 +1,10 @@
+
+#install headers
+install(FILES
+indices.hh
+model.hh
+newtoncontroller.hh
+properties.hh
+propertydefaults.hh
+volumevariables.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/1pnc/implicit)
diff --git a/dumux/porousmediumflow/1pnc/implicit/indices.hh b/dumux/porousmediumflow/1pnc/implicit/indices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b2abc84480ea1e9f4b9ef08f61e6b019af578efc
--- /dev/null
+++ b/dumux/porousmediumflow/1pnc/implicit/indices.hh
@@ -0,0 +1,60 @@
+// -*- 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 Defines the primary variable and equation indices used by
+ *        the 1pnc model
+ */
+#ifndef DUMUX_1PNC_INDICES_HH
+#define DUMUX_1PNC_INDICES_HH
+
+#include "properties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup OnePNCModel
+ * \ingroup ImplicitIndices
+ * \brief The indices for the isothermal one-phase n-component model.
+ *
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int PVOffset = 0>
+class OnePNCIndices
+{
+    //! Set the default phase used by the fluid system to the first one
+    static const int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+
+    //! Component indices
+    static const int phaseCompIdx = phaseIdx;//!< The index of the main component of the considered phase
+
+    //! Equation indices
+    static const int conti0EqIdx = PVOffset + 0; //!< Reference index for mass conservation equation.
+
+    //! Primary variable indices
+    static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
+    static const int firstMoleFracIdx = PVOffset + 1; //!< Index of the either the saturation or the mass fraction of the fluid phase
+
+    //Component indices
+    static const int firstTransportEqIdx = PVOffset + 1; //!< transport equation index
+};
+}
+
+#endif
diff --git a/dumux/porousmediumflow/1pnc/implicit/model.hh b/dumux/porousmediumflow/1pnc/implicit/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..60079c860222731f519991c4837f0e2e72d43283
--- /dev/null
+++ b/dumux/porousmediumflow/1pnc/implicit/model.hh
@@ -0,0 +1,349 @@
+// -*- 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  Base class for all models which use the single-phase, n-component fully implicit model.
+*         Adaption of the fully implicit model to the one-phase n-component flow model.
+*/
+
+#ifndef DUMUX_1PNC_MODEL_HH
+#define DUMUX_1PNC_MODEL_HH
+
+#include <dumux/porousmediumflow/implicit/velocityoutput.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/model.hh>
+
+#include "properties.hh"
+#include "indices.hh"
+// #include "localresidual.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup OnePNCModel
+ * \brief Adaption of the fully implicit scheme to the
+ *        one-phase n-component flow model.
+ *
+* This model implements a one-phase flow of a compressible fluid, that consists of n components,
+ * using a standard Darcy
+ * approach as the equation for the conservation of momentum:
+ \f[
+ v = - \frac{\textbf K}{\mu}
+ \left(\textbf{grad}\, p - \varrho {\textbf g} \right)
+ \f]
+ *
+ * Gravity can be enabled or disabled via the property system.
+ * By inserting this into the continuity equation, one gets
+ \f[
+ \phi\frac{\partial \varrho}{\partial t} - \text{div} \left\{
+   \varrho \frac{\textbf K}{\mu}  \left(\textbf{grad}\, p - \varrho {\textbf g} \right)
+ \right\} = q \;,
+ \f]
+ *
+ * The transport of the components \f$\kappa \in \{ w, a, ... \}\f$ is described by the following equation:
+ \f[
+ \phi \frac{ \partial \varrho X^\kappa}{\partial t}
+ - \text{div} \left\lbrace \varrho X^\kappa \frac{{\textbf K}}{\mu} \left( \textbf{grad}\, p -
+ \varrho {\textbf g} \right)
+ + \varrho D^\kappa_\text{pm} \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa \right\rbrace = q.
+ \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.
+ * The model is able to use either mole or mass fractions. The property useMoles can be set to either true or false in the
+ * problem file. Make sure that the according units are used in the problem setup. useMoles is set to true by default.
+ *
+ * The primary variables are the pressure \f$p\f$ and the mole or mass fraction of dissolved component \f$x\f$.
+ */
+
+template<class TypeTag>
+class OnePNCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
+{
+    using ParentType = typename GET_PROP_TYPE(TypeTag, BaseModel);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using NonIsothermalModel = Dumux::NonIsothermalModel<TypeTag>;
+
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+//    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+
+    static const int phaseIdx = Indices::phaseIdx;
+
+    enum { dim = GridView::dimension };
+    enum { dimWorld = GridView::dimensionworld };
+//
+//     enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+//
+    enum {  numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+//
+//      enum {
+//             pressureIdx = Indices::pressureIdx,
+//             firstMoleFracIdx = Indices::firstMoleFracIdx,
+//     };
+
+    typedef typename GridView::template Codim<dim>::Entity Vertex;
+    using Element = typename GridView::template Codim<0>::Entity;
+
+//     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+    using CoordScalar = typename GridView::ctype;
+    using Tensor = Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld>;
+
+    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+    enum { dofCodim = isBox ? dim : 0 };
+
+public:
+    /*!
+     * \brief Apply the initial conditions to the model.
+     *
+     * \param problem The object representing the problem which needs to
+     *             be simulated.
+     */
+    void init(Problem &problem)
+    {
+        ParentType::init(problem);
+
+        // register standardized vtk output fields
+
+        auto& vtkOutputModule = problem.vtkOutputModule();
+        vtkOutputModule.addSecondaryVariable("p", [](const VolumeVariables& v){ return v.pressure(phaseIdx); });
+        vtkOutputModule.addSecondaryVariable("rho", [](const VolumeVariables& v){ return v.density(phaseIdx); });
+        vtkOutputModule.addSecondaryVariable("porosity", [](const VolumeVariables& v){ return v.porosity(); });
+        vtkOutputModule.addSecondaryVariable("temperature", [](const VolumeVariables& v){ return v.temperature(); });
+
+        vtkOutputModule.addSecondaryVariable("Kxx",
+                                             [this](const VolumeVariables& v){ return this->perm_(v.permeability())[0][0]; });
+        if (dim >= 2)
+            vtkOutputModule.addSecondaryVariable("Kyy",
+                                                 [this](const VolumeVariables& v){ return this->perm_(v.permeability())[1][1]; });
+        if (dim >= 3)
+            vtkOutputModule.addSecondaryVariable("Kzz",
+                                                 [this](const VolumeVariables& v){ return this->perm_(v.permeability())[2][2]; });
+
+       for (int i = 0; i < numComponents; ++i)
+           vtkOutputModule.addSecondaryVariable("x" + FluidSystem::componentName(i),
+                                                [i](const VolumeVariables& v){ return v.moleFraction(phaseIdx, i); });
+
+       for (int i = 0; i < numComponents; ++i)
+           vtkOutputModule.addSecondaryVariable("m^w_" + FluidSystem::componentName(i),
+                                                 [i](const VolumeVariables& v){ return v.molarity(phaseIdx,i); });
+
+        NonIsothermalModel::maybeAddTemperature(vtkOutputModule);
+    }
+
+    /*!
+     * \brief Compute the total storage of all conservation quantities
+     *
+     * \param storage Contains the storage of each component
+     * \param phaseIdx The phase index
+     */
+//     void globalPhaseStorage(PrimaryVariables &storage, int phaseIdx)
+//     {
+//         storage = 0;
+//         for (const auto& element : elements(this->gridView_()))
+//         {
+//             if(element.partitionType() == Dune::InteriorEntity)
+//             {
+//                 this->localResidual().evalPhaseStorage(element, phaseIdx);
+//
+//                 for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
+//                     storage += this->localResidual().storageTerm()[i];
+//             }
+//         }
+//         this->gridView_().comm().sum(storage);
+//     }
+
+    /*!
+     * \brief Called by the update() method if applying the newton
+     *         method was unsuccessful.
+     */
+    void updateFailed()
+    {
+        ParentType::updateFailed();
+    }
+
+    /*!
+     * \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();
+//     }
+
+    /*!
+     * \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<Scalar, 1> > ScalarField;
+//         typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > VectorField;
+//
+//         // get the number of degrees of freedom
+//         auto numDofs = this->numDofs();
+//
+//         ScalarField *pressure            = writer.allocateManagedBuffer (numDofs);
+//         ScalarField *rho          = writer.allocateManagedBuffer (numDofs);
+//         ScalarField *temperature   = writer.allocateManagedBuffer (numDofs);
+//         ScalarField *poro          = writer.allocateManagedBuffer (numDofs);
+//         VectorField *velocity = 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)
+//             {
+//                 (*velocity)[i] = Scalar(0);
+//             }
+//         }
+//
+//         ScalarField *moleFraction[numComponents];
+// //         for (int i = 0; i < numPhases; ++i)
+//             for (int j = 0; j < numComponents; ++j)
+//                 moleFraction[j] = writer.allocateManagedBuffer(numDofs);
+//
+//         ScalarField *molarity[numComponents];
+//         for (int j = 0; j < numComponents ; ++j)
+//             molarity[j] = writer.allocateManagedBuffer(numDofs);
+//
+//         ScalarField *Perm[dim];
+//         for (int j = 0; j < dim; ++j) //Permeability only in main directions xx and yy
+//             Perm[j] = writer.allocateManagedBuffer(numDofs);
+//
+//         auto numElements = this->gridView_().size(0);
+//         ScalarField *rank = writer.allocateManagedBuffer (numElements);
+//
+//         for (const auto& element : elements(this->gridView_()))
+//         {
+//             auto eIdxGlobal = this->problem_().elementMapper().index(element);
+//             (*rank)[eIdxGlobal] = 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)
+//             {
+//                 auto dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
+//
+//                 GlobalPosition globalPos = fvGeometry.subContVol[scvIdx].global;
+//                 (*pressure)[dofIdxGlobal]             = elemVolVars[scvIdx].pressure(phaseIdx);
+//                 (*rho)[dofIdxGlobal]           = elemVolVars[scvIdx].density(phaseIdx);
+//                 (*poro)[dofIdxGlobal]           = elemVolVars[scvIdx].porosity();
+//                 (*temperature)[dofIdxGlobal]    = elemVolVars[scvIdx].temperature();
+//
+//                 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+//                     (*moleFraction[compIdx])[dofIdxGlobal]= elemVolVars[scvIdx].moleFraction(compIdx);
+//
+//                 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+//                     (*molarity[compIdx])[dofIdxGlobal] = (elemVolVars[scvIdx].molarity(compIdx));
+//
+//                 Tensor K = perm_(this->problem_().spatialParams().intrinsicPermeability(element, fvGeometry, scvIdx));
+//
+//                 for (int j = 0; j<dim; ++j){
+//                     (*Perm[j])[dofIdxGlobal] = K[j][j];
+//                 }
+//             };
+//
+//             // velocity output
+//             if(velocityOutput.enableOutput()){
+//                 velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, element, phaseIdx);
+//             }
+//
+//         } // loop over element
+//
+//         writer.attachDofData(*pressure, "p", isBox);
+//         writer.attachDofData(*rho, "rho", isBox);
+//         writer.attachDofData(*poro, "porosity", isBox);
+//         writer.attachDofData(*temperature, "temperature", isBox);
+//         writer.attachDofData(*Perm[0], "Kxx", isBox);
+//         if (dim >= 2)
+//             writer.attachDofData(*Perm[1], "Kyy", isBox);
+//         if (dim == 3)
+//             writer.attachDofData(*Perm[2], "Kzz", isBox);
+//
+//             for (int j = 0; j < numComponents; ++j)
+//             {
+//                 std::ostringstream oss;
+//                 oss << "x"
+//                     << FluidSystem::componentName(j);
+// //                     << FluidSystem::phaseName(i);
+//                 writer.attachDofData(*moleFraction[j], oss.str(), isBox);
+//             }
+//
+//         for (int j = 0; j < numComponents; ++j)
+//         {
+//             std::ostringstream oss;
+//             oss << "m^w_"
+//                 << FluidSystem::componentName(j);
+//             writer.attachDofData(*molarity[j], oss.str().c_str(), isBox);
+//         }
+//
+//         if (velocityOutput.enableOutput()) // check if velocity output is demanded
+//         {
+//             writer.attachDofData(*velocity,  "velocity", isBox, dim);
+//         }
+//
+//         writer.attachCellData(*rank, "process rank");
+//     }
+
+
+
+private :
+    Tensor perm_(Scalar perm)
+    {
+        Tensor K(0.0);
+
+        for(int i=0; i<dim; i++)
+            K[i][i] = perm;
+
+       return K;
+    }
+
+    Tensor perm_(Tensor perm)
+    {
+       return perm;
+    }
+
+};
+
+}
+
+#include "propertydefaults.hh"
+
+#endif
diff --git a/dumux/porousmediumflow/1pnc/implicit/newtoncontroller.hh b/dumux/porousmediumflow/1pnc/implicit/newtoncontroller.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3308119b0cdc1f43d41b66761267f1701a528f59
--- /dev/null
+++ b/dumux/porousmediumflow/1pnc/implicit/newtoncontroller.hh
@@ -0,0 +1,145 @@
+// -*- 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 one-phase n-component specific controller for the newton solver.
+ */
+#ifndef DUMUX_1PNC_NEWTON_CONTROLLER_HH
+#define DUMUX_1PNC_NEWTON_CONTROLLER_HH
+
+#include "properties.hh"
+
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+namespace Dumux {
+/*!
+ * \ingroup Newton
+ * \ingroup OnePNCModel
+ * \brief A one-phase n-component 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 OnePNCNewtonController : public NewtonController<TypeTag>
+{
+    typedef NewtonController<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+//     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+//     typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+//     typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+//     typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+//     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+//     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+//     typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+//     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+
+public:
+  OnePNCNewtonController(Problem &problem)
+      : ParentType(problem)
+    {};
+
+    /*!
+     * \brief
+     * Suggest a new time step size based either on the number of newton
+     * iterations required or on the variable switch
+     *
+     * \param uCurrentIter The current global solution vector
+     * \param uLastIter The previous global solution vector
+     *
+     */
+//     void newtonEndStep(SolutionVector &uCurrentIter,
+//                        const SolutionVector &uLastIter)
+//     {
+//         int succeeded;
+//         try {
+//             // call the method of the base class
+//             this->method().model().updateStaticData(uCurrentIter, uLastIter);
+//             ParentType::newtonEndStep(uCurrentIter, uLastIter);
+
+
+//TODO
+            // loop over all element analogous to model.hh -> elemVolVars
+            // call solDependentSource
+            // evaluate if source term is admissible using some of the calculation from solDependentSource
+            // if source is too large, throw NumericalProblem like in region2.hh
+//             for (const auto& element : elements(this->problem_().gridView()))
+//             {
+//                 FVElementGeometry fvGeometry;
+//                 fvGeometry.update(this->problem_().gridView(), element);
+//
+//                 ElementVolumeVariables elemVolVars;
+//                 elemVolVars.update(this->problem_(),
+//                                   element,
+//                                   fvGeometry,
+//                                   false /* oldSol? */);
+//
+//                 for (unsigned int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+//                 {
+//                   PrimaryVariables source(0.0);
+//
+//                   this->problem_().solDependentSource(source, element, fvGeometry, scvIdx, elemVolVars);
+//
+//                   const auto& volVars = elemVolVars[scvIdx];
+//                   Scalar moleFracCaO_sPhase = volVars.precipitateVolumeFraction(cPhaseIdx)*volVars.molarDensity(cPhaseIdx)
+//                                             /(volVars.precipitateVolumeFraction(hPhaseIdx)*volVars.molarDensity(hPhaseIdx)
+//                                             + volVars.precipitateVolumeFraction(cPhaseIdx)*volVars.molarDensity(cPhaseIdx));
+//                   // if (isCharge = true)
+//                   if (- source[CaOIdx]*this->problem_().timeManager().timeStepSize() + moleFracCaO_sPhase* volVars.molarDensity(cPhaseIdx)
+//                       < 0 + 1e-6){
+//                       DUNE_THROW(NumericalProblem,
+//                           "Source term delivers unphysical value");
+//                   }
+//                  }
+//              }
+
+
+
+//             succeeded = 1;
+//             if (this->gridView_().comm().size() > 1)
+//                 succeeded = this->gridView_().comm().min(succeeded);
+//         }
+//         catch (Dumux::NumericalProblem &e)
+//         {
+//             std::cout << "rank " << this->problem_().gridView().comm().rank()
+//                       << " caught an exception while updating:" << e.what()
+//                       << "\n";
+//             succeeded = 0;
+//             if (this->gridView_().comm().size() > 1)
+//                 succeeded = this->gridView_().comm().min(succeeded);
+//         }
+//
+//         if (!succeeded) {
+//             DUNE_THROW(NumericalProblem,
+//                        "A process did not succeed in linearizing the system");
+//         }
+//     }
+
+    bool newtonConverged()
+    {
+
+        return ParentType::newtonConverged();
+    }
+};
+}
+
+#endif
diff --git a/dumux/porousmediumflow/1pnc/implicit/properties.hh b/dumux/porousmediumflow/1pnc/implicit/properties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..38a3bbaed8c7083bfe5fdfbb87277ce7bfe75763
--- /dev/null
+++ b/dumux/porousmediumflow/1pnc/implicit/properties.hh
@@ -0,0 +1,82 @@
+// -*- 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 Properties
+ * \ingroup ImplicitProperties
+ * \ingroup OnePNCModel
+ *
+ * \file
+ *
+ * \brief Defines the properties required for the one-phase n-component
+ *        fully implicit model.
+ */
+#ifndef DUMUX_1PNC_PROPERTIES_HH
+#define DUMUX_1PNC_PROPERTIES_HH
+
+#include <dumux/implicit/box/properties.hh>
+#include <dumux/implicit/cellcentered/properties.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/properties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the implicit isothermal one phase n component problems
+NEW_TYPE_TAG(OnePNC);
+NEW_TYPE_TAG(BoxOnePNC, INHERITS_FROM(BoxModel, OnePNC));
+NEW_TYPE_TAG(CCOnePNC, INHERITS_FROM(CCModel, OnePNC));
+
+//! The type tag for the implicit non-isothermal one phase n component problems
+NEW_TYPE_TAG(OnePNCNI, INHERITS_FROM(OnePNC, NonIsothermal));
+NEW_TYPE_TAG(BoxOnePNCNI, INHERITS_FROM(BoxModel, OnePNCNI));
+NEW_TYPE_TAG(CCOnePNCNI, INHERITS_FROM(CCModel, OnePNCNI));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(NumPhases);   //!< Number of fluid phases in the system
+NEW_PROP_TAG(PhaseIdx); //!< A phase index in to allow that a two-phase fluidsystem is used
+NEW_PROP_TAG(NumComponents); //!< Number of fluid components in the system
+NEW_PROP_TAG(Indices); //!< Enumerations for the model
+//NEW_PROP_TAG(OnePNCIndices); //!< Enumerations for the 1pncMin models
+NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters
+NEW_PROP_TAG(EffectiveDiffusivityModel); //!< The employed model for the computation of the effective diffusivity
+NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
+NEW_PROP_TAG(FluidState); //!< Type of the fluid state to be used
+//NEW_PROP_TAG(Chemistry); //!< The chemistry class with which solves equlibrium reactions
+//NEW_PROP_TAG(ReplaceCompEqIdx); //!< The index of the total mass balance equation, if one component balance is replaced (ReplaceCompEqIdx < NumComponents)
+//NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output
+NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
+NEW_PROP_TAG(UseMoles); //!< Defines whether mole (true) or mass (false) fractions are used
+// NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the upwind weight for the mass conservation equations
+// NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< The value of the upwind parameter for the mobility
+//NEW_PROP_TAG(BaseFluxVariables); //! The base flux variables
+NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
+NEW_PROP_TAG(TauTortuosity); //!< Tortuosity value (tau) used in macroscopic diffusion
+
+}
+}
+
+#endif
diff --git a/dumux/porousmediumflow/1pnc/implicit/propertydefaults.hh b/dumux/porousmediumflow/1pnc/implicit/propertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9981a668542056ac5edbf6af1a9c5e4b8e2dc96c
--- /dev/null
+++ b/dumux/porousmediumflow/1pnc/implicit/propertydefaults.hh
@@ -0,0 +1,213 @@
+// -*- 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 Properties
+ * \ingroup ImplicitProperties
+ * \ingroup OnePNCModel
+ * \file
+ *
+ * \brief Defines some default values for the properties required by the
+ *        one-phase n-component fully implicit model.
+ */
+#ifndef DUMUX_1PNC_PROPERTY_DEFAULTS_HH
+#define DUMUX_1PNC_PROPERTY_DEFAULTS_HH
+
+#include "indices.hh"
+#include "model.hh"
+//#include "fluxvariables.hh"
+#include "volumevariables.hh"
+#include "properties.hh"
+#include "newtoncontroller.hh"
+
+#include <dumux/porousmediumflow/compositional/localresidual.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
+#include <dumux/material/spatialparams/implicit1p.hh>
+#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
+//#include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
+#include <dumux/material/fluidmatrixinteractions/1p/thermalconductivityaverage.hh>
+#include <dumux/material/fluidstates/compositional.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+
+/*!
+ * \brief Set the property for the number of equations: For each existing component one equation has to be solved.
+ */
+SET_PROP(OnePNC, NumEq)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+public:
+    static const int value = FluidSystem::numComponents;
+};
+
+
+SET_INT_PROP(OnePNC, NumPhases,1); //!< The number of phases in the 1pnc model is 1
+
+SET_BOOL_PROP(OnePNC, UseMoles, true); //!< Define that mole fractions are used in the balance equations
+
+/*!
+ * \brief Set the property for the number of components.
+ *
+ * We just forward the number from the fluid system
+ *
+ */
+SET_PROP(OnePNC, NumComponents)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+    static const int value = FluidSystem::numComponents;
+;
+
+};
+//! Set as default that no component mass balance is replaced by the total mass balance
+SET_PROP(OnePNC, ReplaceCompEqIdx)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+    static const int value = FluidSystem::numComponents;
+};
+
+
+/*!
+ * \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(OnePNC, FluidState){
+    private:
+        typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+        typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    public:
+        typedef Dumux::CompositionalFluidState<Scalar, FluidSystem> type;
+};
+
+//! Use the 1pnc local residual
+// SET_TYPE_PROP(OnePNC, LocalResidual, OnePNCLocalResidual<TypeTag>);
+
+//! Use the 1pnc newton controller
+SET_TYPE_PROP(OnePNC, NewtonController, OnePNCNewtonController<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(OnePNC, Model, OnePNCModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(OnePNC, VolumeVariables, OnePNCVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+// SET_TYPE_PROP(OnePNC, FluxVariables, OnePNCFluxVariables<TypeTag>);
+
+//! define the base flux variables to realize Darcy flow
+// SET_TYPE_PROP(OnePNC, BaseFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
+
+//! the upwind weight for the mass conservation equations.
+// SET_SCALAR_PROP(OnePNC, ImplicitMassUpwindWeight, 1.0);
+
+// //! Set default mobility upwind weight to 1.0, i.e. fully upwind
+// SET_SCALAR_PROP(OnePNC, ImplicitMobilityUpwindWeight, 1.0);
+
+//! The indices required by the isothermal 2pnc model
+SET_TYPE_PROP(OnePNC, Indices, OnePNCIndices <TypeTag, /*PVOffset=*/0>);
+
+//! Set the phaseIndex per default to zero (important for one-phase fluidsystems).
+SET_INT_PROP(OnePNC, PhaseIdx, 0);
+
+//! Use the ImplicitSpatialParams by default
+SET_TYPE_PROP(OnePNC, SpatialParams, ImplicitSpatialParamsOneP<TypeTag>);
+
+//! Enable gravity by default
+SET_BOOL_PROP(OnePNC, ProblemEnableGravity, false);
+
+//! Disable velocity output by default
+// SET_BOOL_PROP(OnePNC, VtkAddVelocity, false);
+
+//! Somerton is used as default model to compute the effective thermal heat conductivity
+// SET_PROP(OnePNCNI, ThermalConductivityModel)
+// {
+// private:
+//     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+//     typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+// public:
+//     typedef ThermalConductivityAverage<Scalar> type;
+// };
+
+//! average is used as default model to compute the effective thermal heat conductivity
+SET_PROP(OnePNC, ThermalConductivityModel)
+{ private :
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+  public:
+    typedef ThermalConductivityAverage<Scalar> type;
+};
+
+
+//////////////////////////////////////////////////////////////////
+// Property values for isothermal model required for the general non-isothermal model
+//////////////////////////////////////////////////////////////////
+
+// set isothermal Model
+SET_TYPE_PROP(OnePNCNI, IsothermalModel, OnePNCModel<TypeTag>);
+
+//set isothermal VolumeVariables
+SET_TYPE_PROP(OnePNC, IsothermalVolumeVariables, OnePNCVolumeVariables<TypeTag>);
+
+//! temperature is already written by the isothermal model
+// SET_BOOL_PROP(OnePNCNI, NiOutputLevel, 0);
+
+// set isothermal FluxVariables
+// SET_TYPE_PROP(OnePNCNI, IsothermalFluxVariables, OnePNCFluxVariables<TypeTag>);
+
+//set isothermal VolumeVariables
+SET_TYPE_PROP(OnePNCNI, IsothermalVolumeVariables, OnePNCVolumeVariables<TypeTag>);
+
+//set isothermal LocalResidual
+// SET_TYPE_PROP(OnePNCNI, IsothermalLocalResidual, OnePNCLocalResidual<TypeTag>);
+
+//set isothermal Indices
+SET_TYPE_PROP(OnePNCNI, IsothermalIndices, OnePNCIndices<TypeTag, /*PVOffset=*/0>);
+
+//set isothermal NumEq
+SET_PROP(OnePNCNI, IsothermalNumEq)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+    static const int value = FluidSystem::numComponents;
+};
+
+//physical processes to be considered by the non-isothermal model
+SET_BOOL_PROP(OnePNCNI, EnableAdvection, true);
+SET_BOOL_PROP(OnePNCNI, EnableMolecularDiffusion, true);
+SET_BOOL_PROP(OnePNCNI, EnableEnergyBalance, true);
+
+}
+}
+
+#endif
diff --git a/dumux/porousmediumflow/1pnc/implicit/volumevariables.hh b/dumux/porousmediumflow/1pnc/implicit/volumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5614cab0a5435b6fe8af9db6f474537404c764dc
--- /dev/null
+++ b/dumux/porousmediumflow/1pnc/implicit/volumevariables.hh
@@ -0,0 +1,403 @@
+// -*- 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 Quantities required by the single-phase, n-component box
+ *        model defined on a vertex.
+ */
+#ifndef DUMUX_1PNC_VOLUME_VARIABLES_HH
+#define DUMUX_1PNC_VOLUME_VARIABLES_HH
+
+#include <dumux/discretization/volumevariables.hh>
+
+#include "properties.hh"
+
+
+
+namespace Dumux
+{
+
+namespace Properties
+{
+NEW_PROP_TAG(IsothermalVolumeVariables);
+}
+
+/*!
+ * \ingroup OnePNCModel
+ * \ingroup ImplicitVolumeVariables
+ * \brief Contains the quantities which are are constant within a
+ *        finite volume in the one-phase, n-component model.
+ *
+ * \note The return functions for the fluid state variables always forward to the actual
+ *       fluid state using the phaseIdx from the DuMuX property system. Furthermore, the
+ *       default value is not used, but is only here to enable calling these functions
+ *       without handing in a phase index (as in a single-phasic context there is only one phase).
+ *       This way one can use two-phase fluid systems for this one-phasic flow and transport
+ *       model by specifying which phase is present through the DuMuX property system.
+ */
+template <class TypeTag>
+class OnePNCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
+{
+    using ParentType = ImplicitVolumeVariables<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 PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+
+//TODO
+// ? deprecated in next ??
+    using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+//     typedef typename GET_PROP_TYPE(TypeTag, IsothermalVolumeVariables) IsothermalVolumeVariables;
+
+// new and necessary in next ??
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+    using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
+    using PermeabilityType = typename SpatialParams::PermeabilityType;
+
+    enum
+    {
+        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+        numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
+
+        phaseIdx = Indices::phaseIdx,
+        phaseCompIdx = Indices::phaseCompIdx,
+        firstTransportEqIdx = Indices::firstTransportEqIdx,
+
+        // primary variable indices
+        pressureIdx = Indices::pressureIdx,
+        firstMoleFracIdx = Indices::firstMoleFracIdx,
+
+    };
+
+    static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+    static const int dim = GridView::dimension;
+    static const int dimWorld = GridView::dimensionworld;
+
+    using Element = typename GridView::template Codim<0>::Entity;
+
+// new
+    using DimVector = Dune::FieldVector<Scalar,dim>;
+    using GlobalPosition = Dune::FieldVector<Scalar,dimWorld>;
+//     typedef typename Grid::ctype CoordScalar; ?? deprecated ?
+
+//     enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+//     enum { dofCodim = isBox ? dim : 0 };
+public:
+
+    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+    /*!
+     * \copydoc ImplicitVolumeVariables::update
+     * \param priVars The primary Variables
+     */
+    void update(const ElementSolutionVector &elemSol,
+                const Problem &problem,
+                const Element &element,
+                const SubControlVolume &scv
+                /*bool isOldSol ??*/)
+    {
+        ParentType::update(elemSol, problem, element, scv);
+
+        completeFluidState(elemSol, problem, element, scv, fluidState_/*, isOldSol*/);
+
+        porosity_ = problem.spatialParams().porosity(element, scv, elemSol);
+        // dispersivity_ = problem.spatialParams().dispersivity(element, scv, elemSol);
+        permeability_ = problem.spatialParams().permeability(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.updatePhase(fluidState_, phaseIdx);
+
+        int compIIdx = phaseCompIdx;
+
+        for (unsigned int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
+        {
+            diffCoeff_[compJIdx] = 0.0;
+            if(compIIdx!= compJIdx)
+                {
+                diffCoeff_[compJIdx] = FluidSystem::binaryDiffusionCoefficient(fluidState_,
+                                                             paramCache,
+                                                             phaseIdx,
+                                                             compIIdx,
+                                                             compJIdx);
+                }
+//             Valgrind::CheckDefined(diffCoeff_[compJIdx]);
+        }
+
+        // energy related quantities not contained in the fluid state
+//         asImp_().callProtectedUpdateEnergy(elemSol, problem, element, fvGeometry, scvIdx, isOldSol);
+    }
+
+   /*!
+    * \copydoc ImplicitModel::completeFluidState
+    * \param isOldSol Specifies whether this is the previous solution or the current one
+    * \param priVars The primary Variables
+    */
+    static void completeFluidState(const ElementSolutionVector &elemSol,
+                                   const Problem& problem,
+                                   const Element& element,
+                                   const SubControlVolume &scv,
+                                   FluidState& fluidState)
+
+    {
+//    old
+//         Scalar t = IsothermalVolumeVariables::callProtectedTemperature(priVars, problem, element,
+//                                                                        fvGeometry, scvIdx);
+        Scalar t = ParentType::temperature(elemSol, problem, element, scv);
+        fluidState.setTemperature(t);
+        fluidState.setSaturation(phaseIdx, 1.);
+
+        const auto& priVars = ParentType::extractDofPriVars(elemSol, scv); // new in next ??
+        fluidState.setPressure(phaseIdx, priVars[pressureIdx]);
+
+        // calculate the phase composition
+
+        Dune::FieldVector<Scalar, numComponents> moleFrac;
+
+        Scalar sumMoleFracNotWater = 0;
+
+        for (int compIdx=firstMoleFracIdx; compIdx<numComponents; ++compIdx){
+             moleFrac[compIdx] = priVars[compIdx];
+             sumMoleFracNotWater +=moleFrac[compIdx];
+            }
+        moleFrac[0] = 1- sumMoleFracNotWater;
+
+        // Set fluid state mole fractions
+        for (int compIdx=0; compIdx<numComponents; ++compIdx)
+        {
+            fluidState.setMoleFraction(phaseIdx, compIdx, moleFrac[compIdx]);
+        }
+
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updateAll(fluidState);
+
+        Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
+        Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
+
+        fluidState.setDensity(phaseIdx, rho);
+        fluidState.setViscosity(phaseIdx, mu);
+
+        // compute and set the enthalpy
+        Scalar h = Implementation::enthalpy(fluidState, paramCache, phaseIdx);
+        fluidState.setEnthalpy(phaseIdx, h);
+    }
+
+    /*!
+     * \brief Return the fluid configuration at the given primary
+     *        variables
+     */
+    const FluidState &fluidState() const
+    { return fluidState_; }
+
+    /*!
+     * \brief Return density \f$\mathrm{[kg/m^3]}\f$ the of the fluid phase.
+     *
+     * We always forward to the fluid state with the phaseIdx property (see class description).
+     */
+    Scalar density(int phaseIdx) const
+    { return fluidState_.density(phaseIdx); }
+
+    /*!
+     * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ the of the fluid phase.
+     *
+     * We always forward to the fluid state with the phaseIdx property (see class description).
+     */
+    Scalar molarDensity(int phaseIdx = 0) const
+    { return fluidState_.molarDensity(phaseIdx); }
+
+    /*!
+     * \brief Return the saturation
+     *
+     * This method is here for compatibility reasons with other models. The saturation
+     * is always 1.0 in a one-phasic context.
+     */
+    Scalar saturation(int pIdx = 0) const
+    { return 1.0; }
+
+     /*!
+      * \brief Return mole fraction \f$\mathrm{[mol/mol]}\f$ of a component in the phase.
+      *
+      * \param phaseIdx the index of the fluid phase
+      * \param compIdx the index of the component
+      *
+      * We always forward to the fluid state with the phaseIdx property (see class description).
+      */
+     Scalar moleFraction(int phaseIdx, int compIdx) const
+     { return fluidState_.moleFraction(phaseIdx, compIdx); }
+
+    /*!
+     * \brief Return the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
+     *        the control volume.
+     *
+     * \param phaseIdx The phase index
+     *
+     * We always forward to the fluid state with the phaseIdx property (see class description).
+     */
+    Scalar pressure(int pIdx = 0) const
+    { return fluidState_.pressure(phaseIdx); }
+
+    /*!
+     * \brief Return temperature \f$\mathrm{[K]}\f$ inside the sub-control volume.
+     *
+     * Note that we assume thermodynamic equilibrium, i.e. the
+     * temperature of the rock matrix and of all fluid phases are
+     * identical.
+     */
+    Scalar temperature() const
+    { return fluidState_.temperature(); }
+
+    /*!
+     * \brief Returns the mobility \f$\mathrm{[1/(Pa s)]}\f$.
+     *
+     * The term mobility is usually not employed in the one phase context.
+     * The method is here for compatibility reasons with other models.
+     *
+     * We always forward to the fluid state with the phaseIdx property (see class description).
+     */
+    Scalar mobility(int pIdx = 0) const
+    { return 1.0/fluidState_.viscosity(phaseIdx); }
+
+    /*!
+     * \brief Return the average porosity \f$\mathrm{[-]}\f$ within the control volume.
+     */
+    Scalar porosity() const
+    { return porosity_; }
+
+
+//     /*!
+//      * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
+//      */
+//     Scalar diffCoeff(int compIdx) const
+//     { return diffCoeff_[compIdx]; }
+
+    /*!
+     * \brief Return the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ in the fluid.
+     */
+    Scalar diffusionCoefficient(int pIdx, int compIdx) const
+    { return diffCoeff_; }
+
+    /*!
+     * \brief Returns the molarity of a component in the phase
+     *
+     * \param phaseIdx the index of the fluid phase
+     * \param compIdx the index of the component
+     */
+     Scalar molarity(int compIdx) const // [moles/m^3]
+    { return fluidState_.molarity(phaseIdx, compIdx);}
+// old    { return this->fluidState_.molarity(phaseIdx, compIdx);}
+     /*!
+      * \brief Returns the mass fraction of a component in the phase
+      *
+      * \param phaseIdx the index of the fluid phase
+      * \param compIdx the index of the component
+      */
+//      Scalar massFraction(int compIdx) const
+//      {
+//         return this->fluidState_.massFraction(phaseIdx, compIdx);
+//      }
+
+
+   /*!
+    * Circumvents the inheritance architecture of the ninisothermal model
+    */
+//     static Scalar callProtectedTemperature(const PrimaryVariables &priVars,
+//                                            const Problem& problem,
+//                                            const Element &element,
+//                                            const FVElementGeometry &fvGeometry,
+//                                            int scvIdx)
+//     {
+//          return Implementation::temperature_(priVars, problem,element, fvGeometry, scvIdx);
+//     }
+//
+//    /*!
+//     * Circumvents the inheritance architecture of the ninisothermal model
+//     */
+//    void callProtectedUpdateEnergy(const PrimaryVariables &priVars,
+//                                    const Problem &problem,
+//                                    const Element &element,
+//                                    const FVElementGeometry &fvGeometry,
+//                                    const int scvIdx,
+//                                    bool isOldSol)
+//     {
+//         asImp_().updateEnergy_(priVars, problem,element, fvGeometry, scvIdx, isOldSol);
+//     };
+
+    /*!
+     * \brief Returns the permeability within the control volume in \f$[m^2]\f$.
+     */
+    const PermeabilityType& permeability() const
+    { return permeability_; }
+
+protected:
+
+    static Scalar temperature_(const PrimaryVariables &priVars,
+                               const Problem& problem,
+                               const Element &element,
+                               const FVElementGeometry &fvGeometry,
+                               int scvIdx)
+    {
+         return problem.temperatureAtPos(fvGeometry.subContVol[scvIdx].global);
+    }
+
+//     template<class ParameterCache>
+//     static Scalar enthalpy_(const FluidState& fluidState,
+//                             const ParameterCache& paramCache,
+//                             int phaseIdx)
+//     {
+//         return 0;
+//     }
+
+    /*!
+        * \brief Called by update() to compute the energy related quantities
+        */
+//     void updateEnergy_(const PrimaryVariables &priVars,
+//                         const Problem &problem,
+//                         const Element &element,
+//                         const FVElementGeometry &fvGeometry,
+//                         const int scvIdx,
+//                         bool isOldSol)
+//     { };
+//
+    Scalar porosity_;        //!< Effective porosity within the control volume
+    PermeabilityType permeability_;
+    Scalar density_;
+    FluidState fluidState_;
+    Dune::FieldVector<Scalar, numComponents> diffCoeff_;
+
+private:
+    Implementation &asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    const Implementation &asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/porousmediumflow/CMakeLists.txt b/dumux/porousmediumflow/CMakeLists.txt
index ff813a915a6f9c2750cc6438dd91feb96c72473d..6bc244e01f3deeac4fb22350e3fc890736bba5d1 100644
--- a/dumux/porousmediumflow/CMakeLists.txt
+++ b/dumux/porousmediumflow/CMakeLists.txt
@@ -1,5 +1,6 @@
 add_subdirectory("1p")
 add_subdirectory("1p2c")
+add_subdirectory("1pnc")
 add_subdirectory("2p")
 add_subdirectory("2p1c")
 add_subdirectory("2p2c")