diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index 02cac8eabdff8d7406cc771d674f0a98a8f39717..f71ee9fc7605e1f06ff347585f606a543ef6945f 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -50,6 +50,7 @@ void BLRegion::SetMesh(std::vector<double> locM, std::vector<double> globM)
   XxExt.resize(nMarkers);
   ReExt.resize(nMarkers);
   Regime.resize(nMarkers);
+  Blowing.resize(nMarkers);
 } // End SetMesh
 
 void BLRegion::SetExtVariables(std::vector<double> _xxExt, std::vector<double> _deltaStarExt)
@@ -138,6 +139,17 @@ void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<d
   }
 }
 
+void BLRegion::ComputeBlwVel()
+{
+  for (size_t iPoint=1; iPoint<nMarkers; ++iPoint){
+
+    DeltaStar[iPoint] = U[iPoint*nVar+0] * U[iPoint*nVar+1]; 
+    Blowing[iPoint-1] = (Rhoe[iPoint]*VtExt[iPoint]*DeltaStar[iPoint]
+    -Rhoe[iPoint-1]*VtExt[iPoint-1]*DeltaStar[iPoint-1])/(Rhoe[iPoint]*(LocMarkers[iPoint]-LocMarkers[iPoint-1]));
+    
+  }
+}
+
 unsigned int BLRegion::GetnVar(){return nVar;}
 size_t BLRegion::GetnMarkers(){return nMarkers;}
 double BLRegion::GetnCrit(){return nCrit;}
diff --git a/dart/src/wBLRegion.h b/dart/src/wBLRegion.h
index 8d818757e311db4e8d37cb1aaa39ac897819a624..d3383c4ab8f7f41c97ce3b243150acc6371371f2 100644
--- a/dart/src/wBLRegion.h
+++ b/dart/src/wBLRegion.h
@@ -59,6 +59,7 @@ public:
   /* Transition related variables */
 
   std::vector<int> Regime;
+  std::vector<double> Blowing;
 
   unsigned int transMarker;
   double xtrF;
@@ -86,7 +87,11 @@ public:
   double GetXxExt(size_t iPoint);
   double GetDeltaStarExt(size_t iPoint);
 
+
   void PrintSolution(size_t iPoint);
+  void ComputeBlwVel();
+
+private:
 };
 } // namespace dart
 #endif //WBLREGION_H
\ No newline at end of file
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index dcdd8877b7099e8f9b7e255831572da1fdfb92e8..7af50a9e0da44aee5a3368f05989dd107ca499d7 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -53,9 +53,9 @@ void TimeSolver::InitialCondition(size_t iPoint, BLRegion *bl)
     bl->Delta[iPoint] = bl->Delta[iPoint - 1];
     bl->DeltaStar[iPoint] = bl->DeltaStar[iPoint - 1];
   }
-    if (bl->Regime[iPoint] == 1 && bl->U[iPoint * nVar + 4] == 0){
-        bl->U[iPoint * nVar + 4] = 0.03;
-    } // End if
+  if (bl->Regime[iPoint] == 1 && bl->U[iPoint * nVar + 4] == 0){
+      bl->U[iPoint * nVar + 4] = 0.03;
+  } // End if
 } // End InitialCondition
 
 int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
@@ -92,27 +92,15 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
   unsigned int nErrors = 0; // Number of errors encountered
   unsigned int innerIters = 0; // Inner (non-linear) iterations
 
