diff --git a/config.h.cmake b/config.h.cmake
index aee39426a82ac0849bad7fdbe7fd05d505298a94..859fa68708b447269b138d5ee9e3f825a91c5db6 100644
--- a/config.h.cmake
+++ b/config.h.cmake
@@ -40,6 +40,9 @@
 /* Define to the revision of dumux */
 #define DUMUX_VERSION_REVISION ${DUMUX_VERSION_REVISION}
 
+/* Define the path to dumux */
+#define DUMUX_SOURCE_DIR "${CMAKE_SOURCE_DIR}"
+
 /* Define to 1 if Valgrind was found */
 #cmakedefine HAVE_VALGRIND 1
 
@@ -55,6 +58,12 @@
 /* Define path to gstat executable */
 #cmakedefine GSTAT_EXECUTABLE "@GSTAT_EXECUTABLE@"
 
+/* Defines whether pvpython has been found */
+#cmakedefine HAVE_PVPYTHON 1
+
+/* Define the path to pvpython */
+#define PVPYTHON_EXECUTABLE "${PVPYTHON_EXECUTABLE}"
+
 /* end dumux
    Everything below here will be overwritten
 */
diff --git a/doc/doxygen/modules.txt b/doc/doxygen/modules.txt
index 956ecf1ad9e342541b8be24986049dc74aacf1b3..a8f3cabf3d0e81de75b9e8c3e4cc10040799e75d 100644
--- a/doc/doxygen/modules.txt
+++ b/doc/doxygen/modules.txt
@@ -170,26 +170,38 @@
 /* ***************** FreeflowModels ******************/
 /*!
  * \defgroup FreeflowModels Free Flow Models
- * \brief Single-phase Navier Stokes / Stokes model
+ * \brief Single-phase Navier-Stokes / Stokes model
  */
 	/*!
      * \ingroup FreeflowModels
-     * \defgroup NavierStokesModel NavierStokes
-     * \brief Single-phase Navier Stokes flow
+     * \defgroup NavierStokesModel Navier-Stokes
+     * \brief Single-phase Navier-Stokes flow
      * \copydetails ./freeflow/navierstokes/model.hh
      */
     /*!
      * \ingroup FreeflowModels
-     * \defgroup NavierStokesNCModel NavierStokes nc
-     * \brief Single-phase multi-component Navier Stokes flow
+     * \defgroup NavierStokesNCModel Navier-Stokes nc
+     * \brief Single-phase multi-component Navier-Stokes flow
      * \copydetails ./freeflow/navierstokesnc/model.hh
      */
     /*!
      * \ingroup FreeflowModels
      * \defgroup NavierStokesNIModel nonisothermal
-     * \brief An energy equation adaptor for isothermal Navier Stokes models
+     * \brief An energy equation adaptor for isothermal Navier-Stokes models
      * \copydetails ./freeflow/nonisothermal/model.hh
      */
