diff --git a/appl/coupling-ff-pm/CMakeLists.txt b/appl/coupling-ff-pm/CMakeLists.txt
index a2de49f3bba1495c15a2ed0b206c8c449905405b..d3e96b83cfd100be34c5e8da7806d8089e6de69c 100644
--- a/appl/coupling-ff-pm/CMakeLists.txt
+++ b/appl/coupling-ff-pm/CMakeLists.txt
@@ -4,3 +4,4 @@ add_custom_target(test_exercises)
 
 add_subdirectory(monolithic)
 add_subdirectory(iterative)
+add_subdirectory(iterative-reversed)
diff --git a/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh
index 5c7fd43b24215c4d550d2cb1277b4bfbeeccef00..97390739a2de1c352e551fae7f35c70dca01f6cf 100644
--- a/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh
+++ b/appl/coupling-ff-pm/iterative-reversed/ffproblem-reversed.hh
@@ -206,7 +206,7 @@ public:
           values.setNeumann(Indices::momentumYBalanceIdx);
           //values.setCouplingNeumann(Indices::conti0EqIdx);
           //values.setCouplingNeumann(Indices::momentumYBalanceIdx);
-          //values.setBJS(Indices::momentumXBalanceIdx);
+          values.setBJS(Indices::momentumXBalanceIdx);
         }
 #endif
 
@@ -258,22 +258,23 @@ public:
           const auto velocity = couplingInterface_.getScalarQuantityOnFace( velocityId_, faceId );
           const auto& volVars = elemVolVars[scvf.insideScvIdx()];
 
-//          const Scalar ccPressure = volVars.pressure();
-//          const Scalar mobility = volVars.mobility();
+          const Scalar ccPressure = volVars.pressure();
+          const Scalar mobility = volVars.mobility();
           const Scalar density = volVars.density();
-//          const auto K = volVars.permeability();
+          const auto K = volVars.permeability();
 
-//          // v = -kr/mu*K * (gradP + rho*g) = -mobility*K * (gradP + rho*g)
-//          const auto alpha = Dumux::vtmv( scvf.unitOuterNormal(), K, this->gravity() );
+          // v = -kr/mu*K * (gradP + rho*g) = -mobility*K * (gradP + rho*g)
+          const auto alpha = Dumux::vtmv( scvf.unitOuterNormal(), K, this->gravity() );
 
-//          auto distanceVector = scvf.center() - element.geometry().center();
-//          distanceVector /= distanceVector.two_norm2();
-//          const Scalar ti = Dumux::vtmv(distanceVector, K, scvf.unitOuterNormal());
+          auto distanceVector = scvf.center() - element.geometry().center();
+          distanceVector /= distanceVector.two_norm2();
+          const Scalar ti = Dumux::vtmv(distanceVector, K, scvf.unitOuterNormal());
 
-//          const Scalar pressure = (1/mobility * (scvf.unitOuterNormal() * velocity) + density * alpha)/ti
-//              + ccPressure;
+          const Scalar pressure = (1/mobility * (scvf.unitOuterNormal() * velocity) + density * alpha)/ti
+              + ccPressure;
 
           values[Indices::conti0EqIdx] = velocity * density;
+          values[Indices::momentumYBalanceIdx] = scvf.directionSign() * ( pressure - initialAtPos(scvf.center())[Indices::pressureIdx]) ;
 //          values[Indices::momentumYBalanceIdx] = scvf.directionSign() * ( pressure - initialAtPos(scvf.center())[Indices::pressureIdx]) ;
         }
 #endif
diff --git a/appl/coupling-ff-pm/iterative-reversed/pmproblem-reversed.hh b/appl/coupling-ff-pm/iterative-reversed/pmproblem-reversed.hh
index 1dd752fc022e2847dbe0f4240002db13938a5d4d..dfca6d5cf12703aeb9fcd6de21d010a8d1fbb65d 100644
--- a/appl/coupling-ff-pm/iterative-reversed/pmproblem-reversed.hh
+++ b/appl/coupling-ff-pm/iterative-reversed/pmproblem-reversed.hh
@@ -170,6 +170,17 @@ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
         // set the coupling boundary condition at the interface
         if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
             values.setAllCouplingNeumann();
