diff --git a/dumux/porousmediumflow/richards/model.hh b/dumux/porousmediumflow/richards/model.hh
index 54683beca98329cb1e0a5ea0d7d4f4f17e082fba..7b89488361bb4f7126fcc56ba21584b1bc5c9796 100644
--- a/dumux/porousmediumflow/richards/model.hh
+++ b/dumux/porousmediumflow/richards/model.hh
@@ -202,7 +202,7 @@ public:
 SET_TYPE_PROP(Richards, PrimaryVariableSwitch, ExtendedRichardsPrimaryVariableSwitch<TypeTag>);
 
 //! The primary variable switch for the richards model
-//SET_BOOL_PROP(Richards, ProblemUsePrimaryVariableSwitch, false);
+// SET_BOOL_PROP(Richards, ProblemUsePrimaryVariableSwitch, false);
 
 //! The spatial parameters to be employed.
 //! Use FVSpatialParams by default.
diff --git a/dumux/porousmediumflow/richards/newtonsolver.hh b/dumux/porousmediumflow/richards/newtonsolver.hh
index 7a8bc044bc2957171265d46fac0041c389764d8a..1cd9e43786968cf36d7ea6ad787d2a19444934f6 100644
--- a/dumux/porousmediumflow/richards/newtonsolver.hh
+++ b/dumux/porousmediumflow/richards/newtonsolver.hh
@@ -25,7 +25,7 @@
 #define DUMUX_RICHARDS_NEWTON_SOLVER_HH
 
 #include <dumux/common/properties.hh>