+    /*!
+     * \ingroup FreeflowModels
+     * \defgroup RANSModel Reynolds-Averaged Navier-Stokes
+     * \brief Single-phase Reynolds-Averaged Navier-Stokes flow
+     * \copydetails ./freeflow/rans/model.hh
+     */
+      /*!
+      * \ingroup RANSModel
+      * \defgroup ZeroEqModel 0-Eq. Models
+      * \brief Zero-equation or algebraic turbulence models
+      * \copydetails ./freeflow/rans/model.hh
+      */
 
 /* ***************** Benchmarks and Tests ******************/
 /*!
@@ -284,17 +296,17 @@
 	/*!  
  	 * \ingroup BenchmarksAndTests
  	 * \defgroup FreeflowTests Free Flow Tests
- 	 * \brief Varios tests for single-phase Navier Stokes / Stokes tests
+ 	 * \brief Varios tests for single-phase Navier-Stokes / Stokes tests
  	 */
  	  	/*!  
  	     * \ingroup FreeflowTests
-	 	 * \defgroup NavierStokesTests Single-phase Navier Stokes tests
-	 	 * \brief Various tests using a Single-phase Navier Stokes flow. The files are listed below.
+	 	 * \defgroup NavierStokesTests Single-phase Navier-Stokes tests
+	 	 * \brief Various tests using a Single-phase Navier-Stokes flow. The files are listed below.
 	 	 */
 	   	/*!  
  	     * \ingroup FreeflowTests
-	 	 * \defgroup NavierStokesNCTests Single-phase Navier Stokes nc tests
-	 	 * \brief Various tests using a Single-phase Navier Stokes flow. The files are listed below.
+	 	 * \defgroup NavierStokesNCTests Single-phase Navier-Stokes nc tests
+	 	 * \brief Various tests using a Single-phase Navier-Stokes flow. The files are listed below.
 	 	 */
 	/*!  
  	 * \ingroup BenchmarksAndTests
diff --git a/doc/handbook/dumux-handbook.bib b/doc/handbook/dumux-handbook.bib
index e25743f6a27d9e410f2d2b71c592367f2da4c3e9..fc2d2dc277b178b9ace61eeb2b06bff36fbf3038 100644
--- a/doc/handbook/dumux-handbook.bib
+++ b/doc/handbook/dumux-handbook.bib
@@ -1576,3 +1576,64 @@ url={http://dx.doi.org/10.1007/s11242-015-0599-1}
   timestamp = {2017.05.11},
   url =       {http://dx.doi.org/10.1023/B:IJOT.0000022327.04529.f3}
 }
+
+@Article{Laufer1954a,
+  author =    {John Laufer},
+  title =     {{The structure of turbulence in fully developed pipe flow}},
+  journal =   {NACA Report},
+  year =      {1954},
+  volume =    {1174},
+  pages =     {417--434},
+  owner =     {fetzer},
+  timestamp = {2012.09.03},
+  url =       {https://ntrs.nasa.gov/search.jsp?R=19930092199}
+}
+
+@Article{Hanna1981a,
+  author =    {O. T. Hanna and O. C. Sandell and P. R. Mazet},
+  title =     {{Heat and Mass Transfer in Turbulent Flow Under Conditions of Drag Reduction}},
+  journal =   {AIChE Journal},
+  year =      {1981},
+  volume =    {27},
+  number =    {4},
+  pages =     {693--697},
+  doi =       {10.1002/aic.690270424},
+  url =       {http://dx.doi.org/10.1002/aic.690270424}
+}
+
+@Book{Oertel2012a,
+  title =     {{Prandtl - F\"uhrer durch die Str\"omungslehre: Grundlagen und Ph\"anomene}},
+  publisher = {Springer Vieweg},
+  year =      {2012},
+  author =    {Oertel, Herbert},
+  isbn =      {978-3-8348-2315.1},
+  pages =     {XII, 764},
+  address =   {Wiesbaden},
+  edition =   {13},
+  editor =    {Oertel, Herbert and Böhle, Martin},
+  doi =       {10.1007/978-3-8348-2315-1},
+  url =       {http://dx.doi.org/10.1007/978-3-8348-2315-1}
+}
+
+@Article{vanDriest1956a,
+  author =    {E. R. {van Driest}},
+  title =     {{On Turbulent Flow Near a Wall}},
+  journal =   {AIAA Journal},
+  year =      {1956},
+  volume =    {23},
+  number =    {11},
+  pages =     {1007--1011},
+  doi =       {10.2514/8.3713},
+  url =       {http://dx.doi.org/10.2514/8.3713}
+}
+
+@Article{Baldwin1978a,
+  author =    {B. S. Baldwin and H. Lomax},
+  title =     {{Thin Layer Approximation and Algebraic Model for Seperated Turbulent Flows}},
+  journal =   {ACM Trans Math Software},
+  year =      {1978},
+  volume =    {78--257},
+  pages =     {1--9},
+  doi =       {10.2514/6.1978-257},
+  url =       {http://dx.doi.org/10.2514/6.1978-257}
+}
diff --git a/dumux/common/parameters.hh b/dumux/common/parameters.hh
index d954c3e98f676aab89fece294269d1e22e1062af..e70273449b32d3ae0779bb1e7dba4e92ee150809 100644
--- a/dumux/common/parameters.hh
+++ b/dumux/common/parameters.hh
@@ -259,6 +259,10 @@ private:
 
         // parameters in the mpfa group
         params["Mpfa.Q"] = "0.0";
+
+        // parameters in the RANS group
+        params["RANS.KarmanConstant"] = "0.41";
+        params["RANS.EddyViscosityModel"] = "1";
     }
 };
 
diff --git a/dumux/freeflow/CMakeLists.txt b/dumux/freeflow/CMakeLists.txt
index dee5ba1349e17c1b5d5c328c6d0a2116b1344f51..6f831c047117ec2d6b8df38d753709f959b2768f 100644
--- a/dumux/freeflow/CMakeLists.txt
+++ b/dumux/freeflow/CMakeLists.txt
@@ -1,3 +1,4 @@
 add_subdirectory("navierstokes")
 add_subdirectory("navierstokesnc")
 add_subdirectory("nonisothermal")
+add_subdirectory("rans")
diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index 1bb579746297a77413dba03cbe9f8a64b39f665d..6df4d798452eb5ee88ce1bf96d0e484b9c343492 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -207,6 +207,7 @@ public:
                 return (upwindWeight * upstreamVelocity + (1.0 - upwindWeight) * downstreamVelocity) * insideVolVars.density();
             };
 
+
             // Get the momentum that is advectively transported and account for the flow direction.
             const Scalar momentum = selfIsUpstream ? computeMomentum(velocitySelf, velocityOpposite)
                                                    : computeMomentum(velocityOpposite, velocitySelf);
@@ -220,7 +221,7 @@ public:
         // The velocity gradient already accounts for the orientation
         // of the staggered face's outer normal vector.
         const Scalar gradV = (velocityOpposite - velocitySelf) / scvf.selfToOppositeDistance();
-        normalFlux -= insideVolVars.viscosity() * 2.0 * gradV;
+        normalFlux -= insideVolVars.effectiveViscosity() * 2.0 * gradV;
 
         // The pressure term.
         // If specified, the pressure can be normalized using the initial value on the scfv of interest.
@@ -419,7 +420,7 @@ private:
         const auto& outsideVolVars = elemVolVars[normalFace.outsideScvIdx()];
 
         // Get the averaged viscosity at the staggered face normal to the current scvf.
-        const Scalar muAvg = (insideVolVars.viscosity() + outsideVolVars.viscosity()) * 0.5;
+        const Scalar muAvg = (insideVolVars.effectiveViscosity() + outsideVolVars.effectiveViscosity()) * 0.5;
 
         // For the normal gradient, get the velocities perpendicular to the velocity at the current scvf.
         // The inner one is located at staggered face within the own element,
diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh
index e37e58a5ee71f76bb572fc4de594d1c10e2408e0..3cf6aa6ede5ae1aa7a42b1cc1875ef629281f25e 100644
--- a/dumux/freeflow/navierstokes/staggered/localresidual.hh
+++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh
@@ -183,22 +183,22 @@ public:
 
 protected:
 
-     /*!
-     * \brief Evaluate boundary conditions
-     */
-    void evalBoundary_(const Element& element,
-                       const FVElementGeometry& fvGeometry,
-                       const ElementVolumeVariables& elemVolVars,
-                       const ElementFaceVariables& elemFaceVars,
-                       const ElementBoundaryTypes& elemBcTypes,
-                       const ElementFluxVariablesCache& elemFluxVarsCache)
-    {
-        evalBoundaryForCellCenter_(element, fvGeometry, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
-        for (auto&& scvf : scvfs(fvGeometry))
-        {
-            evalBoundaryForFace_(element, fvGeometry, scvf, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
-        }
-    }
+//      /*!
+//      * \brief Evaluate boundary conditions
+//      */
+//     void evalBoundary_(const Element& element,
+//                        const FVElementGeometry& fvGeometry,
+//                        const ElementVolumeVariables& elemVolVars,
+//                        const ElementFaceVariables& elemFaceVars,
+//                        const ElementBoundaryTypes& elemBcTypes,
+//                        const ElementFluxVariablesCache& elemFluxVarsCache)
+//     {
+//         evalBoundaryForCellCenter_(element, fvGeometry, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
+//         for (auto&& scvf : scvfs(fvGeometry))
+//         {
+//             evalBoundaryForFace_(element, fvGeometry, scvf, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
+//         }
+//     }
 
      /*!
      * \brief Evaluate boundary conditions for a cell center dof
diff --git a/dumux/freeflow/navierstokes/volumevariables.hh b/dumux/freeflow/navierstokes/volumevariables.hh
index 535dde96f7a2dc73bcc4b30dac1f51a1993616a5..82290cb5f54ca67c2d769998ccb9f29d9739436c 100644
--- a/dumux/freeflow/navierstokes/volumevariables.hh
+++ b/dumux/freeflow/navierstokes/volumevariables.hh
@@ -151,12 +151,6 @@ public:
     Scalar pressure(int phaseIdx = 0) const
     { return fluidState_.pressure(defaultPhaseIdx); }
 
-    /*!
-     * \brief Return the saturation
-     */
-    Scalar saturation(int phaseIdx = 0) const
-    { return 1.0; }
-
     /*!
      * \brief Return the mass density \f$\mathrm{[kg/m^3]}\f$ of a given phase within the
      *        control volume.
@@ -189,6 +183,13 @@ public:
     Scalar viscosity(int phaseIdx = 0) const
     { return fluidState_.viscosity(defaultPhaseIdx); }
 
+    /*!
+     * \brief Return the effective dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
+     *        control volume.
+     */
+    Scalar effectiveViscosity(int phaseIdx = 0) const
+    { return viscosity(defaultPhaseIdx); }
+
     /*!
      * \brief Return the fluid state of the control volume.
      */
diff --git a/dumux/freeflow/rans/CMakeLists.txt b/dumux/freeflow/rans/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b1ccda770bea551460d825be7284ae307a75d379
--- /dev/null
+++ b/dumux/freeflow/rans/CMakeLists.txt
@@ -0,0 +1,9 @@
+add_subdirectory("zeroeq")
+
+#install headers
+install(FILES
+model.hh
+problem.hh
+volumevariables.hh
+vtkoutputfields.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/rans)
diff --git a/dumux/freeflow/rans/model.hh b/dumux/freeflow/rans/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f10fd67c78538ddb1e07d26aa8396a1204165fc6
--- /dev/null
+++ b/dumux/freeflow/rans/model.hh
@@ -0,0 +1,90 @@
+// -*- 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 RANSModel
+ *
+ * \brief A single-phase, isothermal Reynolds-Averaged Navier-Stokes model
+ *
+ * This model implements a single-phase, isothermal Reynolds-Averaged
+ * Navier-Stokes model, solving the <B> momentum balance equation </B>
+ * \f[
+ \frac{\partial (\varrho \textbf{v})}{\partial t} + \nabla \cdot (\varrho \textbf{v} \textbf{v}^{\textup{T}}) = \nabla \cdot (\mu_\textrm{eff} (\nabla \textbf{v} + \nabla \textbf{v}^{\textup{T}}))
+   - \nabla p + \varrho \textbf{g} - \textbf{f}
+ * \f]
+ * The effective viscosity is composed of the fluid and the eddy viscosity:
+ * \f[
+ *    \mu_\textrm{eff} = \mu + \mu_\textrm{t}
+ * \f].
+ */
+
+#ifndef DUMUX_RANS_MODEL_HH
+#define DUMUX_RANS_MODEL_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/freeflow/properties.hh>
+#include <dumux/freeflow/navierstokes/model.hh>
+#include <dumux/freeflow/nonisothermal/model.hh>
+#include <dumux/material/fluidstates/immiscible.hh>
+
+#include "volumevariables.hh"
+#include "vtkoutputfields.hh"
+
+namespace Dumux
+{
+
+// \{
+///////////////////////////////////////////////////////////////////////////
+// properties for the single-phase Reynolds-Averaged Navier-Stokes model
+///////////////////////////////////////////////////////////////////////////
+namespace Properties {
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the single-phase, isothermal Reynolds-Averaged Navier-Stokes model
+NEW_TYPE_TAG(RANS, INHERITS_FROM(NavierStokes));
+
+//! The type tag for the corresponding non-isothermal model
+NEW_TYPE_TAG(RANSNI, INHERITS_FROM(RANS, NavierStokesNonIsothermal));
+
+///////////////////////////////////////////////////////////////////////////
+// default property values for the isothermal single phase model
+///////////////////////////////////////////////////////////////////////////
+SET_BOOL_PROP(RANS, EnableInertiaTerms, true); //!< Explicitly force the consideration of inertia terms by default
+
+//! The volume variables
+SET_TYPE_PROP(RANS, VolumeVariables, RANSVolumeVariables<TypeTag>);
+
+//! The specific vtk output fields
+SET_PROP(RANS, VtkOutputFields)
+{
+private:
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+public:
+     using type = RANSVtkOutputFields<FVGridGeometry>;
+};
+// \}
+}
+
+} // end namespace
+
+#endif // DUMUX_RANS_MODEL_HH
diff --git a/dumux/freeflow/rans/problem.hh b/dumux/freeflow/rans/problem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3eb655df03dffb66f5a5a5a1579d62ccbfe77ac0
--- /dev/null
+++ b/dumux/freeflow/rans/problem.hh
@@ -0,0 +1,339 @@
+// -*- 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 RANSModel
+ * \copydoc Dumux::RANSProblem
+ */
+#ifndef DUMUX_RANS_PROBLEM_HH
+#define DUMUX_RANS_PROBLEM_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/staggeredfvproblem.hh>
+#include <dumux/discretization/localview.hh>
+#include <dumux/discretization/staggered/elementsolution.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/freeflow/navierstokes/problem.hh>
+
+#include "model.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup RANSModel
+ * \brief Reynolds-Averaged Navier-Stokes problem base class.
+ *
+ * This implements some base functionality for RANS models.
+ * Especially vectors containing all wall-relevant properties, which are accessed
+ * by the volumevariables.
+ * \todo inherit all functions (especially gravity and temperature from Navier-Stokes)
+ * This implements gravity (if desired) and a function returning the temperature.
+ * Includes a specialized method used only by the staggered grid discretization.
+ */
+template<class TypeTag>
+class RANSProblem : public NavierStokesProblem<TypeTag>
+{
+    using ParentType = NavierStokesProblem<TypeTag>;
+    using Implementation = typename GET_PROP_TYPE(TypeTag, Problem);
+
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Grid = typename GridView::Grid;
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+    enum {
+        dim = Grid::dimension,
+      };
+    using GlobalPosition = Dune::FieldVector<Scalar, dim>;
+    using DimVector = Dune::FieldVector<Scalar, dim>;
+    using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
+
+    enum {
+        massBalanceIdx = Indices::massBalanceIdx,
+        momentumBalanceIdx = Indices::momentumBalanceIdx
+    };
+
+    using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+    typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+    typename DofTypeIndices::FaceIdx faceIdx;
+
+public:
+    //! The constructor sets the gravity, if desired by the user.
+    RANSProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry)
+    {
+        updateStaticWallProperties();
+    }
+
+    /*!
+     * \brief Update the static (solution independent) relations to the walls
+     *
+     * This function determines all element with a wall intersection,
+     * the wall distances and the relation to the neighboring elements.
+     */
+    void updateStaticWallProperties()
+    {
+        using std::abs;
+        std::cout << "Update static wall properties. ";
+
+        // update size and initial values of the global vectors
+        wallElementIDs_.resize(this->fvGridGeometry().elementMapper().size());
+        wallDistances_.resize(this->fvGridGeometry().elementMapper().size(), std::numeric_limits<Scalar>::max());
+        neighborIDs_.resize(this->fvGridGeometry().elementMapper().size());
+        cellCenters_.resize(this->fvGridGeometry().elementMapper().size(), GlobalPosition(0.0));
+        velocity_.resize(this->fvGridGeometry().elementMapper().size(), DimVector(0.0));
+        velocityMaximum_.resize(this->fvGridGeometry().elementMapper().size(), DimVector(0.0));
+        velocityMinimum_.resize(this->fvGridGeometry().elementMapper().size(), DimVector(std::numeric_limits<Scalar>::max()));
+        velocityGradients_.resize(this->fvGridGeometry().elementMapper().size(), DimMatrix(0.0));
+        flowNormalAxis_.resize(this->fvGridGeometry().elementMapper().size(), 0);
+        wallNormalAxis_.resize(this->fvGridGeometry().elementMapper().size(), 1);
+        kinematicViscosity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+
+        // retrieve all wall intersections and corresponding elements
+        std::vector<unsigned int> wallElements;
+        std::vector<GlobalPosition> wallPositions;
+        std::vector<unsigned int> wallNormalAxisTemp;
+        auto& gridView(this->fvGridGeometry().gridView());
+        for (const auto& element : elements(gridView))
+        {
+            for (const auto& intersection : intersections(gridView, element))
+            {
+                GlobalPosition global = intersection.geometry().center();
+                if (asImp_().isOnWall(global))
+                {
+                    wallElements.push_back(this->fvGridGeometry().elementMapper().index(element));
+                    wallPositions.push_back(global);
+                    for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
+                    {
+                        if (abs(intersection.centerUnitOuterNormal()[dimIdx]) > 1e-10)
+                            wallNormalAxisTemp.push_back(dimIdx);
+                    }
+                }
+            }
+        }
+        std::cout << "NumWallIntersections=" << wallPositions.size() << std::endl;
+
+        // search for shortest distance to wall for each element
+        for (const auto& element : elements(gridView))
+        {
+            unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+            cellCenters_[elementID] = element.geometry().center();
+            for (unsigned int i = 0; i < wallPositions.size(); ++i)
+            {
+                static const int problemWallNormalAxis
+                    = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.WallNormalAxis", -1);
+                int searchAxis = problemWallNormalAxis;
+
+                // search along wall normal axis of the intersection
+                if (problemWallNormalAxis < 0 || problemWallNormalAxis >= dim)
+                {
+                    searchAxis = wallNormalAxisTemp[i];
+                }
+
+                GlobalPosition global = element.geometry().center();
+                global -= wallPositions[i];
+                // second and argument ensures to use only aligned elements
+                if (abs(global[searchAxis]) < wallDistances_[elementID]
+                    && abs(global[searchAxis]) < global.two_norm() + 1e-8
+                    && abs(global[searchAxis]) > global.two_norm() - 1e-8)
+                {
+                    wallDistances_[elementID] = abs(global[searchAxis]);
+                    wallElementIDs_[elementID] = wallElements[i];
+                    wallNormalAxis_[elementID] = searchAxis;
+                }
+            }
+        }
+
+        // search for neighbor IDs
+        for (const auto& element : elements(gridView))
+        {
+            unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+            for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
+            {
+                neighborIDs_[elementID][dimIdx][0] = elementID;
+                neighborIDs_[elementID][dimIdx][1] = elementID;
+            }
+
+            for (const auto& intersection : intersections(gridView, element))
+            {
+                if (intersection.boundary())
+                    continue;
+
+                unsigned int neighborID = this->fvGridGeometry().elementMapper().index(intersection.outside());
+                for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
+                {
+                    if (abs(cellCenters_[elementID][dimIdx] - cellCenters_[neighborID][dimIdx]) > 1e-8)
+                    {
+                        if (cellCenters_[elementID][dimIdx] > cellCenters_[neighborID][dimIdx])
+                        {
+                            neighborIDs_[elementID][dimIdx][0] = neighborID;
+                        }
+                        if (cellCenters_[elementID][dimIdx] < cellCenters_[neighborID][dimIdx])
+                        {
+                            neighborIDs_[elementID][dimIdx][1] = neighborID;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    /*!
+     * \brief Update the dynamic (solution dependent) relations to the walls
+     *
+     * The basic function calcuates the cell-centered velocities and
+     * the respective gradients.
+     * Further, the kinematic viscosity at the wall is stored.
+     *
+     * \param curSol The solution vector.
+     */
+    void updateDynamicWallProperties(const SolutionVector& curSol)
+    {
+        using std::abs;
+        using std::max;
+        using std::min;
+        std::cout << "Update dynamic wall properties." << std::endl;
+
+        static const int flowNormalAxis
+            = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.FlowNormalAxis", -1);
+
+        // calculate cell-center-averaged velocities
+        for (const auto& element : elements(this->fvGridGeometry().gridView()))
+        {
+            auto fvGeometry = localView(this->fvGridGeometry());
+            fvGeometry.bindElement(element);
+            unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+
+            // calculate velocities
+            DimVector velocityTemp(0.0);
+            for (auto&& scvf : scvfs(fvGeometry))
+            {
+                const int dofIdxFace = scvf.dofIndex();
+                const auto numericalSolutionFace = curSol[faceIdx][dofIdxFace][momentumBalanceIdx];
+                velocityTemp[scvf.directionIndex()] += numericalSolutionFace;
+            }
+            for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
+                velocity_[elementID][dimIdx] = velocityTemp[dimIdx] * 0.5; // faces are equidistant to cell center
+        }
+
+        // calculate cell-center-averaged velocity gradients, maximum, and minimum values
+        for (const auto& element : elements(this->fvGridGeometry().gridView()))
+        {
+            unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+            unsigned int wallElementID = wallElementIDs_[elementID];
+
+            Scalar maxVelocity = 0.0;
+            for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
+            {
+                for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
+                {
+                    velocityGradients_[elementID][velIdx][dimIdx]
+                        = (velocity_[neighborIDs_[elementID][dimIdx][1]][velIdx]
+                              - velocity_[neighborIDs_[elementID][dimIdx][0]][velIdx])
+                          / (cellCenters_[neighborIDs_[elementID][dimIdx][1]][dimIdx]
+                              - cellCenters_[neighborIDs_[elementID][dimIdx][0]][dimIdx]);
+                }
+
+                if (abs(velocity_[elementID][dimIdx]) > abs(velocityMaximum_[wallElementID][dimIdx]))
+                {
+                    velocityMaximum_[wallElementID][dimIdx] = velocity_[elementID][dimIdx];
+                }
+                if (abs(velocity_[elementID][dimIdx]) < abs(velocityMinimum_[wallElementID][dimIdx]))
+                {
+                    velocityMinimum_[wallElementID][dimIdx] = velocity_[elementID][dimIdx];
+                }
+
+                if (0 <= flowNormalAxis && flowNormalAxis < dim)
+                {
+                    flowNormalAxis_[elementID] = flowNormalAxis;
+                }
+                else if (abs(maxVelocity) < abs(velocity_[elementID][dimIdx]))
+                {
+                    maxVelocity = abs(velocity_[elementID][dimIdx]);
+                    flowNormalAxis_[elementID] = dimIdx;
+                }
+            }
+        }
+
+        // calculate or call all secondary variables
+        for (const auto& element : elements(this->fvGridGeometry().gridView()))
+        {
+            unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+
+            auto fvGeometry = localView(this->fvGridGeometry());
+            fvGeometry.bindElement(element);
+            for (auto&& scv : scvs(fvGeometry))
+            {
+                const int dofIdx = scv.dofIndex();
+                CellCenterPrimaryVariables priVars(curSol[cellCenterIdx][dofIdx]);
+                auto elemSol = elementSolution<SolutionVector, FVGridGeometry>(std::move(priVars));
+                VolumeVariables volVars;
+                volVars.update(elemSol, asImp_(), element, scv);
+                kinematicViscosity_[elementID] = volVars.viscosity() / volVars.density();
+            }
+        }
+    }
+
+    /*!
+     * \brief Returns whether a given point is on a wall
+     *
+     * \param globalPos The position in global coordinates where the temperature should be specified.
+     */
+    bool isOnWall(const GlobalPosition &globalPos) const
+    {
+        // Throw an exception if no walls are implemented
+        DUNE_THROW(Dune::InvalidStateException,
+                   "The problem does not provide an isOnWall() method.");
+    }
+
+public:
+    std::vector<unsigned int> wallElementIDs_;
+    std::vector<Scalar> wallDistances_;
+    std::vector<std::array<std::array<unsigned int, 2>, dim>> neighborIDs_;
+    std::vector<GlobalPosition> cellCenters_;
+    std::vector<DimVector> velocity_;
+    std::vector<DimVector> velocityMaximum_;
+    std::vector<DimVector> velocityMinimum_;
+    std::vector<DimMatrix> velocityGradients_;
+    std::vector<unsigned int> wallNormalAxis_;
+    std::vector<unsigned int> flowNormalAxis_;
+    std::vector<Scalar> kinematicViscosity_;
+
+private:
+    //! Returns the implementation of the problem (i.e. static polymorphism)
+    Implementation &asImp_()
+    { return *static_cast<Implementation *>(this); }
+
+    //! \copydoc asImp_()
+    const Implementation &asImp_() const
+    { return *static_cast<const Implementation *>(this); }
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/rans/volumevariables.hh b/dumux/freeflow/rans/volumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9ae2e70a8dde7ae4d9fa3869ced247ec89a08738
--- /dev/null
+++ b/dumux/freeflow/rans/volumevariables.hh
@@ -0,0 +1,260 @@
+// -*- 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 RANSModel
+ *
+ * \copydoc Dumux::RANSVolumeVariables
+ */
+#ifndef DUMUX_RANS_VOLUME_VARIABLES_HH
+#define DUMUX_RANS_VOLUME_VARIABLES_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/material/fluidstates/immiscible.hh>
+
+namespace Dumux
+{
+
+// forward declaration
+template <class TypeTag, bool enableEnergyBalance>
+class RANSVolumeVariablesImplementation;
+
+/*!
+ * \ingroup RANSModel
+ * \brief Volume variables for the single-phase Reynolds-Averaged Navier-Stokes models.
+ *        The class is specialized for isothermal and non-isothermal models.
+ */
+template <class TypeTag>
+using RANSVolumeVariables = RANSVolumeVariablesImplementation<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>;
+
+/*!
+ * \ingroup RANSModel
+ * \brief Volume variables for the isothermal single-phase Reynolds-Averaged Navier-Stokes models.
+ */
+template <class TypeTag>
+class RANSVolumeVariablesImplementation<TypeTag, false>
+: public NavierStokesVolumeVariablesImplementation<TypeTag, false>
+{
+    using ParentType = NavierStokesVolumeVariablesImplementation<TypeTag, false>;
+    using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+
+    enum { dimWorld = GridView::dimensionworld };
+    using DimVector = Dune::FieldVector<Scalar, dimWorld>;
+    using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+
+    static const int defaultPhaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+
+public:
+
+    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+    /*!
+     * \brief Update all quantities for a given control volume
+     *
+     * Wall related quantities are stored and the calculateEddyViscosity(...)
+     * function of the turbulence model implementation is called.
+     *
+     * \param elemSol A vector containing all primary variables connected to the element
+     * \param problem The object specifying the problem which ought to
+     *                be simulated
+     * \param element An element which contains part of the control volume
+     * \param scv The sub-control volume
+     */
+    template<class ElementSolution>
+    void update(const ElementSolution &elemSol,
+                const Problem &problem,
+                const Element &element,
+                const SubControlVolume& scv)
+    {
+        using std::abs;
+        using std::max;
+        using std::sqrt;
+        ParentType::update(elemSol, problem, element, scv);
+
+        // calculate characteristic properties of the turbulent flow
+        elementID_ = problem.fvGridGeometry().elementMapper().index(element);
+        wallElementID_ = problem.wallElementIDs_[elementID_];
+        wallDistance_ = problem.wallDistances_[elementID_];
+        velocity_ = problem.velocity_[elementID_];
+        velocityGradients_ = problem.velocityGradients_[elementID_];
+        unsigned int flowNormalAxis = problem.flowNormalAxis_[elementID_];
+        unsigned int wallNormalAxis = problem.wallNormalAxis_[elementID_];
+        Scalar uStar = sqrt(problem.kinematicViscosity_[wallElementID_]
+                            * abs(problem.velocityGradients_[wallElementID_][flowNormalAxis][wallNormalAxis]));
+        yPlus_ = wallDistance_ * uStar / asImp_().kinematicViscosity();
+        yPlus_ = max(yPlus_, 1e-10); // zero values lead to numerical problems in some turbulence models
+        uPlus_ = velocity_[flowNormalAxis] / max(uStar, 1e-10);
+    };
+
+    /*!
+     * \brief Return the velocity vector \f$\mathrm{[m/s]}\f$ at the control volume center.
+     */
+    DimVector velocity() const
+    { return velocity_; }
+
+    /*!
+     * \brief Return the velocity gradients \f$\mathrm{[1/s]}\f$ at the control volume center.
+     */
+    DimMatrix velocityGradients() const
+    { return velocityGradients_; }
+
+    /*!
+     * \brief Return the wall distance \f$\mathrm{[m]}\f$ of the control volume.
+     */
+    Scalar wallDistance() const
+    { return wallDistance_; }
+
+    /*!
+     * \brief Return the dimensionless wall distance \f$\mathrm{[-]}\f$.
+     */
+    Scalar yPlus() const
+    { return yPlus_; }
+
+    /*!
+     * \brief Return the dimensionless velocity \f$\mathrm{[-]}\f$.
+     */
+    Scalar uPlus() const
+    { return uPlus_; }
+
+    /*!
+     * \brief Set the values of the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$ within the
+     *        control volume.
+     */
+    void setDynamicEddyViscosity(Scalar value)
+    { dynamicEddyViscosity_ = value; }
+
+    /*!
+     * \brief Return the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$ of the flow within the
+     *        control volume.
+     */
+    Scalar dynamicEddyViscosity() const
+    { return dynamicEddyViscosity_; }
+
+    /*!
+     * \brief Return the effective dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
+     *        control volume.
+     */
+    Scalar effectiveViscosity() const
+    { return asImp_().viscosity() + dynamicEddyViscosity(); }
+
+    /*!
+     * \brief Return the kinematic eddy viscosity \f$\mathrm{[m^2/s]}\f$ of the fluid within the
+     *        control volume.
+     */
+    Scalar kinematicViscosity() const
+    { return asImp_().viscosity() / asImp_().density(); }
+
+private:
+    //! Returns the implementation of the problem (i.e. static polymorphism)
+    Implementation &asImp_()
+    { return *static_cast<Implementation *>(this); }
+
+    //! \copydoc asImp_()
+    const Implementation &asImp_() const
+    { return *static_cast<const Implementation *>(this); }
+
+protected:
+    DimVector velocity_;
+    DimMatrix velocityGradients_;
+    Scalar dynamicEddyViscosity_;
+    unsigned int elementID_;
+    unsigned int wallElementID_;
+    Scalar wallDistance_;
+    Scalar yPlus_;
+    Scalar uPlus_;
+};
+
+/*!
+ * \ingroup RANSModel
+ * \brief Volume variables for the non-isothermal single-phase Reynolds-Averaged Navier-Stokes models.
+ */
+template <class TypeTag>
+class RANSVolumeVariablesImplementation<TypeTag, true>
+: public NavierStokesVolumeVariablesImplementation<TypeTag, true>,
+  public RANSVolumeVariablesImplementation<TypeTag, false>
+{
+    using ParentTypeNonIsothermal = NavierStokesVolumeVariablesImplementation<TypeTag, true>;
+    using ParentTypeIsothermal = RANSVolumeVariablesImplementation<TypeTag, false>;
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+    static const int defaultPhaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+//     static const int temperatureIdx = Indices::temperatureIdx;
+
+public:
+
+    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+    /*!
+     * \brief Update all quantities for a given control volume
+     *
+     * \param elemSol A vector containing all primary variables connected to the element
+     * \param problem The object specifying the problem which ought to
+     *                be simulated
+     * \param element An element which contains part of the control volume
+     * \param scv The sub-control volume
+     */
+    template<class ElementSolution>
+    void update(const ElementSolution &elemSol,
+                const Problem &problem,
+                const Element &element,
+                const SubControlVolume &scv)
+    {
+        ParentTypeIsothermal::update(elemSol, problem, element, scv);
+        ParentTypeNonIsothermal::update(elemSol, problem, element, scv);
+        // TODO: convert eddyViscosity to eddyThermalConductivity
+    }
+
+    /*!
+     * \brief Returns the eddy thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
+     *        of the flow phase in the sub-control volume.
+     */
+    Scalar eddyThermalConductivity() const
+    { return eddyThermalConductivity_; }
+
+    /*!
+     * \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
+     *        of the fluid-flow in the sub-control volume.
+     */
+    Scalar effectiveThermalConductivity() const
+    {
+        return FluidSystem::thermalConductivity(this->fluidState_)
+               + eddyThermalConductivity();
+    }
+
+protected:
+    Scalar eddyThermalConductivity_;
+};
+}
+
+#endif
diff --git a/dumux/freeflow/rans/vtkoutputfields.hh b/dumux/freeflow/rans/vtkoutputfields.hh
new file mode 100644
index 0000000000000000000000000000000000000000..631bb23006bb565e834c90c69ef3e0e34a25e26b
--- /dev/null
+++ b/dumux/freeflow/rans/vtkoutputfields.hh
@@ -0,0 +1,98 @@
+// -*- 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 RANSModel
+ * \copydoc Dumux::RANSVtkOutputFields
+ */
+#ifndef DUMUX_RANS_VTK_OUTPUT_FIELDS_HH
+#define DUMUX_RANS_VTK_OUTPUT_FIELDS_HH
+
+#include <dune/common/fvector.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup RANSModel
+ * \brief Adds vtk output fields for the Reynolds-Averaged Navier-Stokes model
+ */
+template<class FVGridGeometry>
+class RANSVtkOutputFields
+{
+    using ctype = typename FVGridGeometry::GridView::ctype;
+    using GlobalPosition = Dune::FieldVector<ctype, FVGridGeometry::GridView::dimensionworld>;
+
+    // Helper type used for tag dispatching (to add discretization-specific fields).
+    template<DiscretizationMethod discMethod>
+    using discMethodTag = std::integral_constant<DiscretizationMethod, discMethod>;
+
+public:
+    //! Initialize the Navier-Stokes specific vtk output fields.
+    template <class VtkOutputModule>
+    static void init(VtkOutputModule& vtk)
+    {
+        vtk.addVolumeVariable([](const auto& v){ return v.velocity()[0]; }, "v_x [m/s]");
+        vtk.addVolumeVariable([](const auto& v){ return v.velocityGradients()[0][0]; }, "dv_x/dx [m/s]");
+        vtk.addVolumeVariable([](const auto& v){ return v.velocityGradients()[0][1]; }, "dv_x/dy [m/s]");
+        vtk.addVolumeVariable([](const auto& v){ return v.pressure(); }, "p [Pa]");
+        vtk.addVolumeVariable([](const auto& v){ return v.pressure() - 1e5; }, "p_rel [Pa]");
+        vtk.addVolumeVariable([](const auto& v){ return v.density(); }, "rho [kg/m^3]");
+        vtk.addVolumeVariable([](const auto& v){ return v.viscosity() / v.density(); }, "nu [m^2/s]");
+        vtk.addVolumeVariable([](const auto& v){ return v.dynamicEddyViscosity() / v.density(); }, "nu_t [m^2/s]");
+        vtk.addVolumeVariable([](const auto& v){ return v.wallDistance(); }, "l_w [m]");
+        vtk.addVolumeVariable([](const auto& v){ return v.yPlus(); }, "y^+ [-]");
+        vtk.addVolumeVariable([](const auto& v){ return v.uPlus(); }, "u^+ [-]");
+
+        // add discretization-specific fields
+        additionalOutput_(vtk, discMethodTag<FVGridGeometry::discMethod>{});
+    }
+
+private:
+
+    //! Adds discretization-specific fields (nothing by default).
+    template <class VtkOutputModule, class AnyMethod>
+    static void additionalOutput_(VtkOutputModule& vtk, AnyMethod)
+    { }
+
+    //! Adds discretization-specific fields (velocity vectors on the faces for the staggered discretization).
+    template <class VtkOutputModule>
+    static void additionalOutput_(VtkOutputModule& vtk, discMethodTag<DiscretizationMethod::staggered>)
+    {
+        const bool writeFaceVars = getParamFromGroup<bool>(vtk.paramGroup(), "Vtk.WriteFaceData", false);
+        if(writeFaceVars)
+        {
+            auto faceVelocityVector = [](const typename FVGridGeometry::SubControlVolumeFace& scvf, const auto& faceVars)
+                                      {
+                                          GlobalPosition velocity(0.0);
+                                          velocity[scvf.directionIndex()] = faceVars.velocitySelf();
+                                          return velocity;
+                                      };
+
+            vtk.addFaceVariable(faceVelocityVector, "faceVelocity");
+        }
+    }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/freeflow/rans/zeroeq/CMakeLists.txt b/dumux/freeflow/rans/zeroeq/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..811156a459b9ad0ac3fa9a641e2162c6cfa3aa75
--- /dev/null
+++ b/dumux/freeflow/rans/zeroeq/CMakeLists.txt
@@ -0,0 +1,6 @@
+#install headers
+install(FILES
+indices.hh
+model.hh
+volumevariables.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/rans/zeroeq)
diff --git a/dumux/freeflow/rans/zeroeq/indices.hh b/dumux/freeflow/rans/zeroeq/indices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4cd6b208df6a8989d966816c6452176118928d25
--- /dev/null
+++ b/dumux/freeflow/rans/zeroeq/indices.hh
@@ -0,0 +1,53 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup ZeroEqModel
+ * \copydoc Dumux::ZeroEqIndices
+ */
+#ifndef DUMUX_ZEROEQ_INDICES_HH
+#define DUMUX_ZEROEQ_INDICES_HH
+
+#include <dumux/freeflow/navierstokes/indices.hh>
+
+namespace Dumux
+{
+// \{
+/*!
+ * \ingroup ZeroEqModel
+ * \brief The common indices for the isothermal ZeroEq model.
+ *
+ * \tparam dimension The dimension of the problem
+ * \tparam numEquations the number of model equations
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <int dimension, int numEquations, int PVOffset = 0>
+struct ZeroEqIndices
+    : NavierStokesIndices<dimension, numEquations, PVOffset>
+{
+    static constexpr int noEddyViscosityModel = 0;
+    static constexpr int prandtl = 1;
+    static constexpr int modifiedVanDriest = 2;
+    static constexpr int baldwinLomax = 3;
+};
+
+// \}
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/rans/zeroeq/model.hh b/dumux/freeflow/rans/zeroeq/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e76bf1144222b04911c800125a7eb11b44a2d0cb
--- /dev/null
+++ b/dumux/freeflow/rans/zeroeq/model.hh
@@ -0,0 +1,86 @@
+// -*- 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 ZeroEqModel
+ *
+ * \brief A single-phase, isothermal Reynolds-Averaged Navier-Stokes 0-Eq. model
+ *
+ * \copydoc RANSModel
+ *
+ * These models calculate the eddy viscosity without solving additional PDEs,
+ * only based on the wall distance and the velocity gradient.
+ * The following models are available:
+ *  -# Prandtl's mixing length, e.g. \cite Oertel2012a
+ *  -# Van-Driest modification, \cite vanDriest1956a and \cite Hanna1981a
+ *  -# Baldwin-Lomax, \cite Baldwin1978a
+ */
+
+#ifndef DUMUX_ZEROEQ_MODEL_HH
+#define DUMUX_ZEROEQ_MODEL_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/freeflow/properties.hh>
+#include <dumux/freeflow/rans/model.hh>
+
+#include "indices.hh"
+#include "volumevariables.hh"
+
+namespace Dumux
+{
+
+// \{
+///////////////////////////////////////////////////////////////////////////
+// properties for the single-phase Reynolds-Averaged Navier-Stokes 0-Eq. model
+///////////////////////////////////////////////////////////////////////////
+namespace Properties {
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the single-phase, isothermal Reynolds-Averaged Navier-Stokes 0-Eq. model
+NEW_TYPE_TAG(ZeroEq, INHERITS_FROM(RANS));
+
+//! The type tag for the corresponding non-isothermal model
+NEW_TYPE_TAG(ZeroEqNI, INHERITS_FROM(ZeroEq, RANSNI));
+
+///////////////////////////////////////////////////////////////////////////
+// default property values for the isothermal single phase model
+///////////////////////////////////////////////////////////////////////////
+
+//! The indices
+SET_PROP(ZeroEq, Indices)
+{
+private:
+    static constexpr int numEq = GET_PROP_VALUE(TypeTag, NumEq);
+    static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
+public:
+    using type = ZeroEqIndices<dim, numEq>;
+};
+
+//! The volume variables
+SET_TYPE_PROP(ZeroEq, VolumeVariables, ZeroEqVolumeVariables<TypeTag>);
+
+// \}
+}
+
+} // end namespace
+
+#endif // DUMUX_ZEROEQ_MODEL_HH
diff --git a/dumux/freeflow/rans/zeroeq/problem.hh b/dumux/freeflow/rans/zeroeq/problem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2c0841d90c44ea5a8da050b40a3bc60e4067aa8f
--- /dev/null
+++ b/dumux/freeflow/rans/zeroeq/problem.hh
@@ -0,0 +1,234 @@
+// -*- 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 ZeroEqModel
+ * \copydoc Dumux::ZeroEqProblem
+ */
+#ifndef DUMUX_ZEROEQ_PROBLEM_HH
+#define DUMUX_ZEROEQ_PROBLEM_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/staggeredfvproblem.hh>
+#include <dumux/discretization/localview.hh>
+#include <dumux/discretization/staggered/elementsolution.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/freeflow/rans/problem.hh>
+
+#include "model.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ZeroEqModel
+ * \brief Zero-equation turbulence problem base class.
+ *
+ * This implements some base functionality for zero-equation models
+ * and a routine for the determining the eddy viscosity of the Baldwin-Lomax model.
+ */
+template<class TypeTag>
+class ZeroEqProblem : public RANSProblem<TypeTag>
+{
+    using ParentType = RANSProblem<TypeTag>;
+    using Implementation = typename GET_PROP_TYPE(TypeTag, Problem);
+
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Grid = typename GridView::Grid;
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+    enum {
+        dim = Grid::dimension,
+      };
+    using GlobalPosition = Dune::FieldVector<Scalar, dim>;
+    using DimVector = Dune::FieldVector<Scalar, dim>;
+    using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
+
+    enum {
+        massBalanceIdx = Indices::massBalanceIdx,
+        momentumBalanceIdx = Indices::momentumBalanceIdx
+    };
+
+    using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+    typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+    typename DofTypeIndices::FaceIdx faceIdx;
+
+public:
+    //! The constructor sets the gravity, if desired by the user.
+    ZeroEqProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry)
+    {
+        updateStaticWallProperties();
+    }
+
+    /*!
+     * \brief Correct size of the static (solution independent) wall variables
+     */
+    void updateStaticWallProperties()
+    {
+        // update size and initial values of the global vectors
+        kinematicEddyViscosity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+    }
+
+    /*!
+     * \brief Update the dynamic (solution dependent) relations to the walls
+     *
+     * This calculate all necessary information for the Baldwin-Lomax turbulence model
+     *
+     * \param curSol The solution vector.
+     */
+    void updateDynamicWallProperties(const SolutionVector& curSol)
+    {
+        ParentType::updateDynamicWallProperties(curSol);
+
+        std::vector<Scalar> kinematicEddyViscosityInner(this->fvGridGeometry().elementMapper().size(), 0.0);
+        std::vector<Scalar> kinematicEddyViscosityOuter(this->fvGridGeometry().elementMapper().size(), 0.0);
+        std::vector<Scalar> kinematicEddyViscosityDifference(this->fvGridGeometry().elementMapper().size(), 0.0);
+        std::vector<Scalar> switchingPosition(this->fvGridGeometry().elementMapper().size(), std::numeric_limits<Scalar>::max());
+
+        static const int eddyViscosityModel
+            = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.EddyViscosityModel");
+        if (eddyViscosityModel != Indices::baldwinLomax)
+            return;
+
+        using std::abs;
+        using std::exp;
+        using std::min;
+        using std::pow;
+        using std::sqrt;
+        const Scalar aPlus = 26.0;
+        const Scalar k = 0.0168;
+        const Scalar cCP = 1.6;
+        const Scalar cWake = 0.25;
+        const Scalar cKleb = 0.3;
+
+        static const Scalar karmanConstant
+            = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.KarmanConstant");
+
+        std::vector<Scalar> storedFMax;
+        std::vector<Scalar> storedYFMax;
+        storedFMax.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+        storedYFMax.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+
+        // (1) calculate inner viscosity and Klebanoff function
+        for (const auto& element : elements(this->fvGridGeometry().gridView()))
+        {
+            unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+            unsigned int wallElementID = this->wallElementIDs_[elementID];
+            Scalar wallDistance = this->wallDistances_[elementID];
+            unsigned int flowNormalAxis = this->flowNormalAxis_[elementID];
+            unsigned int wallNormalAxis = this->wallNormalAxis_[elementID];
+
+            Scalar omegaAbs = abs(this->velocityGradients_[elementID][flowNormalAxis][wallNormalAxis]
+                                  - this->velocityGradients_[elementID][wallNormalAxis][flowNormalAxis]);
+            Scalar uStar = sqrt(this->kinematicViscosity_[wallElementID]
+                                * abs(this->velocityGradients_[wallElementID][flowNormalAxis][wallNormalAxis]));
+            Scalar yPlus = wallDistance * uStar / this->kinematicViscosity_[elementID];
+            Scalar mixingLength = karmanConstant * wallDistance * (1.0 - exp(-yPlus / aPlus));
+            kinematicEddyViscosityInner[elementID] = mixingLength * mixingLength * omegaAbs;
+
+            Scalar f = wallDistance * omegaAbs * (1.0 - exp(-yPlus / aPlus));
+            if (f > storedFMax[wallElementID])
+            {
+                storedFMax[wallElementID] = f;
+                storedYFMax[wallElementID] = wallDistance;
+            }
+        }
+
+        // (2) calculate outer viscosity
+        for (const auto& element : elements(this->fvGridGeometry().gridView()))
+        {
+            unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+            unsigned int wallElementID = this->wallElementIDs_[elementID];
+            Scalar wallDistance = this->wallDistances_[elementID];
+
+            Scalar maxVelocityNorm = 0.0;
+            Scalar minVelocityNorm = 0.0;
+            for (unsigned dimIdx = 0; dimIdx < dim; ++dimIdx)
+            {
+                maxVelocityNorm += this->velocityMaximum_[wallElementID][dimIdx]
+                                   * this->velocityMaximum_[wallElementID][dimIdx];
+                minVelocityNorm += this->velocityMinimum_[wallElementID][dimIdx]
+                                   * this->velocityMinimum_[wallElementID][dimIdx];
+            }
+            Scalar deltaU = sqrt(maxVelocityNorm) - sqrt(minVelocityNorm);
+            Scalar yFMax = storedYFMax[wallElementID];
+            Scalar fMax = storedFMax[wallElementID];
+            Scalar fWake = min(yFMax * fMax, cWake * yFMax * deltaU * deltaU / fMax);
+            Scalar fKleb = 1.0 / (1.0 + 5.5 * pow(cKleb * wallDistance / yFMax, 6.0));
+            kinematicEddyViscosityOuter[elementID] = k * cCP * fWake * fKleb;
+
+            kinematicEddyViscosityDifference[elementID]
+              = kinematicEddyViscosityInner[elementID] - kinematicEddyViscosityOuter[elementID];
+        }
+
+        // (3) switching point
+        for (const auto& element : elements(this->fvGridGeometry().gridView()))
+        {
+            unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+            unsigned int wallElementID = this->wallElementIDs_[elementID];
+            Scalar wallDistance = this->wallDistances_[elementID];
+
+            // checks if sign switches, by multiplication
+            Scalar check = kinematicEddyViscosityDifference[wallElementID] * kinematicEddyViscosityDifference[elementID];
+            if (check < 0 // means sign has switched
+                && switchingPosition[wallElementID] > wallDistance)
+            {
+                switchingPosition[wallElementID] = wallDistance;
+            }
+        }
+
+        // (4) finally determine eddy viscosity
+        for (const auto& element : elements(this->fvGridGeometry().gridView()))
+        {
+            unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+            unsigned int wallElementID = this->wallElementIDs_[elementID];
+            Scalar wallDistance = this->wallDistances_[elementID];
+
+            kinematicEddyViscosity_[elementID] = (wallDistance >= switchingPosition[wallElementID])
+                                                 ? kinematicEddyViscosityOuter[elementID]
+                                                 : kinematicEddyViscosityInner[elementID];
+        }
+    }
+
+public:
+    std::vector<Scalar> kinematicEddyViscosity_;
+
+private:
+    //! Returns the implementation of the problem (i.e. static polymorphism)
+    Implementation &asImp_()
+    { return *static_cast<Implementation *>(this); }
+
+    //! \copydoc asImp_()
+    const Implementation &asImp_() const
+    { return *static_cast<const Implementation *>(this); }
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/rans/zeroeq/volumevariables.hh b/dumux/freeflow/rans/zeroeq/volumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..417bde0e9708bffaf843fb203993ef6ba863ad3d
--- /dev/null
+++ b/dumux/freeflow/rans/zeroeq/volumevariables.hh
@@ -0,0 +1,164 @@
+// -*- 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 ZeroEqModel
+ *
+ * \copydoc Dumux::ZeroEqVolumeVariables
+ */
+#ifndef DUMUX_ZEROEQ_VOLUME_VARIABLES_HH
+#define DUMUX_ZEROEQ_VOLUME_VARIABLES_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/material/fluidstates/immiscible.hh>
+
+namespace Dumux
+{
+
+// forward declaration
+template <class TypeTag, bool enableEnergyBalance>
+class ZeroEqVolumeVariablesImplementation;
+
+/*!
+ * \ingroup ZeroEqModel
+ * \brief Volume variables for the single-phase 0-Eq. model.
+ *        The class is specialized for isothermal and non-isothermal models.
+ */
+template <class TypeTag>
+using ZeroEqVolumeVariables = ZeroEqVolumeVariablesImplementation<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>;
+
+/*!
+ * \ingroup ZeroEqModel
+ * \brief Volume variables for the isothermal single-phase 0-Eq. model.
+ */
+template <class TypeTag>
+class ZeroEqVolumeVariablesImplementation<TypeTag, false>
+: public RANSVolumeVariablesImplementation<TypeTag, false>
+{
+    using ParentType = RANSVolumeVariablesImplementation<TypeTag, false>;
+    using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+
+    static const int defaultPhaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+
+public:
+    /*!
+     * \brief Update all quantities for a given control volume
+     *
+     * \param elemSol A vector containing all primary variables connected to the element
+     * \param problem The object specifying the problem which ought to
+     *                be simulated
+     * \param element An element which contains part of the control volume
+     * \param scv The sub-control volume
+     */
+    template<class ElementSolution>
+    void update(const ElementSolution &elemSol,
+                const Problem &problem,
+                const Element &element,
+                const SubControlVolume& scv)
+    {
+        ParentType::update(elemSol, problem, element, scv);
+        calculateEddyViscosity(elemSol, problem, element, scv);
+    }
+
+
+    /*!
+     * \brief Calculate and set the dynamic eddy viscosity.
+     *
+     * \param elemSol A vector containing all primary variables connected to the element
+     * \param problem The object specifying the problem which ought to
+     *                be simulated
+     * \param element An element which contains part of the control volume
+     * \param scv The sub-control volume
+     */
+    template<class ElementSolution>
+    void calculateEddyViscosity(const ElementSolution &elemSol,
+                                const Problem &problem,
+                                const Element &element,
+                                const SubControlVolume& scv)
+    {
+        using std::abs;
+        using std::exp;
+        using std::sqrt;
+        Scalar kinematicEddyViscosity = 0.0;
+        static const Scalar karmanConstant
+            = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.KarmanConstant");
+        static const int eddyViscosityModel
+            = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.EddyViscosityModel");
+        unsigned int elementID = problem.fvGridGeometry().elementMapper().index(element);
+        unsigned int flowNormalAxis = problem.flowNormalAxis_[elementID];
+        unsigned int wallNormalAxis = problem.wallNormalAxis_[elementID];
+        Scalar velGrad = abs(asImp_().velocityGradients()[flowNormalAxis][wallNormalAxis]);
+
+        if (eddyViscosityModel == Indices::noEddyViscosityModel)
+        {
+            // kinematicEddyViscosity = 0.0
+        }
+        else if (eddyViscosityModel == Indices::prandtl)
+        {
+            Scalar mixingLength = karmanConstant * asImp_().wallDistance();
+            kinematicEddyViscosity = mixingLength * mixingLength * velGrad;
+        }
+        else if (eddyViscosityModel == Indices::modifiedVanDriest)
+        {
+            Scalar mixingLength = karmanConstant * asImp_().wallDistance()
+                                  * (1.0 - exp(-asImp_().yPlus() / 26.0))
+                                  / sqrt(1.0 - exp(-0.26 * asImp_().yPlus()));
+            kinematicEddyViscosity = mixingLength * mixingLength * velGrad;
+        }
+        else if (eddyViscosityModel == Indices::baldwinLomax)
+        {
+            kinematicEddyViscosity = problem.kinematicEddyViscosity_[elementID];
+        }
+        else
+        {
+            DUNE_THROW(Dune::NotImplemented,
+                       "This eddy viscosity model is not implemented: " << eddyViscosityModel);
+        }
+        asImp_().setDynamicEddyViscosity(kinematicEddyViscosity * asImp_().density());
+    }
+
+private:
+    //! Returns the implementation of the problem (i.e. static polymorphism)
+    Implementation &asImp_()
+    { return *static_cast<Implementation *>(this); }
+
+    //! \copydoc asImp_()
+    const Implementation &asImp_() const
+    { return *static_cast<const Implementation *>(this); }
+};
+
+/*!
+ * \ingroup ZeroEqModel
+ * \brief Volume variables for the non-isothermal single-phase 0-Eq. model.
+ */
+template <class TypeTag>
+class ZeroEqVolumeVariablesImplementation<TypeTag, true>
+: public RANSVolumeVariablesImplementation<TypeTag, true>
+{ };
+}
+
+#endif
diff --git a/test/freeflow/CMakeLists.txt b/test/freeflow/CMakeLists.txt
index a2cd437fb9d68eb2705f3fbc35ae46184a45b7cf..5fa27535b21e7a14de55bba0527bca6371711eda 100644
--- a/test/freeflow/CMakeLists.txt
+++ b/test/freeflow/CMakeLists.txt
@@ -1,2 +1,3 @@
 add_subdirectory("navierstokes")
 add_subdirectory("navierstokesnc")
+add_subdirectory("rans")
diff --git a/test/freeflow/rans/CMakeLists.txt b/test/freeflow/rans/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..0b3b205ce7aa702287aaeadd9182b8d4d3d3a1df
--- /dev/null
+++ b/test/freeflow/rans/CMakeLists.txt
@@ -0,0 +1,11 @@
+add_input_file_links()
+dune_symlink_to_source_files(FILES laufer_re50000_u+y+.csv)
+
+dune_add_test(NAME test_pipe_laufer
+              SOURCES test_pipe_laufer.cc
+              CMAKE_GUARD HAVE_UMFPACK
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS  --script fuzzy
+                        --files ${CMAKE_SOURCE_DIR}/test/references/pipe_laufer_zeroeq.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/pipe_laufer_zeroeq_reference-00021.vtu
+                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_pipe_laufer test_pipe_laufer_reference.input")
diff --git a/test/freeflow/rans/laufer_re50000_u+y+.csv b/test/freeflow/rans/laufer_re50000_u+y+.csv
new file mode 100644
index 0000000000000000000000000000000000000000..547f79ab0ced283e16dfe40aa6867440cd4b9446
--- /dev/null
+++ b/test/freeflow/rans/laufer_re50000_u+y+.csv
@@ -0,0 +1,46 @@
+###########
+# Dimensionless velocity profile for Re=50,000
+###########
+# Original source:
+# Laufer, J.
+# The structure of turbulence in fully developed pipe flow
+# NACA Report, 1954, 1174, 417-434
+###########
+# Data taken from:
+# Truckenbrodt, E.
+# Grundlagen und elementare Strömungsvorgänge dichtebeständiger Fluide
+# Fluidmechanik, Springer, 2008, XXVII, 364 , doi: 10.1007/978-3-540-79018-1
+###########
+#Y+,[-],U+,[-]
+3.12,3.243
+3.572,3.649
+4.289,4.135
+4.911,5.068
+5.327,5.27
+5.436,4.662
+5.817,5.473
+6.439,6.162
+7.373,7.095
+7.731,7.5
+8.674,7.905
+8.793,8.311
+9.345,7.5
+10.136,8.716
+10.847,9.203
+13.29,10.014
+16.284,11.351
+17.425,11.635
+26.161,13.054
+43.475,14.676
+66.608,15.486
+100,16.703
+166.184,17.838
+233.156,18.932
+295.521,19.541
+400.812,20.554
+508.022,21.568
+609.95,22.378
+712.756,22.986
+816.14,23.311
+915.724,23.716
+1034.441,23.797
diff --git a/test/freeflow/rans/pipelauferproblem.hh b/test/freeflow/rans/pipelauferproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a9b0e0c68bb52a393fdff2cb3bceafdda429b45a
--- /dev/null
+++ b/test/freeflow/rans/pipelauferproblem.hh
@@ -0,0 +1,235 @@
+// -*- 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 RANSTests
+ * \brief Pipe flow test for the staggered grid RANS model
+ *
+ * This test simulates is based on pipe flow experiments by
+ * John Laufers experiments in 1954 \cite Laufer1954a.
+ */
+#ifndef DUMUX_PIPE_LAUFER_PROBLEM_HH
+#define DUMUX_PIPE_LAUFER_PROBLEM_HH
+
+#include <dumux/material/fluidsystems/gasphase.hh>
+#include <dumux/material/components/air.hh>
+
+#include <dumux/freeflow/rans/zeroeq/model.hh>
+#include <dumux/freeflow/rans/zeroeq/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class PipeLauferProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(PipeLauferProblem, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEq));
+
+// the fluid system
+SET_PROP(PipeLauferProblem, FluidSystem)
+{
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using type = FluidSystems::GasPhase<Scalar, Air<Scalar> >;
+};
+
+// Set the grid type
+SET_TYPE_PROP(PipeLauferProblem, Grid,
+              Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+// Set the problem property
+SET_TYPE_PROP(PipeLauferProblem, Problem, Dumux::PipeLauferProblem<TypeTag> );
+
+SET_BOOL_PROP(PipeLauferProblem, EnableFVGridGeometryCache, true);
+
+SET_BOOL_PROP(PipeLauferProblem, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(PipeLauferProblem, EnableGridVolumeVariablesCache, true);
+}
+
+/*!
+ * \ingroup NavierStokesTests
+ * \brief  Test problem for the one-phase (Navier-) Stokes problem in a channel.
+ *
+ * This test simulates is based on pipe flow experiments by
+ * John Laufers experiments in 1954 \cite Laufer1954a.
+ */
+template <class TypeTag>
+class PipeLauferProblem : public ZeroEqProblem<TypeTag>
+{
+    using ParentType = ZeroEqProblem<TypeTag>;
+
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+    // copy some indices for convenience
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+    enum { dimWorld = GridView::dimensionworld };
+    enum {
+        massBalanceIdx = Indices::massBalanceIdx,
+        momentumBalanceIdx = Indices::momentumBalanceIdx,
+        pressureIdx = Indices::pressureIdx,
+        velocityXIdx = Indices::velocityXIdx,
+        velocityYIdx = Indices::velocityYIdx
+    };
+
+    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+
+    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+
+    using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
+
+public:
+    PipeLauferProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry), eps_(1e-6)
+    {
+        inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
+    }
+
+   /*!
+     * \name Problem parameters
+     */
+    // \{
+
+    bool isOnWall(const GlobalPosition &globalPos) const
+    {
+        Scalar localEps_ = 1e-6; // cannot use the epsilon, because ParentType is initialized first
+        return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + localEps_
+              || globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - localEps_;
+    }
+
+    bool shouldWriteRestartFile() const
+    {
+        return false;
+    }
+
+   /*!
+     * \brief Return the temperature within the domain in [K].
+     *
+     * This problem assumes a temperature of 10 degrees Celsius.
+     */
+    Scalar temperature() const
+    { return 273.15 + 10; } // 10C
+
+   /*!
+     * \brief Return the sources within the domain.
+     *
+     * \param globalPos The global position
+     */
+    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
+    {
+        return SourceValues(0.0);
+    }
+    // \}
+   /*!
+     * \name Boundary conditions
+     */
+    // \{
+
+   /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary control volume.
+     *
+     * \param globalPos The position of the center of the finite volume
+     */
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+    {
+        BoundaryTypes values;
+
+        // set Dirichlet values for the velocity everywhere
+        values.setDirichlet(momentumBalanceIdx);
+
+        // set a fixed pressure in one cell
+        if (isOutlet(globalPos))
+        {
+            values.setDirichlet(massBalanceIdx);
+            values.setOutflow(momentumBalanceIdx);
+        }
+        else
+            values.setOutflow(massBalanceIdx);
+
+        return values;
+    }
+
+   /*!
+     * \brief Evaluate the boundary conditions for a dirichlet
+     *        control volume.
+     *
+     * \param globalPos The center of the finite volume which ought to be set.
+     */
+    PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+    {
+        return initialAtPos(globalPos);
+    }
+
+    // \}
+
+   /*!
+     * \name Volume terms
+     */
+    // \{
+
+   /*!
+     * \brief Evaluate the initial value for a control volume.
+     *
+     * \param globalPos The global position
+     */
+    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+    {
+        PrimaryVariables values(0.0);
+        values[pressureIdx] = 1.0e+5;
+        if(!isOnWall(globalPos))
+        {
+            values[velocityXIdx] = inletVelocity_;
+        }
+
+        return values;
+    }
+
+    // \}
+    void setTimeLoop(TimeLoopPtr timeLoop)
+    {
+        timeLoop_ = timeLoop;
+        if(inletVelocity_ > eps_)
+            timeLoop_->setCheckPoint({10.0, 50.0, 100.});
+    }
+
+    Scalar time() const
+    {
+        return timeLoop_->time();
+    }
+
+private:
+    bool isOutlet(const GlobalPosition& globalPos) const
+    {
+        return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_;
+    }
+
+    Scalar eps_;
+    Scalar inletVelocity_;
+    TimeLoopPtr timeLoop_;
+};
+} //end namespace
+
+#endif
diff --git a/test/freeflow/rans/test_pipe_laufer.cc b/test/freeflow/rans/test_pipe_laufer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..afd6e99c4636c8bea64c965ef793d8aa1bcb8b65
--- /dev/null
+++ b/test/freeflow/rans/test_pipe_laufer.cc
@@ -0,0 +1,278 @@
+// -*- 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 Pipe flow test for the staggered grid RANS model
+ *
+ * This test simulates is based on pipe flow experiments by
+ * John Laufers experiments in 1954 \cite Laufer1954a.
+ */
+ #include <config.h>
+
+ #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 "pipelauferproblem.hh"
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include <dumux/assembly/staggeredfvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/gnuplotinterface.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+
+#include <dumux/freeflow/navierstokes/staggered/fluxoverplane.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ *        reading in parameters.
+ *
+ * \param progName  The name of the program, that was tried to be started.
+ * \param errorMsg  The error message that was issued by the start function.
+ *                  Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+    if (errorMsg.size() > 0) {
+        std::string errorMessageOut = "\nUsage: ";
+                    errorMessageOut += progName;
+                    errorMessageOut += " [options]\n";
+                    errorMessageOut += errorMsg;
+                    errorMessageOut += "\nPlease use the provided input files.\n";
+        std::cout << errorMessageOut
+                  << "\n";
+    }
+}
+
+int main(int argc, char** argv) try
+{
+    using namespace Dumux;
+
+    // define the type tag for this problem
+    using TypeTag = TTAG(PipeLauferProblem);
+
+    // 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, usage);
+
+    // 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);
+
+    // 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");
+
+    // check if we are about to restart a previously interrupted simulation
+    Scalar restartTime = 0;
+    if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+        restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+    // instantiate time loop
+    auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
+    timeLoop->setMaxTimeStepSize(maxDt);
+    problem->setTimeLoop(timeLoop);
+
+    // the solution vector
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+    using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+    typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+    typename DofTypeIndices::FaceIdx faceIdx;
+    const auto numDofsCellCenter = leafGridView.size(0);
+    const auto numDofsFace = leafGridView.size(1);
+    SolutionVector x;
+    x[cellCenterIdx].resize(numDofsCellCenter);
+    x[faceIdx].resize(numDofsFace);
+    problem->applyInitialSolution(x);
+    problem->updateDynamicWallProperties(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);
+
+    // intialize the vtk output module
+    using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+    StaggeredVtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+    VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
+    vtkWriter.write(0.0);
+
+    // the assembler with time loop for instationary problem
+    using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
+    auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+    // the linear solver
+    using LinearSolver = Dumux::UMFPackBackend;
+    auto linearSolver = std::make_shared<LinearSolver>();
+
+    // the non-linear solver
+    using NewtonSolver = Dumux::NewtonSolver<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();
+
+        // update wall properties
+        problem->updateDynamicWallProperties(x);
+
+        // 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 solver
+        timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+    } while (!timeLoop->finished());
+
+    timeLoop->finalize(leafGridView.comm());
+
+    ////////////////////////////////////////////////////////////
+    // finalize, print dumux message to say goodbye
+    ////////////////////////////////////////////////////////////
+
+#if HAVE_PVPYTHON
+    bool shouldPlot = getParam<bool>("Output.PlotVelocityProfiles", false);
+    if (shouldPlot)
+    {
+        char fileName[255];
+        std::string fileNameFormat = "%s-%05d";
+        sprintf(fileName, fileNameFormat.c_str(), problem->name().c_str(), timeLoop->timeStepIndex());
+        std::cout << fileName << std::endl;
+        std::string vtuFileName = std::string(fileName) + ".vtu";
+        std::string script = std::string(DUMUX_SOURCE_DIR) + "/bin/postprocessing/extractlinedata.py";
+        std::string syscom;
+
+        // execute the pvpython script
+        std::string command = std::string(PVPYTHON_EXECUTABLE) + " " + script
+                              + " -f " + vtuFileName
+                              + " -v 1"
+                              + " -r 10000";
+        syscom =  command + " -p1 8.0 0.0 0.0"
+                          + " -p2 8.0 0.2469 0.0"
+                          + " -of " + std::string(fileName) + "\n";
+        system(syscom.c_str());
+
+
+        char gnuplotFileName[255];
+        Dumux::GnuplotInterface<Scalar> gnuplot;
+        sprintf(gnuplotFileName, fileNameFormat.c_str(), "velProfiles", timeLoop->timeStepIndex());
+        gnuplot.setOpenPlotWindow(shouldPlot);
+        gnuplot.setDatafileSeparator(',');
+        gnuplot.resetPlot();
+        gnuplot.setXlabel("y^+ [-]");
+        gnuplot.setYlabel("u_+ [-]");
+        gnuplot.setOption("set log x");
+        gnuplot.setOption("set xrange [1:3000]");
+        gnuplot.addFileToPlot("laufer_re50000_u+y+.csv", "u 1:2 w p t 'Laufer1954a, Re=50000'");
+        gnuplot.addFileToPlot(std::string(fileName) + ".csv", "u 13:14 w l");
+        gnuplot.plot(std::string(gnuplotFileName));
+    }
+#endif
+
+    // print dumux end message
+    if (mpiHelper.rank() == 0)
+    {
+        Parameters::print();
+        DumuxMessage::print(/*firstCall=*/false);
+    }
+
+    return 0;
+} // end main
+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/freeflow/rans/test_pipe_laufer.input b/test/freeflow/rans/test_pipe_laufer.input
new file mode 100644
index 0000000000000000000000000000000000000000..f507011c757726e58be6a7a942706b3abf58f5ca
--- /dev/null
+++ b/test/freeflow/rans/test_pipe_laufer.input
@@ -0,0 +1,30 @@
+[TimeLoop]
+DtInitial = 1e-1 # [s]
+TEnd = 100 # [s]
+
+[Grid]
+Verbosity = true
+Positions0 = 0.0 10.0
+Positions1 = 0.0 0.12345 0.2469
+Cells0 = 25
+Cells1 = 25 25
+Grading1 = 1.2 -1.2
+
+[Output]
+PlotVelocityProfiles = true
+
+[Problem]
+Name = pipe_laufer_zeroeq
+InletVelocity = 2.5 # [m/s]
+EnableGravity = false
+
+[RANS]
+EddyViscosityModel = 3
+
+[Newton]
+MaxSteps = 10
+MaxRelativeShift = 1e-5
+
+[Vtk]
+AddVelocity = true
+WriteFaceData = false
diff --git a/test/freeflow/rans/test_pipe_laufer_reference.input b/test/freeflow/rans/test_pipe_laufer_reference.input
new file mode 100644
index 0000000000000000000000000000000000000000..e33bbd4ca2f11f87d25f40a3033eeb3f9862fb88
--- /dev/null
+++ b/test/freeflow/rans/test_pipe_laufer_reference.input
@@ -0,0 +1,27 @@
+[TimeLoop]
+DtInitial = 1e-1 # [s]
+TEnd = 100 # [s]
+
+[Grid]
+Verbosity = true
+Positions0 = 0.0 10.0
+Positions1 = 0.0 0.12345 0.2469
+Cells0 = 16
+Cells1 = 8 8
+Grading1 = 1.4 -1.4
+
+[Problem]
+Name = pipe_laufer_zeroeq_reference
+InletVelocity = 2.5 # [m/s]
+EnableGravity = false
+
+[RANS]
+EddyViscosityModel = 3
+
+[Newton]
+MaxSteps = 10
+MaxRelativeShift = 1e-5
+
+[Vtk]
+AddVelocity = true
+WriteFaceData = false
diff --git a/test/references/pipe_laufer_zeroeq.vtu b/test/references/pipe_laufer_zeroeq.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..c16176ac3710caebb5f7af8731c53ef288c6d943
--- /dev/null
+++ b/test/references/pipe_laufer_zeroeq.vtu
@@ -0,0 +1,578 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+  <UnstructuredGrid>
+    <Piece NumberOfCells="256" NumberOfPoints="289">
+      <CellData Scalars="v_x [m/s]" Vectors="velocity_Air (m/s)">
+        <DataArray type="Float32" Name="v_x [m/s]" NumberOfComponents="1" format="ascii">
+          2.07761 1.47958 1.28367 1.25467 1.2256 1.19748 1.1905 1.18952 1.1839 1.17538 1.16979 1.17192
+          1.17724 1.18028 1.18099 1.18075 2.38628 2.14543 1.94733 1.84353 1.78739 1.76159 1.76636 1.77234
+          1.76679 1.75639 1.74958 1.7531 1.75962 1.76261 1.76304 1.76251 2.48263 2.40526 2.27794 2.16452
+          2.08742 2.04048 2.01824 2.00498 1.99072 1.97588 1.97047 1.97997 1.98992 1.99312 1.99343 1.99278
+          2.51407 2.52122 2.48153 2.41111 2.34207 2.28411 2.23834 2.20558 2.18118 2.16111 2.15029 2.1515
+          2.15494 2.15515 2.15383 2.15234 2.52057 2.5554 2.56969 2.55829 2.53194 2.49793 2.45943 2.42442
+          2.39621 2.37341 2.35442 2.33973 2.33042 2.32449 2.31988 2.31677 2.52129 2.56089 2.59214 2.61312
+          2.6256 2.62815 2.61908 2.6043 2.58925 2.57518 2.56046 2.54427 2.53032 2.52023 2.51276 2.50816
+          2.52126 2.56117 2.59479 2.62329 2.64923 2.67126 2.68737 2.69894 2.70788 2.71458 2.71691 2.71332
+          2.70763 2.70291 2.69919 2.69686 2.52119 2.56104 2.59478 2.62384 2.65108 2.67635 2.69959 2.72192
+          2.74377 2.76466 2.78269 2.79612 2.80658 2.81594 2.82425 2.83 2.52119 2.56104 2.59478 2.62384
+          2.65108 2.67635 2.69959 2.72192 2.74377 2.76466 2.78269 2.79612 2.80658 2.81594 2.82425 2.83
+          2.52126 2.56117 2.59479 2.62329 2.64923 2.67126 2.68737 2.69894 2.70788 2.71458 2.71691 2.71332
+          2.70763 2.70291 2.69919 2.69686 2.52129 2.56089 2.59214 2.61312 2.6256 2.62815 2.61908 2.6043
+          2.58925 2.57518 2.56046 2.54427 2.53032 2.52023 2.51276 2.50816 2.52057 2.5554 2.56969 2.55829
+          2.53194 2.49793 2.45943 2.42442 2.39621 2.37341 2.35442 2.33973 2.33042 2.32449 2.31988 2.31677
+          2.51407 2.52122 2.48153 2.41111 2.34207 2.28411 2.23834 2.20558 2.18118 2.16111 2.15029 2.1515
+          2.15494 2.15515 2.15383 2.15234 2.48263 2.40526 2.27794 2.16452 2.08742 2.04048 2.01824 2.00498
+          1.99072 1.97588 1.97047 1.97997 1.98992 1.99312 1.99343 1.99278 2.38628 2.14543 1.94733 1.84353
+          1.78739 1.76159 1.76636 1.77234 1.76679 1.75639 1.74958 1.7531 1.75962 1.76261 1.76304 1.76251
+          2.07761 1.47958 1.28367 1.25467 1.2256 1.19748 1.1905 1.18952 1.1839 1.17538 1.16979 1.17192
+          1.17724 1.18028 1.18099 1.18075
+        </DataArray>
+        <DataArray type="Float32" Name="dv_x/dx [m/s]" NumberOfComponents="1" format="ascii">
+          -0.95686 -0.635157 -0.179922 -0.0464533 -0.0457548 -0.0280812 -0.00636892 -0.00527861 -0.0113137 -0.0112903 -0.0027612 0.0059661
+          0.00668688 0.00299378 0.00037709 -0.000371439 -0.385348 -0.35116 -0.241522 -0.127953 -0.0655552 -0.0168229 0.00859883 0.000348277
+          -0.0127609 -0.0137686 -0.00262868 0.00803264 0.00760658 0.00273242 -7.92547e-05 -0.00084692 -0.123792 -0.163751 -0.192587 -0.152414
+          -0.0992347 -0.0553436 -0.0283987 -0.0220144 -0.0232791 -0.0162027 0.00326662 0.0155584 0.0105211 0.00280934 -0.000268182 -0.00103498
+          0.0114385 -0.0260362 -0.0880861 -0.111566 -0.101605 -0.0829837 -0.0628222 -0.0457291 -0.0355761 -0.024713 -0.00768423 0.00372169
+          0.00291362 -0.000891487 -0.00224355 -0.00237188 0.0557332 0.0392977 0.00230937 -0.0301996 -0.0482888 -0.0580147 -0.0588047 -0.0505745
+          -0.0408131 -0.0334319 -0.0269454 -0.0191991 -0.0121905 -0.00843178 -0.00617655 -0.00497851 0.0633478 0.056677 0.0417843 0.0267666
+          0.0120249 -0.00521208 -0.0190776 -0.0238655 -0.0232978 -0.0230369 -0.0247241 -0.0241065 -0.0192386 -0.0140473 -0.00964972 -0.00736024
+          0.0638551 0.0588229 0.0496947 0.0435467 0.038378 0.0305152 0.0221386 0.0164115 0.0125123 0.00721903 -0.00100371 -0.00742319
+          -0.0083278 -0.00675228 -0.00484026 -0.00372104 0.0637564 0.0588739 0.0502416 0.0450384 0.0420084 0.038806 0.0364566 0.0353462
+          0.0341902 0.0311362 0.0251649 0.0191133 0.0158566 0.0141348 0.0112495 0.00919434 0.0637564 0.0588739 0.0502416 0.0450384
+          0.0420084 0.038806 0.0364566 0.0353462 0.0341902 0.0311362 0.0251649 0.0191133 0.0158566 0.0141348 0.0112495 0.00919434
+          0.0638551 0.0588229 0.0496947 0.0435467 0.038378 0.0305152 0.0221386 0.0164115 0.0125123 0.00721903 -0.00100371 -0.00742319
+          -0.0083278 -0.00675228 -0.00484026 -0.00372104 0.0633478 0.056677 0.0417843 0.0267666 0.0120249 -0.00521208 -0.0190776 -0.0238655
+          -0.0232978 -0.0230369 -0.0247241 -0.0241065 -0.0192386 -0.0140473 -0.00964972 -0.00736024 0.0557332 0.0392977 0.00230937 -0.0301996
+          -0.0482888 -0.0580147 -0.0588047 -0.0505745 -0.0408131 -0.0334319 -0.0269454 -0.0191991 -0.0121905 -0.00843178 -0.00617655 -0.00497851
+          0.0114385 -0.0260362 -0.0880861 -0.111566 -0.101605 -0.0829837 -0.0628222 -0.0457291 -0.0355761 -0.024713 -0.00768423 0.00372169
+          0.00291362 -0.000891487 -0.00224355 -0.00237188 -0.123792 -0.163751 -0.192587 -0.152414 -0.0992347 -0.0553436 -0.0283987 -0.0220144
+          -0.0232791 -0.0162027 0.00326662 0.0155584 0.0105211 0.00280934 -0.000268182 -0.00103498 -0.385348 -0.35116 -0.241522 -0.127953
+          -0.0655552 -0.0168229 0.00859883 0.000348277 -0.0127609 -0.0137686 -0.00262868 0.00803264 0.00760658 0.00273242 -7.92547e-05 -0.00084692
+          -0.95686 -0.635157 -0.179922 -0.0464533 -0.0457548 -0.0280812 -0.00636892 -0.00527861 -0.0113137 -0.0112903 -0.0027612 0.0059661
+          0.00668688 0.00299378 0.00037709 -0.000371439
+        </DataArray>
+        <DataArray type="Float32" Name="dv_x/dy [m/s]" NumberOfComponents="1" format="ascii">
+          71.6647 154.597 154.087 136.72 130.434 130.973 133.701 135.317 135.334 134.897 134.615 134.935
+          135.215 135.202 135.139 135.07 39.1814 89.551 96.1864 88.0195 83.3731 81.5524 80.0763 78.8884
+          78.0526 77.4415 77.4585 78.1704 78.6185 78.6341 78.5962 78.5561 8.83074 25.9671 36.9134 39.2201
+          38.3288 36.1063 32.6142 29.9373 28.6342 27.9666 27.6889 27.53 27.3165 27.1246 27.0035 26.9377
+          1.87278 7.41083 14.4002 19.4353 21.9405 22.5785 21.7757 20.7025 20.0136 19.6207 18.9506 17.7568
+          16.8062 16.3555 16.1127 15.991 0.254612 1.39841 3.89971 7.12168 9.99591 12.1292 13.4232 14.057
+          14.3867 14.5981 14.4606 13.8472 13.2342 12.8709 12.6545 12.5446 0.0174282 0.145258 0.632041 1.63685
+          2.95342 4.36494 5.74016 6.91287 7.84874 8.59144 9.12835 9.40803 9.49903 9.52962 9.5519 9.57169
+          -0.00186854 0.00272496 0.0475259 0.19288 0.458353 0.867036 1.44807 2.11568 2.7794 3.40825 3.99744 4.52997
+          4.9692 5.31906 5.60285 5.78896 -0.00103646 -0.0019143 -0.000129744 0.00781458 0.0264059 0.0723988 0.173897 0.327112
+          0.51074 0.712754 0.93622 1.17829 1.4083 1.60852 1.77987 1.89475 0.00103646 0.0019143 0.000129744 -0.00781458
+          -0.0264059 -0.0723988 -0.173897 -0.327112 -0.51074 -0.712754 -0.93622 -1.17829 -1.4083 -1.60852 -1.77987 -1.89475
+          0.00186854 -0.00272496 -0.0475259 -0.19288 -0.458353 -0.867036 -1.44807 -2.11568 -2.7794 -3.40825 -3.99744 -4.52997
+          -4.9692 -5.31906 -5.60285 -5.78896 -0.0174282 -0.145258 -0.632041 -1.63685 -2.95342 -4.36494 -5.74016 -6.91287
+          -7.84874 -8.59144 -9.12835 -9.40803 -9.49903 -9.52962 -9.5519 -9.57169 -0.254612 -1.39841 -3.89971 -7.12168
+          -9.99591 -12.1292 -13.4232 -14.057 -14.3867 -14.5981 -14.4606 -13.8472 -13.2342 -12.8709 -12.6545 -12.5446
+          -1.87278 -7.41083 -14.4002 -19.4353 -21.9405 -22.5785 -21.7757 -20.7025 -20.0136 -19.6207 -18.9506 -17.7568
+          -16.8062 -16.3555 -16.1127 -15.991 -8.83074 -25.9671 -36.9134 -39.2201 -38.3288 -36.1063 -32.6142 -29.9373
+          -28.6342 -27.9666 -27.6889 -27.53 -27.3165 -27.1246 -27.0035 -26.9377 -39.1814 -89.551 -96.1864 -88.0195
+          -83.3731 -81.5524 -80.0763 -78.8884 -78.0526 -77.4415 -77.4585 -78.1704 -78.6185 -78.6341 -78.5962 -78.5561
+          -71.6647 -154.597 -154.087 -136.72 -130.434 -130.973 -133.701 -135.317 -135.334 -134.897 -134.615 -134.935
+          -135.215 -135.202 -135.139 -135.07
+        </DataArray>
+        <DataArray type="Float32" Name="p [Pa]" NumberOfComponents="1" format="ascii">
+          100001 100001 100001 100001 100001 100001 100001 100001 100001 100000 100000 100000
+          100000 100000 100000 100000 100001 100001 100001 100001 100001 100001 100001 100001
+          100001 100000 100000 100000 100000 100000 100000 100000 100001 100001 100001 100001
+          100001 100001 100001 100001 100001 100000 100000 100000 100000 100000 100000 100000
+          100001 100001 100001 100001 100001 100001 100001 100001 100001 100000 100000 100000
+          100000 100000 100000 100000 100001 100001 100001 100001 100001 100001 100001 100001
+          100001 100000 100000 100000 100000 100000 100000 100000 100001 100001 100001 100001
+          100001 100001 100001 100001 100001 100000 100000 100000 100000 100000 100000 100000
+          100001 100001 100001 100001 100001 100001 100001 100001 100001 100000 100000 100000
+          100000 100000 100000 100000 100001 100001 100001 100001 100001 100001 100001 100001
+          100001 100000 100000 100000 100000 100000 100000 100000 100001 100001 100001 100001
+          100001 100001 100001 100001 100001 100000 100000 100000 100000 100000 100000 100000
+          100001 100001 100001 100001 100001 100001 100001 100001 100001 100000 100000 100000
+          100000 100000 100000 100000 100001 100001 100001 100001 100001 100001 100001 100001
+          100001 100000 100000 100000 100000 100000 100000 100000 100001 100001 100001 100001
+          100001 100001 100001 100001 100001 100000 100000 100000 100000 100000 100000 100000
+          100001 100001 100001 100001 100001 100001 100001 100001 100001 100000 100000 100000
+          100000 100000 100000 100000 100001 100001 100001 100001 100001 100001 100001 100001
+          100001 100000 100000 100000 100000 100000 100000 100000 100001 100001 100001 100001
+          100001 100001 100001 100001 100001 100000 100000 100000 100000 100000 100000 100000
+          100001 100001 100001 100001 100001 100001 100001 100001 100001 100000 100000 100000
+          100000 100000 100000 100000
+        </DataArray>
+        <DataArray type="Float32" Name="p_rel [Pa]" NumberOfComponents="1" format="ascii">
+          1.30543 1.17114 1.0534 0.957444 0.867397 0.780484 0.700841 0.624126 0.545742 0.467086 0.389893 0.318209
+          0.249301 0.178674 0.107209 0.0356539 1.30533 1.17127 1.05352 0.957484 0.867414 0.780489 0.700836 0.624126
+          0.545745 0.467091 0.389894 0.318207 0.249299 0.178673 0.107209 0.0356541 1.30511 1.17134 1.05363 0.957553
+          0.867444 0.780539 0.700877 0.62414 0.54576 0.467106 0.389892 0.318194 0.249292 0.178671 0.107209 0.0356545
+          1.30481 1.1714 1.05374 0.957619 0.867464 0.780612 0.700949 0.624149 0.545765 0.467108 0.38993 0.318224
+          0.249288 0.178669 0.107209 0.0356557 1.30443 1.17146 1.05384 0.957674 0.867482 0.780689 0.701023 0.624157
+          0.545771 0.46711 0.390012 0.318306 0.249288 0.178669 0.10721 0.0356575 1.30397 1.17153 1.05393 0.957711
+          0.867497 0.780761 0.70109 0.624161 0.545775 0.467112 0.390105 0.318397 0.249287 0.178668 0.107209 0.0356632
+          1.30348 1.17159 1.05403 0.95775 0.867524 0.780819 0.701128 0.624161 0.545776 0.467117 0.390202 0.318491
+          0.249283 0.178665 0.107208 0.0356817 1.30308 1.17165 1.0541 0.957778 0.867543 0.78085 0.701155 0.624172
+          0.54579 0.467135 0.390273 0.31855 0.249282 0.178664 0.107208 0.0357212 1.30308 1.17165 1.0541 0.957778
+          0.867543 0.78085 0.701155 0.624172 0.54579 0.467135 0.390273 0.31855 0.249282 0.178664 0.107208 0.0357212
+          1.30348 1.17159 1.05403 0.95775 0.867524 0.780819 0.701128 0.624161 0.545776 0.467117 0.390202 0.318491
+          0.249283 0.178665 0.107208 0.0356817 1.30397 1.17153 1.05393 0.957711 0.867497 0.780761 0.70109 0.624161
+          0.545775 0.467112 0.390105 0.318397 0.249287 0.178668 0.107209 0.0356632 1.30443 1.17146 1.05384 0.957674
+          0.867482 0.780689 0.701023 0.624157 0.545771 0.46711 0.390012 0.318306 0.249288 0.178669 0.10721 0.0356575
+          1.30481 1.1714 1.05374 0.957619 0.867464 0.780612 0.700949 0.624149 0.545765 0.467108 0.38993 0.318224
+          0.249288 0.178669 0.107209 0.0356557 1.30511 1.17134 1.05363 0.957553 0.867444 0.780539 0.700877 0.62414
+          0.54576 0.467106 0.389892 0.318194 0.249292 0.178671 0.107209 0.0356545 1.30533 1.17127 1.05352 0.957484
+          0.867414 0.780489 0.700836 0.624126 0.545745 0.467091 0.389894 0.318207 0.249299 0.178673 0.107209 0.0356541
+          1.30543 1.17114 1.0534 0.957444 0.867397 0.780484 0.700841 0.624126 0.545742 0.467086 0.389893 0.318209
+          0.249301 0.178674 0.107209 0.0356539
+        </DataArray>
+        <DataArray type="Float32" Name="rho [kg/m^3]" NumberOfComponents="1" format="ascii">
+          1.23014 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23012 1.23012
+          1.23012 1.23012 1.23012 1.23012 1.23014 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013
+          1.23013 1.23013 1.23012 1.23012 1.23012 1.23012 1.23012 1.23012 1.23014 1.23013 1.23013 1.23013
+          1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23012 1.23012 1.23012 1.23012 1.23012 1.23012
+          1.23014 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23012 1.23012
+          1.23012 1.23012 1.23012 1.23012 1.23014 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013
+          1.23013 1.23013 1.23012 1.23012 1.23012 1.23012 1.23012 1.23012 1.23014 1.23013 1.23013 1.23013
+          1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23012 1.23012 1.23012 1.23012 1.23012 1.23012
+          1.23014 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23012 1.23012
+          1.23012 1.23012 1.23012 1.23012 1.23014 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013
+          1.23013 1.23013 1.23012 1.23012 1.23012 1.23012 1.23012 1.23012 1.23014 1.23013 1.23013 1.23013
+          1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23012 1.23012 1.23012 1.23012 1.23012 1.23012
+          1.23014 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23012 1.23012
+          1.23012 1.23012 1.23012 1.23012 1.23014 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013
+          1.23013 1.23013 1.23012 1.23012 1.23012 1.23012 1.23012 1.23012 1.23014 1.23013 1.23013 1.23013
+          1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23012 1.23012 1.23012 1.23012 1.23012 1.23012
+          1.23014 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23012 1.23012
+          1.23012 1.23012 1.23012 1.23012 1.23014 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013
+          1.23013 1.23013 1.23012 1.23012 1.23012 1.23012 1.23012 1.23012 1.23014 1.23013 1.23013 1.23013
+          1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23012 1.23012 1.23012 1.23012 1.23012 1.23012
+          1.23014 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23013 1.23012 1.23012
+          1.23012 1.23012 1.23012 1.23012
+        </DataArray>
+        <DataArray type="Float32" Name="nu [m^2/s]" NumberOfComponents="1" format="ascii">
+          1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05
+          1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05
+          1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05
+          1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05
+          1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05
+          1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05
+          1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05
+          1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05
+          1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05
+          1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05
+          1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05
+          1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05
+          1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05
+          1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05
+          1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05
+          1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05
+          1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05
+          1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05
+          1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05
+          1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05
+          1.43726e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43727e-05 1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05
+          1.43728e-05 1.43728e-05 1.43728e-05 1.43728e-05
+        </DataArray>
+        <DataArray type="Float32" Name="nu_t [m^2/s]" NumberOfComponents="1" format="ascii">
+          7.91604e-07 3.43486e-06 3.41343e-06 2.72149e-06 2.48886e-06 2.50844e-06 2.60859e-06 2.66875e-06 2.66939e-06 2.65304e-06 2.64253e-06 2.65448e-06
+          2.66491e-06 2.66445e-06 2.66207e-06 2.65949e-06 1.59309e-05 4.79202e-05 0.000114341 0.000146157 0.000134052 0.000131497 0.000130955 0.000130076
+          0.000128709 0.00012742 0.000127266 0.000128644 0.000129564 0.000129581 0.000129477 0.000129366 1.28196e-05 3.85614e-05 0.000113892 0.000201837
+          0.00022658 0.000233286 0.000370124 0.000387786 0.000396884 0.000400233 0.000395888 0.000394033 0.000391339 0.000388574 0.000386758 0.000385727
+          2.32091e-06 6.9813e-06 0.000104392 0.000201065 0.000225713 0.000232393 0.000370031 0.000387688 0.000396784 0.000402569 0.00039875 0.000603762
+          0.000609585 0.00061155 0.000612984 0.000614257 1.76256e-07 5.30179e-07 4.66189e-05 0.000190251 0.000213574 0.000219895 0.00036865 0.000386242
+          0.000395303 0.000401067 0.000397262 0.000603572 0.000609393 0.000611357 0.000612791 0.000614064 1.49726e-08 4.50377e-08 6.25642e-06 0.000116937
+          0.000131272 0.000135157 0.000353284 0.000370142 0.000378825 0.000384349 0.000380703 0.000601364 0.000607164 0.000609121 0.000610549 0.000611818
+          1.46502e-09 4.40678e-09 6.43413e-07 2.39486e-05 2.68843e-05 2.768e-05 0.00024879 0.000260662 0.000266777 0.000270667 0.0002681 0.000579989
+          0.000585583 0.00058747 0.000588847 0.000590071 1.57785e-10 4.74619e-10 6.96409e-08 2.88443e-06 3.23803e-06 3.33386e-06 6.69465e-05 7.01411e-05
+          7.17866e-05 7.28334e-05 7.21424e-05 0.000437258 0.000441475 0.000442897 0.000443936 0.000444858 1.57785e-10 4.74619e-10 6.96409e-08 2.88443e-06
+          3.23803e-06 3.33386e-06 6.69465e-05 7.01411e-05 7.17866e-05 7.28334e-05 7.21424e-05 0.000437258 0.000441475 0.000442897 0.000443936 0.000444858
+          1.46502e-09 4.40678e-09 6.43413e-07 2.39486e-05 2.68843e-05 2.768e-05 0.00024879 0.000260662 0.000266777 0.000270667 0.0002681 0.000579989
+          0.000585583 0.00058747 0.000588847 0.000590071 1.49726e-08 4.50377e-08 6.25642e-06 0.000116937 0.000131272 0.000135157 0.000353284 0.000370142
+          0.000378825 0.000384349 0.000380703 0.000601364 0.000607164 0.000609121 0.000610549 0.000611818 1.76256e-07 5.30179e-07 4.66189e-05 0.000190251
+          0.000213574 0.000219895 0.00036865 0.000386242 0.000395303 0.000401067 0.000397262 0.000603572 0.000609393 0.000611357 0.000612791 0.000614064
+          2.32091e-06 6.9813e-06 0.000104392 0.000201065 0.000225713 0.000232393 0.000370031 0.000387688 0.000396784 0.000402569 0.00039875 0.000603762
+          0.000609585 0.00061155 0.000612984 0.000614257 1.28196e-05 3.85614e-05 0.000113892 0.000201837 0.00022658 0.000233286 0.000370124 0.000387786
+          0.000396884 0.000400233 0.000395888 0.000394033 0.000391339 0.000388574 0.000386758 0.000385727 1.59309e-05 4.79202e-05 0.000114341 0.000146157
+          0.000134052 0.000131497 0.000130955 0.000130076 0.000128709 0.00012742 0.000127266 0.000128644 0.000129564 0.000129581 0.000129477 0.000129366
+          7.91604e-07 3.43486e-06 3.41343e-06 2.72149e-06 2.48886e-06 2.50844e-06 2.60859e-06 2.66875e-06 2.66939e-06 2.65304e-06 2.64253e-06 2.65448e-06
+          2.66491e-06 2.66445e-06 2.66207e-06 2.65949e-06
+        </DataArray>
+        <DataArray type="Float32" Name="l_w [m]" NumberOfComponents="1" format="ascii">
+          0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461
+          0.00179461 0.00179461 0.00179461 0.00179461 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166
+          0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.0121315 0.0121315 0.0121315 0.0121315
+          0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315
+          0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734
+          0.0205734 0.0205734 0.0205734 0.0205734 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919
+          0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0489379 0.0489379 0.0489379 0.0489379
+          0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379
+          0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023
+          0.0721023 0.0721023 0.0721023 0.0721023 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532
+          0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532
+          0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532 0.104532
+          0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023 0.0721023
+          0.0721023 0.0721023 0.0721023 0.0721023 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379
+          0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0489379 0.0323919 0.0323919 0.0323919 0.0323919
+          0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919 0.0323919
+          0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734 0.0205734
+          0.0205734 0.0205734 0.0205734 0.0205734 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315
+          0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.0121315 0.00610166 0.00610166 0.00610166 0.00610166
+          0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166 0.00610166
+          0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461 0.00179461
+          0.00179461 0.00179461 0.00179461 0.00179461
+        </DataArray>
+        <DataArray type="Float32" Name="y^+ [-]" NumberOfComponents="1" format="ascii">
+          4.00731 5.88574 5.87602 5.53497 5.40624 5.4174 5.47353 5.5065 5.50685 5.49794 5.49219 5.49872
+          5.50441 5.50415 5.50286 5.50145 13.6249 20.0115 19.9785 18.8189 18.3812 18.4192 18.61 18.7221
+          18.7233 18.693 18.6735 18.6957 18.715 18.7141 18.7097 18.7049 27.0894 39.7876 39.7219 37.4164
+          36.5462 36.6216 37.001 37.2239 37.2263 37.1661 37.1272 37.1714 37.2098 37.2081 37.1993 37.1898
+          45.9399 67.4742 67.3627 63.4529 61.9771 62.105 62.7485 63.1265 63.1305 63.0284 62.9625 63.0373
+          63.1025 63.0996 63.0848 63.0687 72.3304 106.235 106.06 99.904 97.5804 97.7819 98.795 99.3901
+          99.3964 99.2357 99.1319 99.2497 99.3524 99.3478 99.3244 99.299 109.277 160.501 160.236 150.936
+          147.425 147.729 149.26 150.159 150.169 149.926 149.769 149.947 150.102 150.095 150.06 150.022
+          161.003 236.473 236.082 222.38 217.208 217.656 219.911 221.236 221.25 220.892 220.661 220.923
+          221.152 221.142 221.09 221.033 233.418 342.833 342.267 322.402 314.903 315.553 318.823 320.743
+          320.763 320.245 319.91 320.29 320.621 320.607 320.531 320.449 233.418 342.833 342.267 322.402
+          314.903 315.553 318.823 320.743 320.763 320.245 319.91 320.29 320.621 320.607 320.531 320.449
+          161.003 236.473 236.082 222.38 217.208 217.656 219.911 221.236 221.25 220.892 220.661 220.923
+          221.152 221.142 221.09 221.033 109.277 160.501 160.236 150.936 147.425 147.729 149.26 150.159
+          150.169 149.926 149.769 149.947 150.102 150.095 150.06 150.022 72.3304 106.235 106.06 99.904
+          97.5804 97.7819 98.795 99.3901 99.3964 99.2357 99.1319 99.2497 99.3524 99.3478 99.3244 99.299
+          45.9399 67.4742 67.3627 63.4529 61.9771 62.105 62.7485 63.1265 63.1305 63.0284 62.9625 63.0373
+          63.1025 63.0996 63.0848 63.0687 27.0894 39.7876 39.7219 37.4164 36.5462 36.6216 37.001 37.2239
+          37.2263 37.1661 37.1272 37.1714 37.2098 37.2081 37.1993 37.1898 13.6249 20.0115 19.9785 18.8189
+          18.3812 18.4192 18.61 18.7221 18.7233 18.693 18.6735 18.6957 18.715 18.7141 18.7097 18.7049
+          4.00731 5.88574 5.87602 5.53497 5.40624 5.4174 5.47353 5.5065 5.50685 5.49794 5.49219 5.49872
+          5.50441 5.50415 5.50286 5.50145
+        </DataArray>
+        <DataArray type="Float32" Name="u^+ [-]" NumberOfComponents="1" format="ascii">
+          64.7357 31.3883 27.2772 28.3039 28.3064 27.5999 27.1576 26.9727 26.8436 26.6935 26.5944 26.6113
+          26.7045 26.7746 26.7969 26.7984 74.3532 45.5141 41.3797 41.5878 41.2813 40.6017 40.2941 40.1883
+          40.0601 39.8886 39.7756 39.8082 39.915 39.9846 40.0037 40.0019 77.3554 51.0261 48.405 48.829
+          48.2109 47.0297 46.0401 45.4637 45.1375 44.8735 44.7974 44.9598 45.1391 45.2137 45.2314 45.2282
+          78.3352 53.4862 52.7312 54.3918 54.0923 52.6449 51.0609 50.0123 49.4558 49.0801 48.8854 48.855
+          48.8824 48.8893 48.8708 48.8496 78.5377 54.2114 54.6046 57.7119 58.4776 57.5731 56.1043 54.9747
+          54.3313 53.9015 53.5262 53.129 52.8629 52.7308 52.6386 52.5814 78.5602 54.3277 55.0816 58.9487
+          60.6407 60.5744 59.7464 59.0534 58.7084 58.4839 58.2103 57.7737 57.3975 57.1711 57.0152 56.9253
+          78.5592 54.3337 55.138 59.1782 61.1863 61.5682 61.3042 61.1993 61.3983 61.6497 61.7672 61.6124
+          61.4195 61.3153 61.2452 61.208 78.557 54.3309 55.1378 59.1906 61.2292 61.6854 61.5829 61.7205
+          62.212 62.7871 63.2627 63.4924 63.6642 63.8792 64.0829 64.2297 78.557 54.3309 55.1378 59.1906
+          61.2292 61.6854 61.5829 61.7205 62.212 62.7871 63.2627 63.4924 63.6642 63.8792 64.0829 64.2297
+          78.5592 54.3337 55.138 59.1782 61.1863 61.5682 61.3042 61.1993 61.3983 61.6497 61.7672 61.6124
+          61.4195 61.3153 61.2452 61.208 78.5602 54.3277 55.0816 58.9487 60.6407 60.5744 59.7464 59.0534
+          58.7084 58.4839 58.2103 57.7737 57.3975 57.1711 57.0152 56.9253 78.5377 54.2114 54.6046 57.7119
+          58.4776 57.5731 56.1043 54.9747 54.3313 53.9015 53.5262 53.129 52.8629 52.7308 52.6386 52.5814
+          78.3352 53.4862 52.7312 54.3918 54.0923 52.6449 51.0609 50.0123 49.4558 49.0801 48.8854 48.855
+          48.8824 48.8893 48.8708 48.8496 77.3554 51.0261 48.405 48.829 48.2109 47.0297 46.0401 45.4637
+          45.1375 44.8735 44.7974 44.9598 45.1391 45.2137 45.2314 45.2282 74.3532 45.5141 41.3797 41.5878
+          41.2813 40.6017 40.2941 40.1883 40.0601 39.8886 39.7756 39.8082 39.915 39.9846 40.0037 40.0019
+          64.7357 31.3883 27.2772 28.3039 28.3064 27.5999 27.1576 26.9727 26.8436 26.6935 26.5944 26.6113
+          26.7045 26.7746 26.7969 26.7984
+        </DataArray>
+        <DataArray type="Float32" Name="velocity_Air (m/s)" NumberOfComponents="3" format="ascii">
+          2.07761 0.00242566 0 1.47957 0.00100872 0 1.28367 0.000116336 0 1.25468 5.01453e-05 0
+          1.22561 0.000116787 0 1.1975 4.46784e-05 0 1.19052 -4.59293e-06 0 1.18953 1.02362e-05 0
+          1.18391 2.20593e-05 0 1.17538 2.6933e-05 0 1.16978 5.26205e-06 0 1.17188 -1.73307e-05 0
+          1.17714 -1.28878e-05 0 1.18011 -4.15204e-06 0 1.18075 5.06639e-07 0 1.18047 1.08101e-06 0
+          2.38628 0.00576564 0 2.14543 0.00303948 0 1.94733 0.000803422 0 1.84352 0.000364102 0
+          1.78737 0.000421259 0 1.76156 0.000109178 0 1.76632 -6.72728e-05 0 1.7723 3.05429e-05 0
+          1.76675 7.86469e-05 0 1.75635 0.00010299 0 1.74955 1.60399e-05 0 1.75309 -6.85792e-05 0
+          1.75963 -4.45018e-05 0 1.76266 -1.38629e-05 0 1.76313 2.77535e-06 0 1.76263 4.39272e-06 0
+          2.48263 0.00687549 0 2.40526 0.00473685 0 2.27794 0.00213196 0 2.16452 0.00114677 0
+          2.08741 0.000958049 0 2.04046 0.00030835 0 2.01822 -5.431e-05 0 2.00495 0.000118866 0
+          1.99069 0.000195461 0 1.97585 0.000236882 0 1.97045 -2.33479e-06 0 1.97996 -0.000185671 0
+          1.98994 -9.23723e-05 0 1.99318 -2.67383e-05 0 1.99354 7.8219e-06 0 1.99293 1.0217e-05 0
+          2.51407 0.00684927 0 2.52122 0.00552131 0 2.48153 0.00340618 0 2.41111 0.00225882 0
+          2.34207 0.00180202 0 2.2841 0.000906265 0 2.23833 0.000319485 0 2.20557 0.000410678 0
+          2.18116 0.000448787 0 2.16109 0.00046691 0 2.15027 -9.35448e-07 0 2.1515 -0.000313439 0
+          2.15495 -0.000131277 0 2.15518 -2.79203e-05 0 2.15389 2.53645e-05 0 2.15243 2.25666e-05 0
+          2.52057 0.00617368 0 2.5554 0.00531584 0 2.56969 0.00392204 0 2.55829 0.00310417 0
+          2.53195 0.00262602 0 2.49793 0.00174611 0 2.45943 0.00105046 0 2.42443 0.000968236 0
+          2.39621 0.000898395 0 2.37341 0.000836654 0 2.35442 0.000218888 0 2.33973 -0.000228425 0
+          2.33042 -6.53359e-05 0 2.3245 3.32452e-05 0 2.3199 8.61041e-05 0 2.31679 5.33453e-05 0
+          2.52129 0.00506216 0 2.56089 0.00443614 0 2.59214 0.00352137 0 2.61312 0.00310863 0
+          2.6256 0.00281738 0 2.62815 0.00222645 0 2.61909 0.00169962 0 2.60431 0.00154789 0
+          2.58926 0.00140614 0 2.57519 0.00126665 0 2.56046 0.000662732 0 2.54428 0.000151795 0
+          2.53032 0.000190803 0 2.52022 0.000219873 0 2.51276 0.000231655 0 2.50815 0.000118542 0
+          2.52126 0.00348503 0 2.56117 0.00306486 0 2.59479 0.00247387 0 2.62329 0.00227611 0
+          2.64923 0.00214311 0 2.67127 0.00186911 0 2.68737 0.00164056 0 2.69894 0.00156339 0
+          2.70789 0.00146865 0 2.71459 0.00134947 0 2.71692 0.000934078 0 2.71333 0.000535553 0
+          2.70763 0.000484579 0 2.70291 0.000442433 0 2.69917 0.000401328 0 2.69684 0.000192124 0
+          2.52119 0.0012828 0 2.56104 0.00112934 0 2.59478 0.000913247 0 2.62384 0.000845638 0
+          2.65108 0.000803304 0 2.67635 0.000726364 0 2.69959 0.000680281 0 2.72193 0.000671616 0
+          2.74378 0.000650989 0 2.76466 0.000613382 0 2.7827 0.000478057 0 2.79612 0.000334349 0
+          2.80658 0.000298959 0 2.81593 0.000266699 0 2.82423 0.000235966 0 2.82997 0.000111265 0
+          2.52119 -0.0012828 0 2.56104 -0.00112934 0 2.59478 -0.000913247 0 2.62384 -0.000845638 0
+          2.65108 -0.000803304 0 2.67635 -0.000726364 0 2.69959 -0.000680281 0 2.72193 -0.000671616 0
+          2.74378 -0.000650989 0 2.76466 -0.000613382 0 2.7827 -0.000478057 0 2.79612 -0.000334349 0
+          2.80658 -0.000298959 0 2.81593 -0.000266699 0 2.82423 -0.000235966 0 2.82997 -0.000111265 0
+          2.52126 -0.00348503 0 2.56117 -0.00306486 0 2.59479 -0.00247387 0 2.62329 -0.00227611 0
+          2.64923 -0.00214311 0 2.67127 -0.00186911 0 2.68737 -0.00164056 0 2.69894 -0.00156339 0
+          2.70789 -0.00146865 0 2.71459 -0.00134947 0 2.71692 -0.000934078 0 2.71333 -0.000535553 0
+          2.70763 -0.000484579 0 2.70291 -0.000442433 0 2.69917 -0.000401328 0 2.69684 -0.000192124 0
+          2.52129 -0.00506216 0 2.56089 -0.00443614 0 2.59214 -0.00352137 0 2.61312 -0.00310863 0
+          2.6256 -0.00281738 0 2.62815 -0.00222645 0 2.61909 -0.00169962 0 2.60431 -0.00154789 0
+          2.58926 -0.00140614 0 2.57519 -0.00126665 0 2.56046 -0.000662732 0 2.54428 -0.000151795 0
+          2.53032 -0.000190803 0 2.52022 -0.000219873 0 2.51276 -0.000231655 0 2.50815 -0.000118542 0
+          2.52057 -0.00617368 0 2.5554 -0.00531584 0 2.56969 -0.00392204 0 2.55829 -0.00310417 0
+          2.53195 -0.00262602 0 2.49793 -0.00174611 0 2.45943 -0.00105046 0 2.42443 -0.000968236 0
+          2.39621 -0.000898395 0 2.37341 -0.000836654 0 2.35442 -0.000218888 0 2.33973 0.000228425 0
+          2.33042 6.53359e-05 0 2.3245 -3.32452e-05 0 2.3199 -8.61041e-05 0 2.31679 -5.33453e-05 0
+          2.51407 -0.00684927 0 2.52122 -0.00552131 0 2.48153 -0.00340618 0 2.41111 -0.00225882 0
+          2.34207 -0.00180202 0 2.2841 -0.000906265 0 2.23833 -0.000319485 0 2.20557 -0.000410678 0
+          2.18116 -0.000448787 0 2.16109 -0.00046691 0 2.15027 9.35448e-07 0 2.1515 0.000313439 0
+          2.15495 0.000131277 0 2.15518 2.79203e-05 0 2.15389 -2.53645e-05 0 2.15243 -2.25666e-05 0
+          2.48263 -0.00687549 0 2.40526 -0.00473685 0 2.27794 -0.00213196 0 2.16452 -0.00114677 0
+          2.08741 -0.000958049 0 2.04046 -0.00030835 0 2.01822 5.431e-05 0 2.00495 -0.000118866 0
+          1.99069 -0.000195461 0 1.97585 -0.000236882 0 1.97045 2.33479e-06 0 1.97996 0.000185671 0
+          1.98994 9.23723e-05 0 1.99318 2.67383e-05 0 1.99354 -7.8219e-06 0 1.99293 -1.0217e-05 0
+          2.38628 -0.00576564 0 2.14543 -0.00303948 0 1.94733 -0.000803422 0 1.84352 -0.000364102 0
+          1.78737 -0.000421259 0 1.76156 -0.000109178 0 1.76632 6.72728e-05 0 1.7723 -3.05429e-05 0
+          1.76675 -7.86469e-05 0 1.75635 -0.00010299 0 1.74955 -1.60399e-05 0 1.75309 6.85792e-05 0
+          1.75963 4.45018e-05 0 1.76266 1.38629e-05 0 1.76313 -2.77535e-06 0 1.76263 -4.39272e-06 0
+          2.07761 -0.00242566 0 1.47957 -0.00100872 0 1.28367 -0.000116336 0 1.25468 -5.01453e-05 0
+          1.22561 -0.000116787 0 1.1975 -4.46784e-05 0 1.19052 4.59293e-06 0 1.18953 -1.02362e-05 0
+          1.18391 -2.20593e-05 0 1.17538 -2.6933e-05 0 1.16978 -5.26205e-06 0 1.17188 1.73307e-05 0
+          1.17714 1.28878e-05 0 1.18011 4.15204e-06 0 1.18075 -5.06639e-07 0 1.18047 -1.08101e-06 0
+        </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 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 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 0.625 0 0 0 0.00358921 0 0.625 0.00358921 0
+          1.25 0 0 1.25 0.00358921 0 1.875 0 0 1.875 0.00358921 0
+          2.5 0 0 2.5 0.00358921 0 3.125 0 0 3.125 0.00358921 0
+          3.75 0 0 3.75 0.00358921 0 4.375 0 0 4.375 0.00358921 0
+          5 0 0 5 0.00358921 0 5.625 0 0 5.625 0.00358921 0
+          6.25 0 0 6.25 0.00358921 0 6.875 0 0 6.875 0.00358921 0
+          7.5 0 0 7.5 0.00358921 0 8.125 0 0 8.125 0.00358921 0
+          8.75 0 0 8.75 0.00358921 0 9.375 0 0 9.375 0.00358921 0
+          10 0 0 10 0.00358921 0 0 0.00861411 0 0.625 0.00861411 0
+          1.25 0.00861411 0 1.875 0.00861411 0 2.5 0.00861411 0 3.125 0.00861411 0
+          3.75 0.00861411 0 4.375 0.00861411 0 5 0.00861411 0 5.625 0.00861411 0
+          6.25 0.00861411 0 6.875 0.00861411 0 7.5 0.00861411 0 8.125 0.00861411 0
+          8.75 0.00861411 0 9.375 0.00861411 0 10 0.00861411 0 0 0.015649 0
+          0.625 0.015649 0 1.25 0.015649 0 1.875 0.015649 0 2.5 0.015649 0
+          3.125 0.015649 0 3.75 0.015649 0 4.375 0.015649 0 5 0.015649 0
+          5.625 0.015649 0 6.25 0.015649 0 6.875 0.015649 0 7.5 0.015649 0
+          8.125 0.015649 0 8.75 0.015649 0 9.375 0.015649 0 10 0.015649 0
+          0 0.0254978 0 0.625 0.0254978 0 1.25 0.0254978 0 1.875 0.0254978 0
+          2.5 0.0254978 0 3.125 0.0254978 0 3.75 0.0254978 0 4.375 0.0254978 0
+          5 0.0254978 0 5.625 0.0254978 0 6.25 0.0254978 0 6.875 0.0254978 0
+          7.5 0.0254978 0 8.125 0.0254978 0 8.75 0.0254978 0 9.375 0.0254978 0
+          10 0.0254978 0 0 0.0392861 0 0.625 0.0392861 0 1.25 0.0392861 0
+          1.875 0.0392861 0 2.5 0.0392861 0 3.125 0.0392861 0 3.75 0.0392861 0
+          4.375 0.0392861 0 5 0.0392861 0 5.625 0.0392861 0 6.25 0.0392861 0
+          6.875 0.0392861 0 7.5 0.0392861 0 8.125 0.0392861 0 8.75 0.0392861 0
+          9.375 0.0392861 0 10 0.0392861 0 0 0.0585897 0 0.625 0.0585897 0
+          1.25 0.0585897 0 1.875 0.0585897 0 2.5 0.0585897 0 3.125 0.0585897 0
+          3.75 0.0585897 0 4.375 0.0585897 0 5 0.0585897 0 5.625 0.0585897 0
+          6.25 0.0585897 0 6.875 0.0585897 0 7.5 0.0585897 0 8.125 0.0585897 0
+          8.75 0.0585897 0 9.375 0.0585897 0 10 0.0585897 0 0 0.0856148 0
+          0.625 0.0856148 0 1.25 0.0856148 0 1.875 0.0856148 0 2.5 0.0856148 0
+          3.125 0.0856148 0 3.75 0.0856148 0 4.375 0.0856148 0 5 0.0856148 0
+          5.625 0.0856148 0 6.25 0.0856148 0 6.875 0.0856148 0 7.5 0.0856148 0
+          8.125 0.0856148 0 8.75 0.0856148 0 9.375 0.0856148 0 10 0.0856148 0
+          0 0.12345 0 0.625 0.12345 0 1.25 0.12345 0 1.875 0.12345 0
+          2.5 0.12345 0 3.125 0.12345 0 3.75 0.12345 0 4.375 0.12345 0
+          5 0.12345 0 5.625 0.12345 0 6.25 0.12345 0 6.875 0.12345 0
+          7.5 0.12345 0 8.125 0.12345 0 8.75 0.12345 0 9.375 0.12345 0
+          10 0.12345 0 0 0.161285 0 0.625 0.161285 0 1.25 0.161285 0
+          1.875 0.161285 0 2.5 0.161285 0 3.125 0.161285 0 3.75 0.161285 0
+          4.375 0.161285 0 5 0.161285 0 5.625 0.161285 0 6.25 0.161285 0
+          6.875 0.161285 0 7.5 0.161285 0 8.125 0.161285 0 8.75 0.161285 0
+          9.375 0.161285 0 10 0.161285 0 0 0.18831 0 0.625 0.18831 0
+          1.25 0.18831 0 1.875 0.18831 0 2.5 0.18831 0 3.125 0.18831 0
+          3.75 0.18831 0 4.375 0.18831 0 5 0.18831 0 5.625 0.18831 0
+          6.25 0.18831 0 6.875 0.18831 0 7.5 0.18831 0 8.125 0.18831 0
+          8.75 0.18831 0 9.375 0.18831 0 10 0.18831 0 0 0.207614 0
+          0.625 0.207614 0 1.25 0.207614 0 1.875 0.207614 0 2.5 0.207614 0
+          3.125 0.207614 0 3.75 0.207614 0 4.375 0.207614 0 5 0.207614 0
+          5.625 0.207614 0 6.25 0.207614 0 6.875 0.207614 0 7.5 0.207614 0
+          8.125 0.207614 0 8.75 0.207614 0 9.375 0.207614 0 10 0.207614 0
+          0 0.221402 0 0.625 0.221402 0 1.25 0.221402 0 1.875 0.221402 0
+          2.5 0.221402 0 3.125 0.221402 0 3.75 0.221402 0 4.375 0.221402 0
+          5 0.221402 0 5.625 0.221402 0 6.25 0.221402 0 6.875 0.221402 0
+          7.5 0.221402 0 8.125 0.221402 0 8.75 0.221402 0 9.375 0.221402 0
+          10 0.221402 0 0 0.231251 0 0.625 0.231251 0 1.25 0.231251 0
+          1.875 0.231251 0 2.5 0.231251 0 3.125 0.231251 0 3.75 0.231251 0
+          4.375 0.231251 0 5 0.231251 0 5.625 0.231251 0 6.25 0.231251 0
+          6.875 0.231251 0 7.5 0.231251 0 8.125 0.231251 0 8.75 0.231251 0
+          9.375 0.231251 0 10 0.231251 0 0 0.238286 0 0.625 0.238286 0
+          1.25 0.238286 0 1.875 0.238286 0 2.5 0.238286 0 3.125 0.238286 0
+          3.75 0.238286 0 4.375 0.238286 0 5 0.238286 0 5.625 0.238286 0
+          6.25 0.238286 0 6.875 0.238286 0 7.5 0.238286 0 8.125 0.238286 0
+          8.75 0.238286 0 9.375 0.238286 0 10 0.238286 0 0 0.243311 0
+          0.625 0.243311 0 1.25 0.243311 0 1.875 0.243311 0 2.5 0.243311 0
+          3.125 0.243311 0 3.75 0.243311 0 4.375 0.243311 0 5 0.243311 0
+          5.625 0.243311 0 6.25 0.243311 0 6.875 0.243311 0 7.5 0.243311 0
+          8.125 0.243311 0 8.75 0.243311 0 9.375 0.243311 0 10 0.243311 0
+          0 0.2469 0 0.625 0.2469 0 1.25 0.2469 0 1.875 0.2469 0
+          2.5 0.2469 0 3.125 0.2469 0 3.75 0.2469 0 4.375 0.2469 0
+          5 0.2469 0 5.625 0.2469 0 6.25 0.2469 0 6.875 0.2469 0
+          7.5 0.2469 0 8.125 0.2469 0 8.75 0.2469 0 9.375 0.2469 0
+          10 0.2469 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 10 12 13 11
+          12 14 15 13 14 16 17 15 16 18 19 17
+          18 20 21 19 20 22 23 21 22 24 25 23
+          24 26 27 25 26 28 29 27 28 30 31 29
+          30 32 33 31 2 3 35 34 3 5 36 35
+          5 7 37 36 7 9 38 37 9 11 39 38
+          11 13 40 39 13 15 41 40 15 17 42 41
+          17 19 43 42 19 21 44 43 21 23 45 44
+          23 25 46 45 25 27 47 46 27 29 48 47
+          29 31 49 48 31 33 50 49 34 35 52 51
+          35 36 53 52 36 37 54 53 37 38 55 54
+          38 39 56 55 39 40 57 56 40 41 58 57
+          41 42 59 58 42 43 60 59 43 44 61 60
+          44 45 62 61 45 46 63 62 46 47 64 63
+          47 48 65 64 48 49 66 65 49 50 67 66
+          51 52 69 68 52 53 70 69 53 54 71 70
+          54 55 72 71 55 56 73 72 56 57 74 73
+          57 58 75 74 58 59 76 75 59 60 77 76
+          60 61 78 77 61 62 79 78 62 63 80 79
+          63 64 81 80 64 65 82 81 65 66 83 82
+          66 67 84 83 68 69 86 85 69 70 87 86
+          70 71 88 87 71 72 89 88 72 73 90 89
+          73 74 91 90 74 75 92 91 75 76 93 92
+          76 77 94 93 77 78 95 94 78 79 96 95
+          79 80 97 96 80 81 98 97 81 82 99 98
+          82 83 100 99 83 84 101 100 85 86 103 102
+          86 87 104 103 87 88 105 104 88 89 106 105
+          89 90 107 106 90 91 108 107 91 92 109 108
+          92 93 110 109 93 94 111 110 94 95 112 111
+          95 96 113 112 96 97 114 113 97 98 115 114
+          98 99 116 115 99 100 117 116 100 101 118 117
+          102 103 120 119 103 104 121 120 104 105 122 121
+          105 106 123 122 106 107 124 123 107 108 125 124
+          108 109 126 125 109 110 127 126 110 111 128 127
+          111 112 129 128 112 113 130 129 113 114 131 130
+          114 115 132 131 115 116 133 132 116 117 134 133
+          117 118 135 134 119 120 137 136 120 121 138 137
+          121 122 139 138 122 123 140 139 123 124 141 140
+          124 125 142 141 125 126 143 142 126 127 144 143
+          127 128 145 144 128 129 146 145 129 130 147 146
+          130 131 148 147 131 132 149 148 132 133 150 149
+          133 134 151 150 134 135 152 151 136 137 154 153
+          137 138 155 154 138 139 156 155 139 140 157 156
+          140 141 158 157 141 142 159 158 142 143 160 159
+          143 144 161 160 144 145 162 161 145 146 163 162
+          146 147 164 163 147 148 165 164 148 149 166 165
+          149 150 167 166 150 151 168 167 151 152 169 168
+          153 154 171 170 154 155 172 171 155 156 173 172
+          156 157 174 173 157 158 175 174 158 159 176 175
+          159 160 177 176 160 161 178 177 161 162 179 178
+          162 163 180 179 163 164 181 180 164 165 182 181
+          165 166 183 182 166 167 184 183 167 168 185 184
+          168 169 186 185 170 171 188 187 171 172 189 188
+          172 173 190 189 173 174 191 190 174 175 192 191
+          175 176 193 192 176 177 194 193 177 178 195 194
+          178 179 196 195 179 180 197 196 180 181 198 197
+          181 182 199 198 182 183 200 199 183 184 201 200
+          184 185 202 201 185 186 203 202 187 188 205 204
+          188 189 206 205 189 190 207 206 190 191 208 207
+          191 192 209 208 192 193 210 209 193 194 211 210
+          194 195 212 211 195 196 213 212 196 197 214 213
+          197 198 215 214 198 199 216 215 199 200 217 216
+          200 201 218 217 201 202 219 218 202 203 220 219
+          204 205 222 221 205 206 223 222 206 207 224 223
+          207 208 225 224 208 209 226 225 209 210 227 226
+          210 211 228 227 211 212 229 228 212 213 230 229
+          213 214 231 230 214 215 232 231 215 216 233 232
+          216 217 234 233 217 218 235 234 218 219 236 235
+          219 220 237 236 221 222 239 238 222 223 240 239
+          223 224 241 240 224 225 242 241 225 226 243 242
+          226 227 244 243 227 228 245 244 228 229 246 245
+          229 230 247 246 230 231 248 247 231 232 249 248
+          232 233 250 249 233 234 251 250 234 235 252 251
+          235 236 253 252 236 237 254 253 238 239 256 255
+          239 240 257 256 240 241 258 257 241 242 259 258
+          242 243 260 259 243 244 261 260 244 245 262 261
+          245 246 263 262 246 247 264 263 247 248 265 264
+          248 249 266 265 249 250 267 266 250 251 268 267
+          251 252 269 268 252 253 270 269 253 254 271 270
+          255 256 273 272 256 257 274 273 257 258 275 274
+          258 259 276 275 259 260 277 276 260 261 278 277
+          261 262 279 278 262 263 280 279 263 264 281 280
+          264 265 282 281 265 266 283 282 266 267 284 283
+          267 268 285 284 268 269 286 285 269 270 287 286
+          270 271 288 287
+        </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 404 408 412 416 420 424 428 432
+          436 440 444 448 452 456 460 464 468 472 476 480
+          484 488 492 496 500 504 508 512 516 520 524 528
+          532 536 540 544 548 552 556 560 564 568 572 576
+          580 584 588 592 596 600 604 608 612 616 620 624
+          628 632 636 640 644 648 652 656 660 664 668 672
+          676 680 684 688 692 696 700 704 708 712 716 720
+          724 728 732 736 740 744 748 752 756 760 764 768
+          772 776 780 784 788 792 796 800 804 808 812 816
+          820 824 828 832 836 840 844 848 852 856 860 864
+          868 872 876 880 884 888 892 896 900 904 908 912
+          916 920 924 928 932 936 940 944 948 952 956 960
+          964 968 972 976 980 984 988 992 996 1000 1004 1008
+          1012 1016 1020 1024
+        </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 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 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>