+#else
+    // // TODO do preCICE stuff in analogy to heat transfer
+        assert( dataIdsWereSet_ );
+        const auto faceId = scvf.index();
+        if ( couplingInterface_.isCoupledEntity(faceId) )
+        {
+          //TODO What do I want to do here?
+          values.setCouplingNeumann(Indices::conti0EqIdx);
+          values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+          values.setBJS(Indices::momentumXBalanceIdx);
+        }
 #endif
         return values;
     }
@@ -215,13 +226,15 @@ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
         if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
             values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
 #else
-        assert( dataIdsWereSet_ );
+        zassert( dataIdsWereSet_ );
         const auto faceId = scvf.index();
-        if ( couplingInterface_.isCoupledEntity(faceId) )
+        if( couplingInterface_.isCoupledEntity( faceId ) )
         {
-//          const Scalar density = 1000.;
-          //values[Indices::conti0EqIdx] = density * couplingInterface_.getScalarQuantityOnFace( velocityId_, faceId );
-          const Scalar facePressure = couplingInterface_.getScalarQuantityOnFace( pressureId_, faceId );
+          const Scalar densityOnFace = 1000; // TODO how to handle compressible fluids?
+
+          const Scalar pressure = couplingInterface_.getScalarQuantityOnFace( pressureId_, faceId );
+
+          const auto& volVars = elemVolVars[scvf.insideScvIdx()];
 
           const Scalar ccPressure = volVars.pressure();
           const Scalar mobility = volVars.mobility();
@@ -229,9 +242,18 @@ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
           const auto K = volVars.permeability();
 
           // v = -kr/mu*K * (gradP + rho*g) = -mobility*K * (gradP + rho*g)
-          const auto alpha = Dumux::vtmv( scvf.unitOuterNormal(), K, problem.gravity() );
+          const auto alpha = Dumux::vtmv( scvf.unitOuterNormal(), K, this->gravity() );
+
+          auto distanceVector = scvf.center() - element.geometry().center();
+          distanceVector /= distanceVector.two_norm2();
+          const Scalar ti = Dumux::vtmv(distanceVector, K, scvf.unitOuterNormal());
+
+//          const Scalar pressure = (1/mobility * (scvf.unitOuterNormal() * velocity) + density * alpha)/ti + ccPressure;
+//          mobility * ( (pressure - ccPressure) * ti - + density * alpha ) / scvf.unitOuterNormal()= velocity  ;
+          const Scalar velocity = mobility * ( (pressure - ccPressure) * ti - + density * alpha ) / scvf.unitOuterNormal();
 
-          values[Indices::conti0EqIdx] = density * ;
+          values[Indices::conti0EqIdx] = velocity * scvf.directionSign() * densityOnFace;
+          values[Indices::momentumYBalanceIdx] = scvf.directionSign() * (couplingInterface_.getScalarQuantityOnFace( pressureId_, faceId ) - initialAtPos(scvf.center())[Indices::pressureIdx]) ;
         }
 #endif
         return values;
diff --git a/appl/coupling-ff-pm/iterative/main_ff.cc b/appl/coupling-ff-pm/iterative/main_ff.cc
index 161e15a4512cf1c51690978d4d2c911649d621cf..8b32d2dcd87dc414cea3ed70f27e0573ab053a43 100644
--- a/appl/coupling-ff-pm/iterative/main_ff.cc
+++ b/appl/coupling-ff-pm/iterative/main_ff.cc
@@ -62,6 +62,62 @@ auto velocityAtInterface(const ElementFaceVariables& elemFaceVars,
     return velocity;
  }
 