-  if(bl->name == "upper" && iPoint == 44){
-      std::cout << "IC 44 " << std::endl;
-      bl->PrintSolution(42);
-      }
-
   while (innerIters < maxIt){
 
     if(innerIters % itFrozenJac == 0){
       JacMatrix = SysMatrix->ComputeJacobian(iPoint, bl, SysRes, 1e-8);
-      if(bl->name == "upper" && iPoint == 44){
-      std::cout << "computeJac at iter " << innerIters  << JacMatrix << std::endl;
-      std::cout << "SysRes " << SysRes << std::endl;
-      }
     }
-    
+
     // Ici on peut faire autrement : .asDiagonal()
     slnIncr = (JacMatrix + IMatrix).partialPivLu().solve(-SysRes);
 
-    if(bl->name == "upper" && iPoint == 44){
-    std::cout << "slnIncr " << slnIncr << std::scientific;}
-
     for(size_t k=0; k<nVar; ++k){
       bl->U[iPoint*nVar+k] += slnIncr(k);
     }
@@ -120,14 +108,9 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
     SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
     normSysRes = (SysRes + slnIncr/timeStep).norm();
 
-    if(bl->name == "upper" && iPoint == 44){
-      std::cout << "normsysRes/nromsysRes0" << normSysRes/normSysRes0 << " CFL"  << CFL << "TimeStep" << timeStep << std::endl;
-      bl->PrintSolution(iPoint);
-    }
-
     if (normSysRes <= tol * normSysRes0){
       /* Successfull exit */
-      std::cout << "Point" << iPoint << " Convergence successfull. Niter: " << innerIters << " NormsysRes" << normSysRes << std::endl;
+      //std::cout << "Point" << iPoint << " Convergence successfull. Niter: " << innerIters << " NormsysRes" << normSysRes << std::endl;
       return 0;
     }
 
@@ -141,7 +124,7 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
 
     innerIters++;
   } // End while (innerIters < maxIt)
-  std::cout << "Point" << iPoint << " Max iterations number reached. normSysRes/normSysRes0 " << normSysRes/normSysRes0 << std::endl;
+  //std::cout << "Point" << iPoint << " Max iterations number reached. normSysRes/normSysRes0 " << normSysRes/normSysRes0 << std::endl;
   return innerIters;
 } // End Integration
 
diff --git a/dart/src/wTimeSolver.h b/dart/src/wTimeSolver.h
index e89628937f3aad81f6d0d9409c130c400653f4b7..e676590293fedd701dafa11dcc67553f6fbe12db 100644
--- a/dart/src/wTimeSolver.h
+++ b/dart/src/wTimeSolver.h
@@ -18,14 +18,14 @@ private:
   unsigned long maxIt;
   double tol;
   unsigned int itFrozenJac0;
-  bool initSln = true;
+  bool initSln;
   double gamma = 1.4;
   double Minf;
 
   ViscFluxes *SysMatrix;
 
 public:
-  TimeSolver(double _CFL0, double _Minf, double Re, unsigned long _maxIt = 10, double _tol = 1e-8, unsigned int _itFrozenJac = 5);
+  TimeSolver(double _CFL0, double _Minf, double Re, unsigned long _maxIt = 50, double _tol = 1e-8, unsigned int _itFrozenJac = 5);
   ~TimeSolver();
   void InitialCondition(size_t iPoint, BLRegion *bl);
   int Integration(size_t iPoint, BLRegion *bl);
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index 3a66e559520b11cc8d5000a89bd36467c08b1814..b99a85ece9388de047c40bacfef4d4748b859134 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -157,6 +157,19 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
 
   /* Space part */
 
+  /*if(bl->name=="lower" && iPoint==1){
+    std::cout << "due_dx" << due_dx << std::endl;
+    std::cout << "c" << c << std::endl;
+    std::cout << "u[1]" << u[1] << std::endl;
+    std::cout << "dT_dx" << dT_dx << std::endl;
+    std::cout << "u[0]" << u[0] << std::endl;
+    std::cout << "dH_dx" << dH_dx << std::endl;
+    std::cout << "dueExt_dx" << dueExt_dx << std::endl;
+    std::cout << "cExt" << cExt << std::endl;
+    std::cout << "ddStarExt_dx" << ddeltaStarExt_dx << std::endl;
+    std::cout << "dHstar_dx" << dHstar_dx << std::endl;
+  }*/
+
   spaceVector(0) = dT_dx + (2 + u[1] - Mea * Mea) * (u[0]/u[3]) * due_dx - Cfa/2;
   spaceVector(1) = u[0] * dHstar_dx + (2*Hstar2a + Hstara * (1-u[1])) * u[0]/u[3] * due_dx - 2*Cda + Hstara * Cfa/2;
   spaceVector(2) = dN_dx - Ax;
@@ -209,7 +222,7 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
     timeMatrix(4,4) = 1.0;
   }
 