-#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/porousmediumflow/richards/privarswitchnewtonsolver.hh>
 #include <dumux/discretization/elementsolution.hh>
 
 namespace Dumux {
@@ -40,10 +40,10 @@ namespace Dumux {
  *       or from possible ModelTraits.
  */
 template <class TypeTag, class Assembler, class LinearSolver>
-class RichardsNewtonSolver : public NewtonSolver<Assembler, LinearSolver>
+class RichardsNewtonSolver : public RichardsPrivarSwitchNewtonSolver<TypeTag, Assembler, LinearSolver>
 {
     using Scalar = typename Assembler::Scalar;
-    using ParentType = NewtonSolver<Assembler, LinearSolver>;
+    using ParentType = RichardsPrivarSwitchNewtonSolver<TypeTag, Assembler, LinearSolver>;
     using SolutionVector = typename Assembler::ResidualType;
 
     using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
diff --git a/dumux/porousmediumflow/richards/primaryvariableswitch.hh b/dumux/porousmediumflow/richards/primaryvariableswitch.hh
index 3f94d659937bc6cad5cb3b5ef2a783b40cc5db22..dc28bac64af1b7897d90d85accfaec521674a730 100644
--- a/dumux/porousmediumflow/richards/primaryvariableswitch.hh
+++ b/dumux/porousmediumflow/richards/primaryvariableswitch.hh
@@ -73,6 +73,9 @@ class ExtendedRichardsPrimaryVariableSwitch
     static constexpr bool useKelvinVaporPressure
         = GET_PROP_VALUE(TypeTag, UseKelvinEquation);
 
+public:
+    using ParentType::ParentType;
+
 protected:
 
     // perform variable switch at a degree of freedom location
diff --git a/dumux/porousmediumflow/richards/privarswitchnewtonsolver.hh b/dumux/porousmediumflow/richards/privarswitchnewtonsolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7b08bf21d1927ed6d7bb8389fde6800bb16f05c7
--- /dev/null
+++ b/dumux/porousmediumflow/richards/privarswitchnewtonsolver.hh
@@ -0,0 +1,69 @@
+// -*- 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
+ * \ingroup RichardsModel
+ * \brief A Richards model newton solver.
+ */
+#ifndef DUMUX_RICHARDS_PRIVAR_SWITCH_NEWTON_SOLVER_HH
+#define DUMUX_RICHARDS_PRIVAR_SWITCH_NEWTON_SOLVER_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/discretization/elementsolution.hh>
+
+namespace Dumux {
+
+template <class TypeTag, class Assembler, class LinearSolver, bool enableWaterDiffusionInAir>
+class RichardsPrivarSwitchNewtonSolverImplementation;
+/*!
+ * \ingroup RichardsModel
+ * \brief A base for the richards newton solver which derives from the right base newton solver.
+  */
+template <class TypeTag, class Assembler, class LinearSolver>
+using RichardsPrivarSwitchNewtonSolver = RichardsPrivarSwitchNewtonSolverImplementation <TypeTag, Assembler, LinearSolver, GET_PROP_VALUE(TypeTag, EnableWaterDiffusionInAir)>;
+
+/*!
+ * \ingroup RichardsModel
+ * \brief the case without a primary variables switch
+  */
+template <class TypeTag, class Assembler, class LinearSolver>
+class RichardsPrivarSwitchNewtonSolverImplementation<TypeTag, Assembler, LinearSolver, false> : public NewtonSolver<Assembler, LinearSolver>
+{
+    using ParentType = NewtonSolver<Assembler, LinearSolver>;
+public:
+    using ParentType::ParentType;
+};
+/*!
+ * \ingroup RichardsModel
+ * \brief the case with switchable primary variables
+  */
+template <class TypeTag, class Assembler, class LinearSolver>
+class RichardsPrivarSwitchNewtonSolverImplementation<TypeTag, Assembler, LinearSolver, true> : public PriVarSwitchNewtonSolver<Assembler, LinearSolver, typename GET_PROP_TYPE(TypeTag, PrimaryVariableSwitch)>
+{
+    using PrimaryVariableSwitch = typename GET_PROP_TYPE(TypeTag, PrimaryVariableSwitch);
+    using ParentType = PriVarSwitchNewtonSolver<Assembler, LinearSolver, PrimaryVariableSwitch>;
+public:
+    using ParentType::ParentType;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/richards/implicit/CMakeLists.txt b/test/porousmediumflow/richards/implicit/CMakeLists.txt
index 8f664a66270d74eae4479e918a7407dc4774af77..34d608c47193abd7f6bf8bb4ebf5fa15d3e01be4 100644
--- a/test/porousmediumflow/richards/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/richards/implicit/CMakeLists.txt
@@ -88,6 +88,27 @@ dune_add_test(SOURCES test_richardsniconduction_fv.cc
                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_ccrichardsniconduction test_richardsniconduction.input -Problem.Name test_ccrichardsniconduction"
                       )
 
+dune_add_test(SOURCES test_richardsnievaporation_fv.cc
+              NAME test_richardsnievaporation_tpfa
+              COMPILE_DEFINITIONS TYPETAG=RichardsNIEvaporationCCTypeTag
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS --script fuzzy
+                       --files ${CMAKE_SOURCE_DIR}/test/references/test_ccrichardsevaporation-reference.vtu
+                               ${CMAKE_CURRENT_BINARY_DIR}/test_ccrichardsnievaporation-00043.vtu
+                       --command "${CMAKE_CURRENT_BINARY_DIR}/test_richardsnievaporation_tpfa test_richardsnievaporation.input -Problem.Name test_ccrichardsnievaporation"
+                      )
+
+dune_add_test(SOURCES test_richardsnievaporation_fv.cc
+              NAME test_richardsnievaporation_box
+              COMPILE_DEFINITIONS TYPETAG=RichardsNIEvaporationBoxTypeTag
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS --script fuzzy
+                       --files ${CMAKE_SOURCE_DIR}/test/references/test_boxrichardsevaporation-reference.vtu
+                               ${CMAKE_CURRENT_BINARY_DIR}/test_boxrichardsnievaporation-00043.vtu
+                       --command "${CMAKE_CURRENT_BINARY_DIR}/test_richardsnievaporation_box test_richardsnievaporation.input -Problem.Name test_boxrichardsnievaporation"
+                      )
+
+
 #install sources
 install(FILES
 richardsanalyticalproblem.hh
@@ -96,9 +117,11 @@ richardslensproblem.hh
 richardslensspatialparams.hh
 richardsniconductionproblem.hh
 richardsniconvectionproblem.hh
+richardsnievaporationproblem.hh
 richardsnispatialparams.hh
 test_richardslens_fv.cc
 test_richardsniconduction_fv.cc
 test_richardsniconvection_fv.cc
+test_richardsnievaporation_fv.cc
 test_ccrichardsanalytical.cc
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/implicit/richards)
diff --git a/test/porousmediumflow/richards/implicit/richardsnievaporationproblem.hh b/test/porousmediumflow/richards/implicit/richardsnievaporationproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1ccf6ac5802862d6515420c48dc04d14e56d486a
--- /dev/null
+++ b/test/porousmediumflow/richards/implicit/richardsnievaporationproblem.hh
@@ -0,0 +1,275 @@
+// -*- 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
+ * \ingroup RichardsTests
+ * \brief Test for the extended richards problem:
+ * The simulation domain is a tube a constant evaporation rate is set at the top and the soil gradually dries out.
+ */
+#ifndef DUMUX_RICHARDS_EVAPORATION_PROBLEM_HH
+#define DUMUX_RICHARDS_EVAPORATION_PROBLEM_HH
+
+#include <math.h>
+
+#include <dumux/discretization/elementsolution.hh>
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/discretization/box/properties.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/porousmediumflow/richards/model.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
+#include <dumux/material/fluidsystems/h2on2.hh>
+#include "richardsnispatialparams.hh"
+
+namespace Dumux {
+
+/**
+ * \ingroup RichardsTests
+ * \brief Test for the RichardsModel in combination with the NI model for an evaporation.
+ */
+template <class TypeTag>
+class RichardsNIEvaporationProblem;
+
+namespace Properties {
+NEW_TYPE_TAG(RichardsNIEvaporationTypeTag, INHERITS_FROM(RichardsNI, RichardsNISpatialParams));
+NEW_TYPE_TAG(RichardsNIEvaporationBoxTypeTag, INHERITS_FROM(BoxModel, RichardsNIEvaporationTypeTag));
+NEW_TYPE_TAG(RichardsNIEvaporationCCTypeTag, INHERITS_FROM(CCTpfaModel, RichardsNIEvaporationTypeTag));
+
+// Set the grid type
+SET_TYPE_PROP(RichardsNIEvaporationTypeTag, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(RichardsNIEvaporationTypeTag, Problem,
+              RichardsNIEvaporationProblem<TypeTag>);
+
+// Set the fluid system
+SET_TYPE_PROP(RichardsNIEvaporationTypeTag, FluidSystem, FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), false>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(RichardsNIEvaporationTypeTag,
+              SpatialParams,
+              RichardsNISpatialParams<TypeTag>);
+
+SET_BOOL_PROP(RichardsNIEvaporationTypeTag, EnableWaterDiffusionInAir, true);
+
+} // end namespace Dumux
+
+/*!
+ * \ingroup RichardsTests
+ *
+ * \brief Test for the RichardsModel in combination with the NI model for evaporation
+
+ * The result of the analytical solution is written into the vtu files.
+ * This problem uses the \ref RichardsModel and \ref NIModel model.
+ *
+ * To run the simulation execute the following line in shell: <br>
+ * <tt>./test_boxRichardsnievaporation -ParameterFile ./test_boxRichardsnievaporation.input</tt> or <br>
+ * <tt>./test_ccRichardsnievaporation -ParameterFile ./test_ccRichardsnievaporation.input</tt>
+ */
+template <class TypeTag>
+class RichardsNIEvaporationProblem : public PorousMediumFlowProblem<TypeTag>
+{
+    using ParentType = PorousMediumFlowProblem<TypeTag>;
+
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FVElementGeometry = typename FVGridGeometry::LocalView;
+    using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
+    using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables);
+    using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
+    using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+    using IapwsH2O = H2O<Scalar>;
+
+    // copy some indices for convenience
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    enum { dimWorld = GridView::dimensionworld };
+
+    enum {
+        pressureIdx = Indices::pressureIdx,
+        conti0EqIdx = Indices::conti0EqIdx,
+        temperatureIdx = Indices::temperatureIdx,
+        energyEqIdx = Indices::energyEqIdx
+    };
+
+    using Element = typename GridView::template Codim<0>::Entity;
+    using GlobalPosition = Dune::FieldVector<typename GridView::ctype, dimWorld>;
+
+public:
+    RichardsNIEvaporationProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry)
+    {
+        //initialize fluid system
+        FluidSystem::init();
+
+        name_ =  getParam<std::string>("Problem.Name");
+        pressure_ = 9.9e4;
+        temperatureInitial_ = 291;
+    }
+
+    /*!
+     * \name Problem parameters
+     */
+    // \{
+
+    /*!
+     * \brief The problem name.
+     *
+     * This is used as a prefix for files generated by the simulation.
+     */
+    const std::string& name() const
+    {
+        return name_;
+    }
+
+    // \}
+
+    /*!
+     * \name Boundary conditions
+     */
+    // \{
+
+    /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary segment.
+     *
+     * \param globalPos The position for which the boundary type is set
+     */
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+    {
+        BoundaryTypes values;
+        if(globalPos[1] < eps_)
+        {
+            values.setAllDirichlet();
+        }
+        else
+        {
+            values.setAllNeumann();
+        }
+        return values;
+    }
+
+    /*!
+     * \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
+     *
+     * For this method, the \a values parameter stores primary variables.
+     */
+    PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+    {
+        return initial_(globalPos);
+    }
+
+    /*!
+     * \brief Evaluate the boundary conditions for a neumann
+     *        boundary segment.
+     *
+     * \param element The finite element
+     * \param fvGeometry The finite-volume geometry in the box scheme
+     * \param elemVolVars The element volume variables
+     * \param scvf The subcontrolvolume face
+     *  Negative values mean influx.
+     */
+    NumEqVector neumann(const Element &element,
+                        const FVElementGeometry& fvGeometry,
+                        const ElementVolumeVariables& elemVolVars,
+                        const SubControlVolumeFace& scvf) const
+    {
+        NumEqVector values(0.0);
+        const auto globalPos = scvf.ipGlobal();
+        const auto& volVars = elemVolVars[scvf.insideScvIdx()];
+        Scalar boundaryLayerThickness = 0.0016;
+
+        if(globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_)
+        {
+             values[conti0EqIdx] = 1e-3;
+             values[energyEqIdx] = FluidSystem::enthalpy( volVars.fluidState(),Indices::nPhaseIdx) * values[conti0EqIdx];
+             values[energyEqIdx] += FluidSystem::thermalConductivity(volVars.fluidState(), Indices::nPhaseIdx)
+                                    * (volVars.temperature() - temperatureInitial_)/boundaryLayerThickness;
+        }
+        return values;
+    }
+
+    // \}
+
+    /*!
+     * \name Volume terms
+     */
+    // \{
+
+
+    /*!
+     * \brief Returns the reference pressure [Pa] of the non-wetting
+     *        fluid phase within a finite volume
+     *
+     * This problem assumes a constant reference pressure of 1 bar.
+     *
+     * \param element The DUNE Codim<0> entity which intersects with
+     *                the finite volume in question
+     * \param fvGeometry The finite volume geometry of the element
+     * \param scvIdx The sub control volume index inside the finite
+     *               volume geometry
+     */
+    Scalar nonWettingReferencePressure() const
+    { return 1e5; };
+
+    /*!
+     * \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.
+     */
+   PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+    {
+        return initial_(globalPos);
+    }
+
+    // \}
+
+private:
+    PrimaryVariables initial_(const GlobalPosition &globalPos) const
+    {
+        PrimaryVariables priVars(0.0);
+        priVars.setState(Indices::bothPhases);
+        priVars[pressureIdx] = pressure_; // initial condition for the pressure
+        priVars[temperatureIdx] = temperatureInitial_;
+        return priVars;
+    }
+
+    Scalar temperatureInitial_;
+    Scalar pressure_;
+    static constexpr Scalar eps_ = 1e-6;
+    std::string name_;
+};
+
+} //end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/richards/implicit/richardsnispatialparams.hh b/test/porousmediumflow/richards/implicit/richardsnispatialparams.hh
index 06f1bf117366b34acc2faf9483d253351c1cdee1..54a85884cf8e5414c8961d8063424318f3ba4ce0 100644
--- a/test/porousmediumflow/richards/implicit/richardsnispatialparams.hh
+++ b/test/porousmediumflow/richards/implicit/richardsnispatialparams.hh
@@ -30,8 +30,7 @@
 #include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
 #include <dumux/material/spatialparams/fv.hh>
 
-namespace Dumux
-{
+namespace Dumux {
 
 /*!
  * \ingroup RichardsTests
@@ -42,8 +41,7 @@ namespace Dumux
 template<class TypeTag>
 class RichardsNISpatialParams;
 
-namespace Properties
-{
+namespace Properties {
 // The spatial parameters TypeTag
 NEW_TYPE_TAG(RichardsNISpatialParams);
 
@@ -62,8 +60,7 @@ public:
     // define the material law parameterized by absolute saturations
     using type = EffToAbsLaw<EffectiveLaw>;
 };
-}
-
+} // end namespace Properties
 
 template<class TypeTag>
 class RichardsNISpatialParams : public FVSpatialParams<TypeTag>
@@ -107,10 +104,6 @@ public:
         materialParams_.setVgn(4.7);
     }
 
-    ~RichardsNISpatialParams()
-    {}
-
-
     /*!
      * \brief Define the intrinsic permeability \f$\mathrm{[m^2]}\f$.
      *
@@ -183,7 +176,6 @@ public:
         return lambdaSolid_;
     }
 
-
 private:
 
     MaterialLawParams materialParams_;
@@ -192,6 +184,6 @@ private:
     Scalar lambdaSolid_;
 };
 
-}
+} // end namespace Dumux
 
 #endif
diff --git a/test/porousmediumflow/richards/implicit/test_richardsnievaporation.input b/test/porousmediumflow/richards/implicit/test_richardsnievaporation.input
new file mode 100644
index 0000000000000000000000000000000000000000..5c8a107941d04b4ff22a817353c21279d746a65d
--- /dev/null
+++ b/test/porousmediumflow/richards/implicit/test_richardsnievaporation.input
@@ -0,0 +1,17 @@
+[TimeLoop]
+DtInitial = 1 # [s]
+TEnd = 3e4 # [s]
+MaxTimeStepSize = 1e3 # [s]
+
+[Grid]
+UpperRight = 10 20
+Cells = 5 20
+
+[Problem]
+Name = test_richardsevaporation # name passed to the output routines
+OutputInterval = 5 # every 5th timestep an output file is written
+EnableGravity = 0 # disable gravity
+UsePrimaryVariableSwitch = true
+
+[Newton]
+EnableChop = false  # chop for better convergence
diff --git a/test/porousmediumflow/richards/implicit/test_richardsnievaporation_fv.cc b/test/porousmediumflow/richards/implicit/test_richardsnievaporation_fv.cc
new file mode 100644
index 0000000000000000000000000000000000000000..bf53af86c7a712d480f61f29bda2c1684696206a
--- /dev/null
+++ b/test/porousmediumflow/richards/implicit/test_richardsnievaporation_fv.cc
@@ -0,0 +1,194 @@
+// -*- 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 Richards box model.
+ */
+#include <config.h>
+
+#include "richardsnievaporationproblem.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/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/porousmediumflow/richards/newtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+
+////////////////////////
+// the main function
+////////////////////////
+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);
+
+    // parse command line arguments and input file
+    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);
+    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+    auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+    // intialize the vtk output module
+    using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+    VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+    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 = Dumux::ILU0BiCGSTABBackend;
+    auto linearSolver = std::make_shared<LinearSolver>();
+
+    // the non-linear solver
+    using NewtonSolver = Dumux::RichardsNewtonSolver<TypeTag, Assembler, LinearSolver>;
+    NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+    // time loop
+    timeLoop->start(); do
+    {
+        // set previous solution for storage evaluations
+        assembler->setPreviousSolution(xOld);
+
+        // solve the non-linear system with time step control
+        nonLinearSolver.solve(x, *timeLoop);
+
+        // 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(nonLinearSolver.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)
+    {
+        Parameters::print();
+        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/references/test_boxrichardsevaporation-reference.vtu b/test/references/test_boxrichardsevaporation-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..ac8cb7b8b07f22e4a35ed3cecc40e78c942001ad
--- /dev/null
+++ b/test/references/test_boxrichardsevaporation-reference.vtu
@@ -0,0 +1,287 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+  <UnstructuredGrid>
+    <Piece NumberOfCells="100" NumberOfPoints="126">
+      <PointData Scalars="Sw">
+        <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571695 0.0571695 0.0571695 0.0571695 0.0571695 0.0571695
+          0 0 0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.94283 0.94283 0.94283 0.94283 0.94283 0.94283
+          1 1 1 1 1 1
+        </DataArray>
+        <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 98999.9 98999.9 98999.9 98999.9 98999.9 98999.9
+          97470.4 97470.4 97470.4 97470.4 97470.4 97470.4
+        </DataArray>
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000
+        </DataArray>
+        <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000.07 1000.07 1000.07 1000.07 1000.07 1000.07
+          2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
+        </DataArray>
+        <DataArray type="Float32" Name="density" NumberOfComponents="1" format="ascii">
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.632 998.632 998.632 998.632 998.632 998.632
+          999.788 999.788 999.788 999.788 999.788 999.788
+        </DataArray>
+        <DataArray type="Float32" Name="mobility" NumberOfComponents="1" format="ascii">
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207045 0.000207045 0.000207045 0.000207045 0.000207045 0.000207045 0.000206541 0.000206541 0.000206541 0.000206541 0.000206541 0.000206541
+          0 0 0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="kr" NumberOfComponents="1" format="ascii">
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18843e-07 2.18843e-07 2.18843e-07 2.18843e-07 2.18843e-07 2.18843e-07 2.18606e-07 2.18606e-07 2.18606e-07 2.18606e-07 2.18606e-07 2.18606e-07
+          0 0 0 0 0 0
+        </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
+          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
+          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
+          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
+          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
+          0.4 0.4 0.4 0.4 0.4 0.4
+        </DataArray>
+        <DataArray type="Float32" Name="x^w_air" NumberOfComponents="1" format="ascii">
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204739 0.0204739 0.0204739 0.0204739 0.0204739 0.0204739 0.0203936 0.0203936 0.0203936 0.0203936 0.0203936 0.0203936
+          -120.791 -120.791 -120.791 -120.791 -120.791 -120.791
+        </DataArray>
+        <DataArray type="Float32" Name="water content" NumberOfComponents="1" format="ascii">
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228678 0.0228678 0.0228678 0.0228678 0.0228678 0.0228678
+          0 0 0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="phasePresence" NumberOfComponents="1" format="ascii">
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          2 2 2 2 2 2
+        </DataArray>
+        <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 290.943 290.943 290.943 290.943 290.943 290.943
+          281.758 281.758 281.758 281.758 281.758 281.758
+        </DataArray>
+      </PointData>
+      <CellData Scalars="process rank">
+        <DataArray type="Float32" Name="process rank" 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
+        </DataArray>
+      </CellData>
+      <Points>
+        <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+          0 0 0 2 0 0 0 1 0 2 1 0
+          4 0 0 4 1 0 6 0 0 6 1 0
+          8 0 0 8 1 0 10 0 0 10 1 0
+          0 2 0 2 2 0 4 2 0 6 2 0
+          8 2 0 10 2 0 0 3 0 2 3 0
+          4 3 0 6 3 0 8 3 0 10 3 0
+          0 4 0 2 4 0 4 4 0 6 4 0
+          8 4 0 10 4 0 0 5 0 2 5 0
+          4 5 0 6 5 0 8 5 0 10 5 0
+          0 6 0 2 6 0 4 6 0 6 6 0
+          8 6 0 10 6 0 0 7 0 2 7 0
+          4 7 0 6 7 0 8 7 0 10 7 0
+          0 8 0 2 8 0 4 8 0 6 8 0
+          8 8 0 10 8 0 0 9 0 2 9 0
+          4 9 0 6 9 0 8 9 0 10 9 0
+          0 10 0 2 10 0 4 10 0 6 10 0
+          8 10 0 10 10 0 0 11 0 2 11 0
+          4 11 0 6 11 0 8 11 0 10 11 0
+          0 12 0 2 12 0 4 12 0 6 12 0
+          8 12 0 10 12 0 0 13 0 2 13 0
+          4 13 0 6 13 0 8 13 0 10 13 0
+          0 14 0 2 14 0 4 14 0 6 14 0
+          8 14 0 10 14 0 0 15 0 2 15 0
+          4 15 0 6 15 0 8 15 0 10 15 0
+          0 16 0 2 16 0 4 16 0 6 16 0
+          8 16 0 10 16 0 0 17 0 2 17 0
+          4 17 0 6 17 0 8 17 0 10 17 0
+          0 18 0 2 18 0 4 18 0 6 18 0
+          8 18 0 10 18 0 0 19 0 2 19 0
+          4 19 0 6 19 0 8 19 0 10 19 0
+          0 20 0 2 20 0 4 20 0 6 20 0
+          8 20 0 10 20 0
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+          0 1 3 2 1 4 5 3 4 6 7 5
+          6 8 9 7 8 10 11 9 2 3 13 12
+          3 5 14 13 5 7 15 14 7 9 16 15
+          9 11 17 16 12 13 19 18 13 14 20 19
+          14 15 21 20 15 16 22 21 16 17 23 22
+          18 19 25 24 19 20 26 25 20 21 27 26
+          21 22 28 27 22 23 29 28 24 25 31 30
+          25 26 32 31 26 27 33 32 27 28 34 33
+          28 29 35 34 30 31 37 36 31 32 38 37
+          32 33 39 38 33 34 40 39 34 35 41 40
+          36 37 43 42 37 38 44 43 38 39 45 44
+          39 40 46 45 40 41 47 46 42 43 49 48
+          43 44 50 49 44 45 51 50 45 46 52 51
+          46 47 53 52 48 49 55 54 49 50 56 55
+          50 51 57 56 51 52 58 57 52 53 59 58
+          54 55 61 60 55 56 62 61 56 57 63 62
+          57 58 64 63 58 59 65 64 60 61 67 66
+          61 62 68 67 62 63 69 68 63 64 70 69
+          64 65 71 70 66 67 73 72 67 68 74 73
+          68 69 75 74 69 70 76 75 70 71 77 76
+          72 73 79 78 73 74 80 79 74 75 81 80
+          75 76 82 81 76 77 83 82 78 79 85 84
+          79 80 86 85 80 81 87 86 81 82 88 87
+          82 83 89 88 84 85 91 90 85 86 92 91
+          86 87 93 92 87 88 94 93 88 89 95 94
+          90 91 97 96 91 92 98 97 92 93 99 98
+          93 94 100 99 94 95 101 100 96 97 103 102
+          97 98 104 103 98 99 105 104 99 100 106 105
+          100 101 107 106 102 103 109 108 103 104 110 109
+          104 105 111 110 105 106 112 111 106 107 113 112
+          108 109 115 114 109 110 116 115 110 111 117 116
+          111 112 118 117 112 113 119 118 114 115 121 120
+          115 116 122 121 116 117 123 122 117 118 124 123
+          118 119 125 124
+        </DataArray>
+        <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+          4 8 12 16 20 24 28 32 36 40 44 48
+          52 56 60 64 68 72 76 80 84 88 92 96
+          100 104 108 112 116 120 124 128 132 136 140 144
+          148 152 156 160 164 168 172 176 180 184 188 192
+          196 200 204 208 212 216 220 224 228 232 236 240
+          244 248 252 256 260 264 268 272 276 280 284 288
+          292 296 300 304 308 312 316 320 324 328 332 336
+          340 344 348 352 356 360 364 368 372 376 380 384
+          388 392 396 400
+        </DataArray>
+        <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>
diff --git a/test/references/test_ccrichardsevaporation-reference.vtu b/test/references/test_ccrichardsevaporation-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..2047c2d055404a8b69ebfe8e2efe5498f41cfe09
--- /dev/null
+++ b/test/references/test_ccrichardsevaporation-reference.vtu
@@ -0,0 +1,259 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+  <UnstructuredGrid>
+    <Piece NumberOfCells="100" NumberOfPoints="126">
+      <CellData Scalars="Sw">
+        <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721
+          0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571721 0.0571703 0.0571703 0.0571703 0.0571703 0.0571703 0
+          0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.942828
+          0.942828 0.942828 0.942828 0.942828 0.942828 0.942828 0.94283 0.94283 0.94283 0.94283 0.94283 1
+          1 1 1 1
+        </DataArray>
+        <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000
+          99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 99000 97470.4
+          97470.4 97470.4 97470.4 97470.4
+        </DataArray>
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000
+        </DataArray>
+        <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
+          1000 1000 1000 1000 1000 1000 1000.05 1000.05 1000.05 1000.05 1000.05 2529.62
+          2529.62 2529.62 2529.62 2529.62
+        </DataArray>
+        <DataArray type="Float32" Name="density" NumberOfComponents="1" format="ascii">
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621 998.621
+          998.621 998.621 998.621 998.621 998.621 998.621 998.627 998.627 998.627 998.627 998.627 999.369
+          999.369 999.369 999.369 999.369
+        </DataArray>
+        <DataArray type="Float32" Name="mobility" NumberOfComponents="1" format="ascii">
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047 0.000207047
+          0.000207047 0.000207046 0.000207046 0.000207046 0.000207046 0.000207046 0.000206733 0.000206733 0.000206733 0.000206733 0.000206733 0
+          0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="kr" NumberOfComponents="1" format="ascii">
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07 2.18844e-07
+          2.18844e-07 2.18843e-07 2.18843e-07 2.18843e-07 2.18843e-07 2.18843e-07 2.18673e-07 2.18673e-07 2.18673e-07 2.18673e-07 2.18673e-07 0
+          0 0 0 0
+        </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
+          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
+          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
+          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
+          0.4 0.4 0.4 0.4
+        </DataArray>
+        <DataArray type="Float32" Name="x^w_air" NumberOfComponents="1" format="ascii">
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743 0.0204743
+          0.0204743 0.0204741 0.0204741 0.0204741 0.0204741 0.0204741 0.0204304 0.0204304 0.0204304 0.0204304 0.0204304 -23.6214
+          -23.6214 -23.6214 -23.6214 -23.6214
+        </DataArray>
+        <DataArray type="Float32" Name="water content" NumberOfComponents="1" format="ascii">
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688
+          0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228688 0.0228681 0.0228681 0.0228681 0.0228681 0.0228681 0
+          0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="phasePresence" NumberOfComponents="1" format="ascii">
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 3
+          3 3 3 3 3 3 3 3 3 3 3 2
+          2 2 2 2
+        </DataArray>
+        <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 290.969 290.969 290.969 290.969 290.969 286.074
+          286.074 286.074 286.074 286.074
+        </DataArray>
+        <DataArray type="Float32" Name="process rank" 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
+        </DataArray>
+      </CellData>
+      <Points>
+        <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+          0 0 0 2 0 0 0 1 0 2 1 0
+          4 0 0 4 1 0 6 0 0 6 1 0
+          8 0 0 8 1 0 10 0 0 10 1 0
+          0 2 0 2 2 0 4 2 0 6 2 0
+          8 2 0 10 2 0 0 3 0 2 3 0
+          4 3 0 6 3 0 8 3 0 10 3 0
+          0 4 0 2 4 0 4 4 0 6 4 0
+          8 4 0 10 4 0 0 5 0 2 5 0
+          4 5 0 6 5 0 8 5 0 10 5 0
+          0 6 0 2 6 0 4 6 0 6 6 0
+          8 6 0 10 6 0 0 7 0 2 7 0
+          4 7 0 6 7 0 8 7 0 10 7 0
+          0 8 0 2 8 0 4 8 0 6 8 0
+          8 8 0 10 8 0 0 9 0 2 9 0
+          4 9 0 6 9 0 8 9 0 10 9 0
+          0 10 0 2 10 0 4 10 0 6 10 0
+          8 10 0 10 10 0 0 11 0 2 11 0
+          4 11 0 6 11 0 8 11 0 10 11 0
+          0 12 0 2 12 0 4 12 0 6 12 0
+          8 12 0 10 12 0 0 13 0 2 13 0
+          4 13 0 6 13 0 8 13 0 10 13 0
+          0 14 0 2 14 0 4 14 0 6 14 0
+          8 14 0 10 14 0 0 15 0 2 15 0
+          4 15 0 6 15 0 8 15 0 10 15 0
+          0 16 0 2 16 0 4 16 0 6 16 0
+          8 16 0 10 16 0 0 17 0 2 17 0
+          4 17 0 6 17 0 8 17 0 10 17 0
+          0 18 0 2 18 0 4 18 0 6 18 0
+          8 18 0 10 18 0 0 19 0 2 19 0
+          4 19 0 6 19 0 8 19 0 10 19 0
+          0 20 0 2 20 0 4 20 0 6 20 0
+          8 20 0 10 20 0
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+          0 1 3 2 1 4 5 3 4 6 7 5
+          6 8 9 7 8 10 11 9 2 3 13 12
+          3 5 14 13 5 7 15 14 7 9 16 15
+          9 11 17 16 12 13 19 18 13 14 20 19
+          14 15 21 20 15 16 22 21 16 17 23 22
+          18 19 25 24 19 20 26 25 20 21 27 26
+          21 22 28 27 22 23 29 28 24 25 31 30
+          25 26 32 31 26 27 33 32 27 28 34 33
+          28 29 35 34 30 31 37 36 31 32 38 37
+          32 33 39 38 33 34 40 39 34 35 41 40
+          36 37 43 42 37 38 44 43 38 39 45 44
+          39 40 46 45 40 41 47 46 42 43 49 48
+          43 44 50 49 44 45 51 50 45 46 52 51
+          46 47 53 52 48 49 55 54 49 50 56 55
+          50 51 57 56 51 52 58 57 52 53 59 58
+          54 55 61 60 55 56 62 61 56 57 63 62
+          57 58 64 63 58 59 65 64 60 61 67 66
+          61 62 68 67 62 63 69 68 63 64 70 69
+          64 65 71 70 66 67 73 72 67 68 74 73
+          68 69 75 74 69 70 76 75 70 71 77 76
+          72 73 79 78 73 74 80 79 74 75 81 80
+          75 76 82 81 76 77 83 82 78 79 85 84
+          79 80 86 85 80 81 87 86 81 82 88 87
+          82 83 89 88 84 85 91 90 85 86 92 91
+          86 87 93 92 87 88 94 93 88 89 95 94
+          90 91 97 96 91 92 98 97 92 93 99 98
+          93 94 100 99 94 95 101 100 96 97 103 102
+          97 98 104 103 98 99 105 104 99 100 106 105
+          100 101 107 106 102 103 109 108 103 104 110 109
+          104 105 111 110 105 106 112 111 106 107 113 112
+          108 109 115 114 109 110 116 115 110 111 117 116
+          111 112 118 117 112 113 119 118 114 115 121 120
+          115 116 122 121 116 117 123 122 117 118 124 123
+          118 119 125 124
+        </DataArray>
+        <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+          4 8 12 16 20 24 28 32 36 40 44 48
+          52 56 60 64 68 72 76 80 84 88 92 96
+          100 104 108 112 116 120 124 128 132 136 140 144
+          148 152 156 160 164 168 172 176 180 184 188 192
+          196 200 204 208 212 216 220 224 228 232 236 240
+          244 248 252 256 260 264 268 272 276 280 284 288
+          292 296 300 304 308 312 316 320 324 328 332 336
+          340 344 348 352 356 360 364 368 372 376 380 384
+          388 392 396 400
+        </DataArray>
+        <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>