+template<class FluxVariables,
+         class Problem,
+         class Element,
+         class SubControlVolumeFace,
+         class FVElementGeometry,
+         class ElementVolumeVariables,
+         class ElementFaceVariables,
+         class ElementFluxVariablesCache>
+auto pressureAtInterface(const Problem& problem,
+                         const Element& element,
+                         const SubControlVolumeFace& scvf,
+                         const FVElementGeometry& fvGeometry,
+                         const ElementVolumeVariables& elemVolVars,
+                         const ElementFaceVariables& elemFaceVars,
+                         const ElementFluxVariablesCache& elemFluxVarsCache)
+{
+  FluxVariables fluxVars;
+  return fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache()) /scvf.area();
+}
+
+template<class FluxVariables, class Problem, class GridVariables, class SolutionVector>
+void setInterfacePressures(const Problem& problem,
+                            const GridVariables& gridVars,
+                            const SolutionVector& sol)
+{
+  const auto& fvGridGeometry = problem.fvGridGeometry();
+  auto fvGeometry = localView(fvGridGeometry);
+  auto elemVolVars = localView(gridVars.curGridVolVars());
+  auto elemFaceVars = localView(gridVars.curGridFaceVars());
+  auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
+
+  auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
+  const auto pressureId = couplingInterface.getIdFromName( "Pressure" );
+
+  for (const auto& element : elements(fvGridGeometry.gridView()))
+  {
+    fvGeometry.bind(element);
+    elemVolVars.bind(element, fvGeometry, sol);
+    elemFaceVars.bind(element, fvGeometry, sol);
+    elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
+
+    for (const auto& scvf : scvfs(fvGeometry))
+    {
+
+      if ( couplingInterface.isCoupledEntity( scvf.index() ) )
+      {
+        //TODO: What to do here?
+        const auto p = pressureAtInterface<FluxVariables>(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
+        couplingInterface.writeScalarQuantityOnFace( pressureId, scvf.index(), p );
+      }
+    }
+  }
+
+}
+
+
 template<class Problem, class GridVariables, class SolutionVector>
 void setInterfaceVelocities(const Problem& problem,
                             const GridVariables& gridVars,
@@ -207,12 +263,23 @@ int main(int argc, char** argv) try
     freeFlowVtkWriter.addField(freeFlowProblem->getAnalyticalVelocityX(), "analyticalV_x");
     freeFlowVtkWriter.write(0.0);
 
+     using FluxVariables = GetPropType<FreeFlowTypeTag, Properties::FluxVariables>;
+
     if ( couplingInterface.hasToWriteInitialData() )
     {
       //TODO
 //      couplingInterface.writeQuantityVector( pressureId );
-      setInterfaceVelocities( *freeFlowProblem, *freeFlowGridVariables, sol );
-      couplingInterface.writeScalarQuantityToOtherSolver( velocityId );
+
+      setInterfacePressures<FluxVariables>( *freeFlowProblem, *freeFlowGridVariables, sol );
+      //For testing
+      {
+        std::cout << "Pressures to be sent to pm" << std::endl;
+        const auto p = couplingInterface.getQuantityVector( pressureId );
+        for (size_t i = 0; i < p.size(); ++i) {
+          std::cout << "p[" << i << "]=" <<p[i] << std::endl;
+        }
+      }
+      couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
       couplingInterface.announceInitialDataWritten();
     }
     couplingInterface.initializeData();
@@ -233,6 +300,8 @@ int main(int argc, char** argv) try
     auto dt = preciceDt;
     auto sol_checkpoint = sol;
 
+    double vtkTime = 1.0;
+
     while ( couplingInterface.isCouplingOngoing() )
     {
         if ( couplingInterface.hasToWriteIterationCheckpoint() )
@@ -243,27 +312,43 @@ int main(int argc, char** argv) try
         }
 
         // TODO
-        couplingInterface.readScalarQuantityFromOtherSolver( pressureId );
-//        // For testing
-//        {
-//          const auto v = couplingInterface.getQuantityVector( pressureId );
-//          const double sum = std::accumulate( v.begin(), v.end(), 0. );
-//          std::cout << "Sum of pressures over boundary to pm: \n" << sum << std::endl;
-//        }
+        couplingInterface.readScalarQuantityFromOtherSolver( velocityId );
+        // For testing
+        {
+          const auto v = couplingInterface.getQuantityVector( velocityId );
+          const double sum = std::accumulate( v.begin(), v.end(), 0. );
+          std::cout << "Sum of pressures over boundary to pm: \n" << sum << std::endl;
+        }
 
         // solve the non-linear system
         nonLinearSolver.solve(sol);
 
         // TODO
-        setInterfaceVelocities( *freeFlowProblem, *freeFlowGridVariables, sol );
-        couplingInterface.writeScalarQuantityToOtherSolver( velocityId );
+        setInterfacePressures<FluxVariables>( *freeFlowProblem, *freeFlowGridVariables, sol );
+        // For testing
+        {
+          const auto p = couplingInterface.getQuantityVector( pressureId );
+          const double sum = std::accumulate( p.begin(), p.end(), 0. );
+          std::cout << "Pressures to be sent to pm" << std::endl;
+          for (size_t i = 0; i < p.size(); ++i) {
+            std::cout << "p[" << i << "]=" << p[i] << std::endl;
+          }
+          std::cout << "Sum of pressures over boundary to pm: \n" << sum << std::endl;
+        }
+        couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
+
 
+        //Read checkpoint
+        freeFlowVtkWriter.write(vtkTime);
+        vtkTime += 1.;
         const double preciceDt = couplingInterface.advance( dt );
         dt = std::min( preciceDt, dt );
 
         if ( couplingInterface.hasToReadIterationCheckpoint() )
         {
-            //Read checkpoint
+//            //Read checkpoint
+//            freeFlowVtkWriter.write(vtkTime);
+//            vtkTime += 1.;
             sol = sol_checkpoint;
             freeFlowGridVariables->update(sol);
             freeFlowGridVariables->advanceTimeStep();
@@ -273,7 +358,7 @@ int main(int argc, char** argv) try
         else // coupling successful
         {
           // write vtk output
-          freeFlowVtkWriter.write(1.0);
+          freeFlowVtkWriter.write(vtkTime);
         }
 
     }
diff --git a/appl/coupling-ff-pm/iterative/main_pm.cc b/appl/coupling-ff-pm/iterative/main_pm.cc
index 2094dd455debf9940e52157ea1f3d9639094f1a1..232201dffd5d4749b7232f148b17a82aaf7f3ee4 100644
--- a/appl/coupling-ff-pm/iterative/main_pm.cc
+++ b/appl/coupling-ff-pm/iterative/main_pm.cc
@@ -69,14 +69,15 @@
 
      const Scalar boundaryFlux = problem.neumann(element, fvGeometry, elemVolVars, scvf);
 
+     const auto K = volVars.permeability();
+     const Scalar ccPressure = volVars.pressure();
+     const Scalar mobility = volVars.mobility();
+     const Scalar density = volVars.density();
+
      // v = -K/mu * (gradP + rho*g)
      auto velocity = scvf.unitOuterNormal();
      velocity *= boundaryFlux; // TODO check sign
-
-    const Scalar ccPressure = volVars.pressure();
-    const Scalar mobility = volVars.mobility();
-    const Scalar density = volVars.density();
-    const auto K = volVars.permeability();
+     velocity /= density;
 
     // v = -kr/mu*K * (gradP + rho*g) = -mobility*K * (gradP + rho*g)
     const auto alpha = Dumux::vtmv( scvf.unitOuterNormal(), K, problem.gravity() );
@@ -120,6 +121,60 @@
 
  }
 
+ /*!
+  * \brief Returns the pressure at the interface using Darcy's law for reconstruction
+  */
+ template<class FluxVariables, class Problem, class Element, class FVElementGeometry,
+          class ElementVolumeVariables, class SubControlVolumeFace,
+          class ElementFluxVariablesCache>
+ auto velocityAtInterface(const Problem& problem,
+                            const Element& element,
+                            const FVElementGeometry& fvGeometry,
+                            const ElementVolumeVariables& elemVolVars,
+                            const SubControlVolumeFace& scvf,
+                            const ElementFluxVariablesCache& elemFluxVarsCache)
+ {
+   const int phaseIdx = 0;
+   FluxVariables fluxVars;
+   fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+   auto upwindTerm = [phaseIdx](const auto& volVars) { return volVars.mobility(phaseIdx); };
+   const auto scalarVelocity = fluxVars.advectiveFlux(phaseIdx, upwindTerm)/scvf.area();
+   return scalarVelocity;
+ }
+
+ template<class FluxVariables, class Problem, class GridVariables, class SolutionVector>
+ void setInterfaceVelocities(const Problem& problem,
+                            const GridVariables& gridVars,
+                            const SolutionVector& sol)
+ {
+   const auto& fvGridGeometry = problem.fvGridGeometry();
+   auto fvGeometry = localView(fvGridGeometry);
+   auto elemVolVars = localView(gridVars.curGridVolVars());
+   auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
+
+   auto& couplingInterface = precice_adapter::PreciceAdapter::getInstance();
+   const auto velocityId = couplingInterface.getIdFromName( "Velocity" );
+
+   for (const auto& element : elements(fvGridGeometry.gridView()))
+   {
+     fvGeometry.bind(element);
+     elemVolVars.bind(element, fvGeometry, sol);
+     elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
+
+     for (const auto& scvf : scvfs(fvGeometry))
+     {
+
+       if ( couplingInterface.isCoupledEntity( scvf.index() ) )
+       {
+         //TODO: What to do here?
+         const double v = velocityAtInterface<FluxVariables>(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+         couplingInterface.writeScalarQuantityOnFace( velocityId, scvf.index(), v );
+       }
+     }
+   }
+
+ }
+
 int main(int argc, char** argv) try
 {
     using namespace Dumux;
@@ -228,12 +283,20 @@ int main(int argc, char** argv) try
     GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
     darcyVtkWriter.write(0.0);
 
-
+    using FluxVariables = GetPropType<DarcyTypeTag, Properties::FluxVariables>;
     if ( couplingInterface.hasToWriteInitialData() )
     {
       //TODO
       //couplingInterface.writeQuantityVector(velocityId);
-      setInterfacePressures( *darcyProblem, *darcyGridVariables, sol );
+      setInterfaceVelocities<FluxVariables>( *darcyProblem, *darcyGridVariables, sol );
+      // For testing
+      {
+        const auto v = couplingInterface.getQuantityVector( velocityId );
+        std::cout << "velocities to be sent to ff" << std::endl;
+        for (size_t i = 0; i < v.size(); ++i) {
+          std::cout << "v[" << i << "]=" << v[i] << std::endl;
+        }
+      }
       couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
       couplingInterface.announceInitialDataWritten();
     }
@@ -254,6 +317,8 @@ int main(int argc, char** argv) try
     auto dt = preciceDt;
     auto sol_checkpoint = sol;
 
+    double vtkTime = 1.0;
+
     while ( couplingInterface.isCouplingOngoing() )
     {
         if ( couplingInterface.hasToWriteIterationCheckpoint() )
@@ -264,18 +329,32 @@ int main(int argc, char** argv) try
         }
 
         // TODO
-        couplingInterface.readScalarQuantityFromOtherSolver( velocityId );
-//        // For testing
-//        {
-//          const auto v = couplingInterface.getQuantityVector( velocityId );
-//          const double sum = std::accumulate( v.begin(), v.end(), 0. );
-//          std::cout << "Sum of fluxes over boundary to pm: \n" << sum << std::endl;
-//        }
+        couplingInterface.readScalarQuantityFromOtherSolver( pressureId );
+        // For testing
+        {
+          const auto p = couplingInterface.getQuantityVector( pressureId );
+          const double sum = std::accumulate( p.begin(), p.end(), 0. );
+          std::cout << "Sum of pressures over boundary to ff: \n" << sum << std::endl;
+          std::cout << "Pressure received from ff" << std::endl;
+          for (size_t i = 0; i < p.size(); ++i) {
+            std::cout << "p[" << i << "]=" << p[i] << std::endl;
+          }
+        }
 
         // solve the non-linear system
         nonLinearSolver.solve(sol);
-        setInterfacePressures( *darcyProblem, *darcyGridVariables, sol );
-        couplingInterface.writeScalarQuantityToOtherSolver( pressureId );
+        setInterfaceVelocities<FluxVariables>( *darcyProblem, *darcyGridVariables, sol );
+        // For testing
+        {
+          const auto v = couplingInterface.getQuantityVector( velocityId );
+          const double sum = std::accumulate( v.begin(), v.end(), 0. );
+          std::cout << "Velocities to be sent to ff" << std::endl;
+          for (size_t i = 0; i < v.size(); ++i) {
+            std::cout << "v[" << i << "]=" << v[i] << std::endl;
+          }
+          std::cout << "Sum of velocities over boundary to ff: \n" << sum << std::endl;
+        }
+        couplingInterface.writeScalarQuantityToOtherSolver( velocityId );
 
         const double preciceDt = couplingInterface.advance( dt );
         dt = std::min( preciceDt, dt );
@@ -283,6 +362,8 @@ int main(int argc, char** argv) try
         if ( couplingInterface.hasToReadIterationCheckpoint() )
         {
             //Read checkpoint
+            darcyVtkWriter.write(vtkTime);
+            vtkTime += 1.;
             sol = sol_checkpoint;
             darcyGridVariables->update(sol);
             darcyGridVariables->advanceTimeStep();
@@ -292,7 +373,7 @@ int main(int argc, char** argv) try
         else // coupling successful
         {
           // write vtk output
-          darcyVtkWriter.write(1.0);
+          darcyVtkWriter.write(vtkTime);
         }
 
     }
diff --git a/appl/coupling-ff-pm/iterative/params.input b/appl/coupling-ff-pm/iterative/params.input
index ab27df78f79ec9e5e3de6ecdc9e402ce778f7a8d..7f2845d53b4ac0368fa7240592d6bce4f22e8f17 100644
--- a/appl/coupling-ff-pm/iterative/params.input
+++ b/appl/coupling-ff-pm/iterative/params.input
@@ -18,8 +18,8 @@ Grading1 = 1
 
 [FreeFlow.Problem]
 Name = stokes
-PressureDifference = 1e-9
-EnableInertiaTerms = false
+PressureDifference = 1e-5
+EnableInertiaTerms = true
 
 [Darcy.Problem]
 Name = darcy
@@ -30,9 +30,13 @@ Permeability = 1e-6 # m^2
 Porosity = 0.4
 AlphaBeaversJoseph = 1.0
 
+[Newton]
+MaxSteps = 20
+
 [Problem]
 Name = ex_ff-pm-interface
 EnableGravity = false
 
 [Vtk]
 AddVelocity = 1
+WriteFaceData = 1
diff --git a/appl/coupling-ff-pm/monolithic/CMakeLists.txt b/appl/coupling-ff-pm/monolithic/CMakeLists.txt
index 4d6a867cb8db5d2ec8339c0bb18a7ccbc980a70e..0c64961e1a7bf0f66fb47e5daf6b5e20d3012117 100644
--- a/appl/coupling-ff-pm/monolithic/CMakeLists.txt
+++ b/appl/coupling-ff-pm/monolithic/CMakeLists.txt
@@ -1,6 +1,8 @@
 # executables for the monolithic case
 dune_add_test(NAME ff-pm_monolithic
               SOURCES main.cc)
+              
+target_compile_definitions(ff-pm_monolithic PUBLIC "ENABLEMONOLITHIC=1")
 
 # add a symlink for each input file
 add_input_file_links()
diff --git a/appl/coupling-ff-pm/monolithic/ffproblem.hh b/appl/coupling-ff-pm/monolithic/ffproblem.hh
index f1b4c4a86954e053b94bf40a213a9d0a308eed5d..668305d1fc8ac325a80f2bad43e17d986ad99842 100644
--- a/appl/coupling-ff-pm/monolithic/ffproblem.hh
+++ b/appl/coupling-ff-pm/monolithic/ffproblem.hh
@@ -202,8 +202,12 @@ public:
         if ( couplingInterface_.isCoupledEntity(faceId) )
         {
           //TODO What do I want to do here?
-          values.setCouplingNeumann(Indices::conti0EqIdx);
-          values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+        //  values.setCouplingNeumann(Indices::conti0EqIdx);
+        //  values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+          values.setDirichlet(Indices::velocityYIdx);
+
+//          values.setNeumann(Indices::conti0EqIdx);
+//          values.setNeumann(Indices::momentumYBalanceIdx);
           values.setBJS(Indices::momentumXBalanceIdx);
         }
 #endif
@@ -216,10 +220,21 @@ public:
      *
      * \param globalPos The global position
      */
-    PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+    using ParentType::dirichlet;
+    PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
     {
         PrimaryVariables values(0.0);
-        values = initialAtPos(globalPos);
+        values = initialAtPos(scvf.center());
+
+        const auto faceId = scvf.index();
+        if( couplingInterface_.isCoupledEntity( faceId ) )
+        {
+          values[Indices::velocityYIdx] =
+              couplingInterface_.getScalarQuantityOnFace( velocityId_, faceId );
+        }
+
+
+
         return values;
     }
 
@@ -253,7 +268,9 @@ public:
         if( couplingInterface_.isCoupledEntity( faceId ) )
         {
           const Scalar density = 1000; // TODO how to handle compressible fluids?
-          values[Indices::conti0EqIdx] = elemFaceVars[scvf].velocitySelf() * scvf.directionSign() * density;
+          const auto& volVars = elemVolVars[scvf.insideScvIdx()];
+          const Scalar density_ = volVars.density();
+          values[Indices::conti0EqIdx] = density * elemFaceVars[scvf].velocitySelf() * scvf.directionSign();
           values[Indices::momentumYBalanceIdx] = scvf.directionSign() * (couplingInterface_.getScalarQuantityOnFace( pressureId_, faceId ) - initialAtPos(scvf.center())[Indices::pressureIdx]) ;
         }
 #endif
diff --git a/appl/coupling-ff-pm/monolithic/params.input b/appl/coupling-ff-pm/monolithic/params.input
index 4117b7f17c501817b738dfdb6836cc2c3223e4a8..455dd8fbc09300fdc031e51ed0e4a2cb35e34f96 100644
--- a/appl/coupling-ff-pm/monolithic/params.input
+++ b/appl/coupling-ff-pm/monolithic/params.input
@@ -27,11 +27,12 @@ Grading1 = 1
 
 [Stokes.Problem]
 Name = stokes
-PressureDifference = 1e-9
-EnableInertiaTerms = false
+PressureDifference = 1e2
+EnableInertiaTerms = true
 
 [Darcy.Problem]
 Name = darcy
+InitialP = 0.5e-9
 
 [Darcy.SpatialParams]
 Permeability = 1e-6 # m^2
diff --git a/appl/coupling-ff-pm/monolithic/pmproblem.hh b/appl/coupling-ff-pm/monolithic/pmproblem.hh
index 56a1a4b9a0897c32bb75b905ee2962d407463075..40e7d7f153f3f4f50fa84ee6cb9a67c98e073694 100644
--- a/appl/coupling-ff-pm/monolithic/pmproblem.hh
+++ b/appl/coupling-ff-pm/monolithic/pmproblem.hh
@@ -170,6 +170,10 @@ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
         // set the coupling boundary condition at the interface
         if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
             values.setAllCouplingNeumann();
+#else
+        const auto faceId = scvf.index();
+        if ( couplingInterface_.isCoupledEntity(faceId) )
+          values.setAllDirichlet();
 #endif
         return values;
     }
@@ -188,6 +192,10 @@ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
         PrimaryVariables values(0.0);
         values = initial(element);
 
+        const auto faceId = scvf.index();
+        if ( couplingInterface_.isCoupledEntity(faceId) )
+          values = couplingInterface_.getScalarQuantityOnFace( pressureId_, faceId );
+
         return values;
     }
 
@@ -221,6 +229,7 @@ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
         {
           const Scalar density = 1000.;
           values[Indices::conti0EqIdx] = density * couplingInterface_.getScalarQuantityOnFace( velocityId_, faceId );
+          std::cout << "pm: values[Indices::conti0EqIdx] = " << values << std::endl;
         }
 #endif
         return values;
diff --git a/doc/mkdocs/dumux_precice/docs/dumux_precice_coupling.md b/doc/mkdocs/dumux_precice/docs/dumux_precice_coupling.md
index ae330280e54240d1238ed9776d0d3e5ae96a32b7..11337ac8cff194e157129f1b7d856b07fbc88239 100644
--- a/doc/mkdocs/dumux_precice/docs/dumux_precice_coupling.md
+++ b/doc/mkdocs/dumux_precice/docs/dumux_precice_coupling.md
@@ -17,12 +17,118 @@ The data is then incorporated into the BJS boundary conditions. Hopefully in a s
 
 - Let's try a `shared_ptr` concept instead of the singleton approach!
     - Does not prevent user to create more than one instance what is somewhat bad!
-    - Does not really work. So forget about that now. Cannot pass it properly to our hacky helper functions.
+    - **Does not really work.** So forget about that now. Cannot pass it properly to our hacky helper functions.
 
 - We need to set the initial pressure in the Darcy problem back to 0.
 
 #### Questions
 
+Flow is reversed after 1.5 iterations!
+
+**Initial velocity**
+![Initial velocity](./figs/initial_data_pm-ff-coupling.png)
+**Velocity after 1 ff step and 1 pm step (=1 coupling step)**
+![Velocity after 1 coupling step](./figs/1step_pm-ff-coupling.png)
+**Velocity after 2 ff steps and 1 pm step (=1.5 coupling step)**
+![Velocity after 2 ff steps and 1 pm step](./figs/1andhalf_step_pm-ff-coupling.png)
+
+**FF solver**
+```c++
+# 1st iteration
+velocities to be sent to ff
+v[0]=-6.61414e-08
+v[1]=-6.19115e-08
+v[2]=-5.6572e-08
+v[3]=-5.04326e-08
+v[4]=-4.36579e-08
+v[5]=-3.63742e-08
+v[6]=-2.86902e-08
+v[7]=-2.07036e-08
+v[8]=-1.25056e-08
+v[9]=-4.18234e-09
+v[10]=4.18234e-09
+v[11]=1.25056e-08
+v[12]=2.07036e-08
+v[13]=2.86902e-08
+v[14]=3.63742e-08
+v[15]=4.36579e-08
+v[16]=5.04326e-08
+v[17]=5.6572e-08
+v[18]=6.19115e-08
+v[19]=6.61414e-08
+
+# 2nd iteration
+velocities to be sent to ff
+v[0]=-2.41576e-05
+v[1]=-2.28153e-05
+v[2]=-2.1064e-05
+v[3]=-1.897e-05
+v[4]=-1.65748e-05
+v[5]=-1.39209e-05
+v[6]=-1.10525e-05
+v[7]=-8.01602e-06
+v[8]=-4.8588e-06
+v[9]=-1.62892e-06
+v[10]=1.62521e-06
+v[11]=4.85508e-06
+v[12]=8.01227e-06
+v[13]=1.10487e-05
+v[14]=1.3917e-05
+v[15]=1.65709e-05
+v[16]=1.8966e-05
+v[17]=2.106e-05
+v[18]=2.28111e-05
+v[19]=2.41533e-05
+```
+
+**PM solver**
+```c++
+# initial data
+p[0]=5e-10
+p[1]=5e-10
+p[2]=5e-10
+p[3]=5e-10
+p[4]=5e-10
+p[5]=5e-10
+p[6]=5e-10
+p[7]=5e-10
+p[8]=5e-10
+p[9]=5e-10
+p[10]=5e-10
+p[11]=5e-10
+p[12]=5e-10
+p[13]=5e-10
+p[14]=5e-10
+p[15]=5e-10
+p[16]=5e-10
+p[17]=5e-10
+p[18]=5e-10
+p[19]=5e-10
+
+# 1st iteration
+p[0]=-1.6348e-06
+p[1]=-1.52962e-06
+p[2]=-1.39717e-06
+p[3]=-1.24513e-06
+p[4]=-1.07754e-06
+p[5]=-8.97521e-07
+p[6]=-7.07721e-07
+p[7]=-5.10538e-07
+p[8]=-3.08192e-07
+p[9]=-1.02794e-07
+p[10]=1.03616e-07
+p[11]=3.09015e-07
+p[12]=5.1136e-07
+p[13]=7.08544e-07
+p[14]=8.98344e-07
+p[15]=1.07837e-06
+p[16]=1.24595e-06
+p[17]=1.39799e-06
+p[18]=1.53044e-06
+p[19]=1.63562e-06
+```
+
+
 In FreeFlow `ffproblem.hh`
 
 - Should alphaBJ and permeability be transfered? I thought they would be constant?
@@ -33,6 +139,20 @@ In Darcy
 - How to reconstruct certain things?!
 - What happens when `values.setBJS` is called? Is it callable with preCICE running?
 
+**NEW**
+
+- I even get negative pressures. Should that be expected? I would guess pressures would be in between 0 and 1e-9 as prescribed on the domain boundaries. Although I understand that is not total pressure so negative values are not forbidden.
+- What are the pressures in the Darcy domain normally? Is there a jump between free flow and Darcy flow?
+- Change direction of data transfer to change the Dirichlet/Neumann coupling? The idea is to give the Darcy-solver more freedom when solving its problem. 
+
+| Solver | writes | reads |
+| ------ | ------ | ----- |
+| FreeFlow | Pressure | Velocity |
+| Darcy | Velocity | Pressure |
+
+- Monolithic solver is broken?!
+
+
 
 ### Conjugate heat transfer problem 
 
@@ -40,9 +160,8 @@ See [test case description](dumux_test_case.md#conjugate-heat-transfer)
 
 
 We do the following now:
-
 | Solver | writes | reads |
-| ------ | ------ | ----- |
+| ------------ | ------ | ----- |
 | SolidEnergy | Temperature | Heat flux |
 | FreeFlow | Heat flux | Temperature |
 
diff --git a/doc/mkdocs/dumux_precice/docs/index.md b/doc/mkdocs/dumux_precice/docs/index.md
index e3d9e0a5651af241840ab6ed7b3f6585e00f0bdc..28816cd64e2ad2688476133e007a82dd2a2cf244 100644
--- a/doc/mkdocs/dumux_precice/docs/index.md
+++ b/doc/mkdocs/dumux_precice/docs/index.md
@@ -4,3 +4,5 @@
 
 - GitLab: [https://git.iws.uni-stuttgart.de/dumux-appl/dumux-precice](https://git.iws.uni-stuttgart.de/dumux-appl/dumux-precice)
 - Logs/Minutes of last meetings can be found [here](logs/logs.md)
+- The DuMuX adapter moved into [another repository](git@gitlab-sgs.informatik.uni-stuttgart.de:jaustar/dumux-precice-wrapper.git) and was added as a submodule.
+