-  /*if(bl->name=="upper" && iPoint==1){
+  /*if(bl->name=="lower" && iPoint==1){
     std::cout << "timeMatrix" << timeMatrix << std::endl;
     std::cout << "spaceVector" << spaceVector << std::endl;
 
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 09db6df85b8de976eb8bc4f9ac4187bb4858052a..633378a5b050f8b20c09b8be663fcd3bf596b3f6 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -86,6 +86,12 @@ int ViscSolver::Run(unsigned int couplIter){
   if (couplIter == 0){
     tSolver->SetCFL0(1);
     tSolver->SetitFrozenJac(1);
+    tSolver->SetinitSln(true);
+  }
+  else{
+    tSolver->SetCFL0(1);
+    tSolver->SetitFrozenJac(1);
+    tSolver->SetinitSln(false);
   }
 
   /* Set boundary condition */
@@ -152,9 +158,14 @@ int ViscSolver::Run(unsigned int couplIter){
       } // End if (!lockTrans).
 
     } // End for iPoints
-  } // End for iRegions
+  } // End for iRegion
 
   ComputeDrag();
+
+  for (size_t iRegion=0; iRegion<bl.size(); ++iRegion){
+    bl[iRegion]->ComputeBlwVel();
+  } // End for iRegion
+
   return 0; // Successfull exit
 }
 
@@ -241,6 +252,23 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
   bl->Us[iPoint] = turbParam[7];
 }
 
+std::vector<double> ViscSolver::ExtractBlowingVel(size_t iRegion)
+{
+  return bl[iRegion]->Blowing;
+}
+std::vector<double> ViscSolver::ExtractDeltaStar(size_t iRegion)
+{
+  return bl[iRegion]->DeltaStar;
+}
+std::vector<double> ViscSolver::ExtractXx(size_t iRegion)
+{
+  std::vector<double> LocMarkersOut;
+  for (size_t k=0; k<bl[iRegion]->GetnMarkers(); ++k){
+    LocMarkersOut.push_back(bl[iRegion]->GetLocMarkers(k));
+  }
+  return LocMarkersOut;
+}
+
 
 void ViscSolver::ComputeDrag()
 {
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
index 0338f9fe7dd866f4bfbbf5d1952e096de3680c2b..0cf54920e1e85d6eff0b5d78dfad3b2cce9630cb 100644
--- a/dart/src/wViscSolver.h
+++ b/dart/src/wViscSolver.h
@@ -33,13 +33,17 @@ public:
                                std::vector<double> _Rhoe);
   void SendExtVariables(size_t region, std::vector<double> _xxExt, std::vector<double> _deltaStarExt);
 
-
   int Run(unsigned int couplIter);
+  std::vector<double> ExtractBlowingVel(size_t iRegion);
+  std::vector<double> ExtractDeltaStar(size_t iRegion);
+  std::vector<double> ExtractXx(size_t iRegion);
+
 
 private:
   void CheckWaves(size_t iPoint, BLRegion *bl);
   void AverageTransition(size_t iPoint, BLRegion *bl, int forced);
   void ComputeDrag();
+  void ComputeBlwVel();
 
 };
 } // namespace dart
diff --git a/dart/vCoupler.py b/dart/vCoupler.py
index 75253b5598709e1cdde173133219624561c8348e..18a07c0a00f8ea27612f39bb488d4e39a501285a 100644
--- a/dart/vCoupler.py
+++ b/dart/vCoupler.py
@@ -26,19 +26,23 @@ from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
 from fwk.coloring import ccolors
 
 import numpy as np
+from matplotlib import pyplot as plt
 
 class Coupler:
     def __init__(self, iSolver, vSolver) -> None:
         self.iSolver=iSolver
         self.vSolver=vSolver
 
-        self.maxCouplIter = 10
+        self.maxCouplIter = 100
         self.couplTol = 1e-4
         pass
 
     def CreateProblem(self, _boundaryAirfoil, _boundaryWake):
 
         self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)] # airfoil and wake python objects
+        self.blwVel = [None, None, None]
+        self.deltaStar = [None, None, None]
+        self.xx = [None, None, None]
         pass
 
     def Run(self):
@@ -70,6 +74,8 @@ class Coupler:
         for n in range(0,len(self.group)):
           self.group[n].connectListNodes, self.group[n].connectListElems, dataIn[n]  = self.group[n].connectList()
         fMeshDict, cMeshDict, data = GroupConstructor().mergeVectors(dataIn, 0, 1)
+        #print('group[1].deltaStar', self.group[1].deltaStar)
+        #print('group[1].xx', self.group[1].xx)
         if couplIter == 0:
           # Initialize mesh          
           self.vSolver.InitMeshes(0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]], fMeshDict['globCoord'][:fMeshDict['bNodes'][1]])
@@ -108,12 +114,31 @@ class Coupler:
         .format(couplIter, self.vSolver.Cdt, np.log10(error), np.log10(self.couplTol)), ccolors.ANSI_RESET)
         print('')
 
-        # Set inviscid boundary condition
+        # Retreive and set blowing velocities.
 
-        for n in range(0, len(self.group)):
+        for iRegion in range(3):
+          self.blwVel[iRegion] = np.asarray(self.vSolver.ExtractBlowingVel(iRegion))
+          self.deltaStar[iRegion] = np.asarray(self.vSolver.ExtractDeltaStar(iRegion))
+          self.xx[iRegion] = np.asarray(self.vSolver.ExtractXx(iRegion))
 
-          for i in range(0, self.group[n].nE):
+        blwVelAF = np.concatenate((self.blwVel[0], self.blwVel[1][1:])) # Remove stagnation point doublon.
+        deltaStarAF = np.concatenate((self.deltaStar[0], self.deltaStar[1][1:])) # Remove stagnation point doublon.
+        xxAF = np.concatenate((self.xx[0], self.xx[1][1:])) # Remove stagnation point doublon.
+        blwVelWK = self.blwVel[2]
+        deltaStarWK = self.deltaStar[2]
+        xxWK = self.xx[2]
+
+        self.group[0].u = blwVelAF[self.group[0].connectListElems.argsort()]
+        self.group[1].u = blwVelWK[self.group[1].connectListElems.argsort()]
+
+        self.group[0].deltaStar = deltaStarAF[self.group[0].connectListNodes.argsort()]
+        self.group[1].deltaStar = deltaStarWK[self.group[1].connectListNodes.argsort()]
 
+        self.group[0].xx = xxAF[self.group[0].connectListNodes.argsort()]
+        self.group[1].xx = xxWK[self.group[1].connectListNodes.argsort()]
+
+        for n in range(0, len(self.group)):
+          for i in range(0, self.group[n].nE):
             self.group[n].boundary.setU(i, self.group[n].u[i])
 
         couplIter += 1
diff --git a/dart/viscous/Drivers/DDriver.py b/dart/viscous/Drivers/DDriver.py
index eb4e4b2a88aa3ca4caf2f9f4f4154aa5b608dc8b..99367faee76015089bc2c727e1e014630f0d8730 100644
--- a/dart/viscous/Drivers/DDriver.py
+++ b/dart/viscous/Drivers/DDriver.py
@@ -14,7 +14,7 @@
 # ------------------------------------------------------------------------------
 #  Imports
 # ------------------------------------------------------------------------------
-from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
+from dart.viscous.Drivers.DGroupConstructor2 import GroupConstructor
 from dart.viscous.Drivers.DPostProcess import PostProcess
 
 from dart.viscous.Solvers.DIntegration import Integration
diff --git a/dart/viscous/Drivers/DGroupConstructor2.py b/dart/viscous/Drivers/DGroupConstructor2.py
new file mode 100755
index 0000000000000000000000000000000000000000..44877c7d62e24180c383d8319156baa755c36ed9
--- /dev/null
+++ b/dart/viscous/Drivers/DGroupConstructor2.py
@@ -0,0 +1,190 @@
+################################################################################
+#                                                                              #
+#                             Flow viscous module file                         #
+#                        Constructor : Merges and sorts data                   #
+#                   of the 3 groups (upper, lower sides and wake)              #
+#                                                                              #
+# Author: Paul Dechamps                                                        #
+# Université de Liège                                                          #
+# Date: 2021.09.10                                                             #
+#                                                                              #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+#  Imports
+# ------------------------------------------------------------------------------
+from matplotlib import pyplot as plt
+import numpy as np
+
+class GroupConstructor():
+    """ ---------- Group Constructor ----------
+        
+        Builds single vectors elements out of the 3 groups given in input.
+    """
+    def __init__(self) -> None:
+        pass
+
+    def mergeVectors(self,dataIn, mapMesh, nFactor):
+        """Merges 3 groups upper/lower sides and wake.
+        """
+        
+        for k in range(len(dataIn)):
+
+            if len(dataIn[k]) == 2: # Airfoil case.
+            
+                xx_up, dv_up, vtExt_up, _, alpha_up = self.__getBLcoordinates(dataIn[k][0])
+                xx_lw, dv_lw, vtExt_lw, _, alpha_lw = self.__getBLcoordinates(dataIn[k][1])
+
+            else: # Wake case
+
+                xx_wk, dv_wk, vtExt_wk, _, alpha_wk = self.__getBLcoordinates(dataIn[k])
+
+        if mapMesh:
+            # Save initial mesh.
+            cMesh        = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
+            cMeshbNodes  = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
+            cxx          = np.concatenate((xx_up, xx_lw[1:], xx_wk))
+            cvtExt       = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk))
+            cxxExt       = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][1:,8],dataIn[1][:,8]))
+
+            # Create fine mesh.
+            fMeshUp      = self.createFineMesh(dataIn[0][0][:,0], nFactor)
+            fMeshLw      = self.createFineMesh(dataIn[0][1][:,0], nFactor)
+            fMeshWk      = self.createFineMesh(dataIn[1][:,0]   , nFactor)
+            fMesh        = np.concatenate((fMeshUp, fMeshLw[1:], fMeshWk))
+            fMeshbNodes  = [0, len(fMeshUp), len(fMeshUp) + len(fMeshLw[1:])]
+
+            fxx_up       = self.interpolateToFineMesh(xx_up, fMeshUp, nFactor)
+            fxx_lw       = self.interpolateToFineMesh(xx_lw, fMeshLw, nFactor)
+            fxx_wk       = self.interpolateToFineMesh(xx_wk, fMeshWk, nFactor)
+            fxx          = np.concatenate((fxx_up, fxx_lw[1:], fxx_wk))
+
+            fvtExt_up    = self.interpolateToFineMesh(vtExt_up, fMeshUp, nFactor)
+            fvtExt_lw    = self.interpolateToFineMesh(vtExt_lw, fMeshLw, nFactor)
+            fvtExt_wk    = self.interpolateToFineMesh(vtExt_wk, fMeshWk, nFactor)
+            fvtExt       = np.concatenate((fvtExt_up, fvtExt_lw[1:], fvtExt_wk))
+
+            fMe_up       = self.interpolateToFineMesh(dataIn[0][0][:,5], fMeshUp, nFactor)
+            fMe_lw       = self.interpolateToFineMesh(dataIn[0][1][:,5], fMeshLw, nFactor)
+            fMe_wk       = self.interpolateToFineMesh(dataIn[1][:,5]   , fMeshWk, nFactor)
+            fMe          = np.concatenate((fMe_up, fMe_lw[1:], fMe_wk))
+
+            frho_up      = self.interpolateToFineMesh(dataIn[0][0][:,6], fMeshUp, nFactor)
+            frho_lw      = self.interpolateToFineMesh(dataIn[0][1][:,6], fMeshLw, nFactor)
+            frho_wk      = self.interpolateToFineMesh(dataIn[1][:,6]   , fMeshWk, nFactor)
+            frho         = np.concatenate((frho_up, frho_lw[1:], frho_wk))
+
+            fdStarExt_up = self.interpolateToFineMesh(dataIn[0][0][:,7], fMeshUp, nFactor)
+            fdStarExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,7], fMeshLw, nFactor)
+            fdStarExt_wk = self.interpolateToFineMesh(dataIn[1][:,7]   , fMeshWk, nFactor)
+            
+            fxxExt_up    = self.interpolateToFineMesh(dataIn[0][0][:,8], fMeshUp, nFactor)
+            fxxExt_lw    = self.interpolateToFineMesh(dataIn[0][1][:,8], fMeshLw, nFactor)
+            fxxExt_wk    = self.interpolateToFineMesh(dataIn[1][:,8]   , fMeshWk, nFactor)
+            
+            fdStarext    = np.concatenate((fdStarExt_up, fdStarExt_lw[1:], fdStarExt_wk))
+            fxxext       = np.concatenate((fxxExt_up, fxxExt_lw[1:], fxxExt_wk))
+
+            fdv          = np.concatenate((dv_up, dv_lw[1:], dv_wk))
+
+            fMeshDict    = {'globCoord': fMesh, 'bNodes': fMeshbNodes, 'locCoord': fxx,
+                            'vtExt': fvtExt, 'Me': fMe, 'rho': frho, 'deltaStarExt': fdStarext, 'xxExt': fxxext, 'dv': fdv}
+
+            cMe          = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][1:,5], dataIn[1][:,5]))
+            cRho         = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][1:,6], dataIn[1][:,6]))
+            cdStarExt    = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][1:,7], dataIn[1][:,7]))
+            cMeshDict    = {'globCoord': cMesh, 'bNodes': cMeshbNodes, 'locCoord': cxx,
+                            'vtExt': cvtExt, 'Me': cMe, 'rho': cRho, 'deltaStarExt': cdStarExt, 'xxExt': cxxExt, 'dv': fdv}
+
+            dataUpper = np.empty((len(fMeshUp), 0))
+            dataLower = np.empty((len(fMeshLw), 0))
+            dataWake  = np.empty((len(fMeshWk), 0))
+            for iParam in range(len(dataIn[0][0][0,:])):
+                interpParamUp = self.interpolateToFineMesh(dataIn[0][0][:,iParam], fMeshUp, nFactor)
+                dataUpper     = np.column_stack((dataUpper, interpParamUp))
+                interpParamLw = self.interpolateToFineMesh(dataIn[0][1][:,iParam], fMeshLw, nFactor)
+                dataLower     = np.column_stack((dataLower, interpParamLw))
+                interpParamWk = self.interpolateToFineMesh(dataIn[1][:,iParam]   , fMeshWk, nFactor)
+                dataWake      = np.column_stack((dataWake, interpParamWk))
+            
+            dataLower = np.delete(dataLower, 0, axis = 0) # Remove stagnation point doublon. 
+
+        else:
+
+            Mesh         = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
+            MeshbNodes   = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
+
+            xx           = np.concatenate((xx_up, xx_lw[1:], xx_wk))
+            vtExt        = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk))
+
+            alpha        = np.concatenate((alpha_up, alpha_lw, alpha_wk))
+            dv           = np.concatenate((dv_up, dv_lw[1:], dv_wk))
+
+            Me           = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][1:,5], dataIn[1][:,5]))
+            rho          = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][1:,6], dataIn[1][:,6]))
+            dStarext     = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][1:,7], dataIn[1][:,7]))
+            xxExt        = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][1:,8], dataIn[1][:,8]))
+
+            fMeshDict    = {'globCoord': Mesh, 'bNodes': MeshbNodes, 'locCoord': xx,
+                            'vtExt': vtExt, 'Me': Me, 'rho': rho, 'deltaStarExt': dStarext, 'xxExt': xxExt, 'dv': dv}
+
+            cMeshDict    = fMeshDict.copy()
+
+            dataUpper    = dataIn[0][0]
+            dataLower    = dataIn[0][1][1:, :]
+            dataWake     = dataIn[1]
+        
+        data = np.concatenate((dataUpper, dataLower, dataWake), axis = 0)
+
+        return fMeshDict, cMeshDict, data
+
+    def __getBLcoordinates(self,data):
+        """Transform the reference coordinates into airfoil coordinates
+        """
+        nN    = len(data[:,0])
+        nE    = nN-1            # Nbr of element along the surface.
+        vt    = np.zeros(nN)    # Outer flow velocity. 
+        xx    = np.zeros(nN)    # Position along the surface of the airfoil.
+        dx    = np.zeros(nE)    # Distance along two nodes.
+        dv    = np.zeros(nE)    # Speed difference according to element.
+        alpha = np.zeros(nE)    # Angle of the element for the rotation matrix.
+
+        for i in range(0, nE):
+
+            alpha[i] = np.arctan2((data[i+1, 1] - data[i, 1]), (data[i+1, 0]-data[i, 0]))
+            dx[i]    = np.sqrt((data[i+1, 0] - data[i,0])**2 + (data[i+1, 1]-data[i, 1])**2)
+            xx[i+1]  = dx[i] + xx[i]
+            vt[i]    = (data[i, 3] * np.cos(alpha[i]) + data[i ,4] * np.sin(alpha[i]))          # Tangential speed at node 1 element i
+            vt[i+1]  = (data[i+1, 3] * np.cos(alpha[i]) + data[i+1, 4] * np.sin(alpha[i]))      # Tangential speed at node 2 element i
+            dv[i]    = (vt[i+1] - vt[i]) / dx[i]                                                # Velocity gradient according to element i.
+                                                                                                #             central scheme with half point
+        
+        return xx, dv, vt, nN, alpha
+
+    def interpolateToFineMesh(self, cVector, fMesh, nFactor):
+        
+        fVector = np.empty(len(fMesh)-1)
+        for cPoint in range(len(cVector) - 1):
+            dv = cVector[cPoint + 1] - cVector[cPoint]
+
+            for fPoint in range(nFactor):
+                fVector[nFactor * cPoint + fPoint] = cVector[cPoint] + fPoint * dv / nFactor
+        
+        fVector = np.append(fVector, cVector[-1])
+
+        return fVector
+
+    def createFineMesh(self, cMesh, nFactor):
+
+        fMesh = []
+        for iPoint in range(len(cMesh) - 1):
+
+            dx = cMesh[iPoint + 1] - cMesh[iPoint]
+
+            for iRefinedPoint in range(nFactor):
+
+                fMesh = np.append(fMesh, cMesh[iPoint] + iRefinedPoint * dx / nFactor)
+
+        fMesh = np.append(fMesh, cMesh[-1])
+        return fMesh
\ No newline at end of file
diff --git a/dart/viscous/Master/DCoupler.py b/dart/viscous/Master/DCoupler.py
index 4d83869c771a6eec9dbaf8dca1f3a9adbbe7eefa..780e4cc676527ab8d8b2918a523ae65ec706d5ea 100755
--- a/dart/viscous/Master/DCoupler.py
+++ b/dart/viscous/Master/DCoupler.py
@@ -54,7 +54,7 @@ class Coupler:
         print('+---------------------------- VII Solver begin --------------------------------+')
         print('')
 
-        while converged == 0:
+        while it < 100:
 
             # Run inviscid solver
             print('+----------------------------- Inviscid solver --------------------------------+')
@@ -90,7 +90,7 @@ class Coupler:
 
             # Check convergence and output coupling iteration.
 
-            if error <= self.tol:
+            if 1 <= 0:
               print('')
               print(ccolors.ANSI_GREEN, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} < {:<2.3f}'
               .format(it, self.vsolver.Cd, np.log10(error), np.log10(self.tol)), ccolors.ANSI_RESET)
diff --git a/dart/viscous/Solvers/DIntegration.py b/dart/viscous/Solvers/DIntegration.py
index 4c97ee0622c8ca7bf34d4f6b135b575d1252651a..624b9551b08ec6c7c78d5bc325c651d04bd39e82 100644
--- a/dart/viscous/Solvers/DIntegration.py
+++ b/dart/viscous/Solvers/DIntegration.py
@@ -41,7 +41,7 @@ class Integration:
     self.sysMatrix = _sysMatrix
 
     self.__CFL0=_Cfl0
-    self.__maxIt = 1000
+    self.__maxIt = 50
     self.__tol = 1e-8
     self.itFrozenJac0 = 5
     
@@ -109,8 +109,6 @@ class Integration:
     innerIters = 0                                                # Number of inner iterations to solver the non-linear system.
     while innerIters <= self.__maxIt:
       try:
-        if iMarker == 2:
-          quit()
         # Jacobian computation.
 
         if innerIters % itFrozenJac == 0:
@@ -129,7 +127,6 @@ class Integration:
 
         sysRes = self.sysMatrix.buildResiduals(cfgFlow, iMarker)
         normSysRes = np.linalg.norm(sysRes + slnIncr/timeStep)
-        print('CFL', CFL, 'timeStep', timeStep)
         convPlot.append(normSysRes/normSysRes0)
 
         # Check convergence.