diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index 8d001e8b4c5c4357a039c6840727aad942786380..2caf52a4c3c78e502ee26ff1c27ed61d0a00a75b 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -48,7 +48,7 @@ void BLRegion::SetMesh(std::vector<double> locM, std::vector<double> globM)
   XxExt.resize(nMarkers);
   ReExt.resize(nMarkers);
   Regime.resize(nMarkers);
-  Blowing.resize(nMarkers);
+  Blowing.resize(nMarkers-1);
 } // End SetMesh
 
 void BLRegion::SetExtVariables(std::vector<double> _xxExt, std::vector<double> _deltaStarExt)
diff --git a/dart/src/wBLRegion.h b/dart/src/wBLRegion.h
index d3383c4ab8f7f41c97ce3b243150acc6371371f2..244fde997b9d0d02d4cd7f459b4a9430c7c22307 100644
--- a/dart/src/wBLRegion.h
+++ b/dart/src/wBLRegion.h
@@ -40,6 +40,9 @@ public:
   Closures *closSolver;
 
   std::string name;
+  
+  size_t coarseSize; // Size of the coarse mesh on the corresponding region (only used if mesh mapping is enabled).
+  size_t fineSize; // Size of the coarse mesh on the corresponding region (only used if mesh mapping is enabled).
 
    /* Boundary layer variables */
 
diff --git a/dart/src/wSolutionInterp.cpp b/dart/src/wSolutionInterp.cpp
index ca13de009747b8501ca8d2d696b80dda27d2abb7..ba32bfb24687c21f45bb3283ad0b88cbe234cbf4 100644
--- a/dart/src/wSolutionInterp.cpp
+++ b/dart/src/wSolutionInterp.cpp
@@ -33,7 +33,7 @@ void SolutionInterp::InputFineMesh(std::vector<double> _fineLocMesh, std::vector
   FineGlobMesh=_fineGlobMesh;
 }
 
-std::vector<double> SolutionInterp::InterpolateToFineMesh(std::vector<double> cVector)
+std::vector<double> SolutionInterp::InterpolateToFineMesh(std::vector<double> cVector, size_t fineSize)
 {
   if (meshFactor<=1)
   {
@@ -42,31 +42,55 @@ std::vector<double> SolutionInterp::InterpolateToFineMesh(std::vector<double> cV
   }
   else
   {
-    std::vector<double> fVector;
+    std::vector<double> fVector(fineSize-1);
     double dv;
-    for (size_t cPoint=0; cPoint<cVector.size()-1; ++cPoint)
+    for (size_t cPoint=0; cPoint < cVector.size()-1; ++cPoint)
     {
-      dv = cVector[cPoint+1] - cVector[cPoint];
-
+      dv = cVector[cPoint + 1] - cVector[cPoint];
       for (size_t fPoint=0; fPoint<meshFactor; ++fPoint)
       {
-        fVector.push_back(cVector[cPoint] + fPoint * dv / meshFactor);
+        fVector[meshFactor * cPoint + fPoint] = cVector[cPoint] + fPoint * dv / meshFactor;
       }
     }
-    std::cout << cVector.back() << std::endl;
     fVector.push_back(cVector.back());
-    return fVector;
+    return fVector; 
   }
 }
 
+std::vector<double> SolutionInterp::BlowingFineToCoarse(std::vector<double> fVector, size_t coarseSize)
+{
+  std::vector<double> cVector(coarseSize);
+
+  for (size_t cCell=0; cCell<cVector.size(); ++cCell)
+  {
+    for (size_t k=0; k<meshFactor; ++k)
+    {
+      cVector[cCell] += fVector[cCell*meshFactor+k];
+    }
+    cVector[cCell]/=meshFactor;
+  }  
+  cVector.back() = fVector.back();
+  return cVector;
+}
+
+std::vector<double> SolutionInterp::VarsFineToCoarse(std::vector<double> fVector, size_t coarseSize)
+{
+  std::vector<double> cVector(coarseSize);
+
+  for (size_t cPoint=0; cPoint<cVector.size(); ++cPoint)
+  {
+    cVector[cPoint] = fVector[cPoint * meshFactor];
+  }
+  cVector.back() = fVector.back();
+  return cVector;
+}
+
+
 std::vector<double> SolutionInterp::CreateFineMesh(std::vector<double> cMesh)
 {
   if (meshFactor<1)
   {
-    std::cout << "SolutionInterp::CreateFineMesh: Mesh factor should be > 1." << std::endl;
-    /*std::stringstream err;
-    err << "SolutionInterp::CreateFineMesh: Mesh factor should be > 1. Value: " << meshFactor << "\n";
-    throw std::runtime_error(err.str());*/
+    throw std::runtime_error("SolutionInterp::CreateFineMesh: Mesh factor should be greater than 1.");
   }
 
   else
diff --git a/dart/src/wSolutionInterp.h b/dart/src/wSolutionInterp.h
index d023d57fed99e71dbceba2d6c625af526a514bd9..a16b2444ce057bdc54e09fa054e91c3c30865819 100644
--- a/dart/src/wSolutionInterp.h
+++ b/dart/src/wSolutionInterp.h
@@ -24,9 +24,11 @@ public:
   ~SolutionInterp();
 
   void InputFineMesh(std::vector<double> _fineLocMesh, std::vector<double> _fineGlobMesh);
+  std::vector<double> BlowingFineToCoarse(std::vector<double> fVector, size_t coarseSize);
+  std::vector<double> VarsFineToCoarse(std::vector<double> fVector, size_t coarseSize);
   std::vector<double> CreateFineMesh(std::vector<double> cMesh);
   void ResetmeshFactor(unsigned int newmeshFactor);
-  std::vector<double> InterpolateToFineMesh(std::vector<double> cVector);
+  std::vector<double> InterpolateToFineMesh(std::vector<double> cVector, size_t fineSize);
   void InputCoarseMesh(std::vector<double> _CoarseLocMesh, std::vector<double> _CoarseGlobMesh);
 };
 } // namespace dart
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index e0d17f706a6bb296dbd435f86d01478dc1d61a4e..66402b83734a2f42ed021d62ae7dc65e5eed8fe3 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -68,9 +68,9 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
 
   unsigned int nVar = bl->GetnVar();
 
-  std::vector<double> temp_U;
+  std::vector<double> Sln0;
   for (size_t i=0; i<nVar; ++i){
-    temp_U.push_back(bl->U[iPoint*nVar+i]);
+    Sln0.push_back(bl->U[iPoint*nVar+i]);
   } // End for
 
   // Initialize time integration variables
@@ -99,11 +99,12 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
 
     if(innerIters % itFrozenJac == 0)
     {
-      if (nErrors > 4)
+      if (nErrors >= 5)
       {
+        ResetSolution(iPoint, bl, Sln0);
         return -1; // Convergence failed
       }
-      CFL = ControlIntegration(iPoint, bl, temp_U, CFL);
+      CFL = ControlIntegration(iPoint, bl, Sln0, CFL);
       itFrozenJac = itFrozenJac0;
       JacMatrix = SysMatrix->ComputeJacobian(iPoint, bl, SysRes, 1e-8);
     }
@@ -114,11 +115,6 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
     {
       bl->U[iPoint*nVar+k] += slnIncr(k);
     }
-    if (itFrozenJac0 == 1)
-    {
-      bl->U[iPoint*nVar + 0] = std::max(bl->U[iPoint*nVar + 0], 1e-7);
-      bl->U[iPoint*nVar + 1] = std::max(bl->U[iPoint*nVar + 1], 1.0005);
-    }
 
     SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
     normSysRes = (SysRes + slnIncr/timeStep).norm();
@@ -142,7 +138,7 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
   {
     if (std::isnan(bl->U[iPoint*nVar + k]))
     {
-      ResetSolution(iPoint, bl, temp_U);
+      ResetSolution(iPoint, bl, Sln0);
       return innerIters; // New CFL.
     }
   }
@@ -156,7 +152,7 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
   }
    
 
-  ResetSolution(iPoint, bl, temp_U);
+  ResetSolution(iPoint, bl, Sln0);
   return innerIters;
 } // End Integration
 
@@ -178,46 +174,50 @@ double TimeSolver::ComputeSoundSpeed(double gradU2)
   return sqrt(1 / (Minf * Minf) + 0.5 * (gamma - 1) * (1 - gradU2));
 }
 
-double TimeSolver::ControlIntegration(size_t iPoint, BLRegion *bl, std::vector<double> temp_U, double CFLIn){
+double TimeSolver::ControlIntegration(size_t iPoint, BLRegion *bl, std::vector<double> Sln0, double CFLIn){
 
   /* Check if there are NaN in the solution */
   unsigned int nVar = bl->GetnVar();
+
+  /* Check for eventual errors */ 
+
+  bl->U[iPoint*nVar + 0] = std::max(bl->U[iPoint*nVar + 0], 1e-7);
+  bl->U[iPoint*nVar + 1] = std::max(bl->U[iPoint*nVar + 1], 1.0005);
+
   for (size_t k=0; k<nVar; ++k)
   {
     if (std::isnan(bl->U[iPoint*nVar + k]))
     {
-      ResetSolution(iPoint, bl, temp_U);
+      ResetSolution(iPoint, bl, Sln0, true);
       itFrozenJac0 = 1; // Impose continuous Jacobian evaluation.
       nErrors+=1;
-      return 0.3; // New CFL.
+      return 0.1; // New CFL.
     }
   }
   return CFLIn; // Keep CFL.
 }
 
-void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> temp_U){
+void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> Sln0, bool usePrevPoint)
+{
 
   unsigned int nVar = bl->GetnVar();
 
   /* Reset solution */
 
-  for(size_t k=0; k<nVar; ++k)
+  if (usePrevPoint)
   {
-    if (std::isnan(bl->U[iPoint*nVar+k]))
-    {
-      for (size_t j = 0; k<bl->GetnVar(); ++j)
-      {
-        if(!std::isnan(temp_U[j]))
-        {
-        bl->U[iPoint*bl->GetnVar()+j] = temp_U[j];
-        }
-        else
-        {
-           bl->U[iPoint*nVar+j] = bl->U[(iPoint-1)*nVar+j];
-        }
-      break;
-    } // End if
+    for(size_t k=0; k<nVar; ++k)
+  {
+    bl->U[iPoint*nVar + k] = bl->U[(iPoint-1)*nVar + k];
   } // End for
+  } // End if (usePrevPoint)
+  else
+  {
+    for(size_t k=0; k<nVar; ++k)
+    {
+      bl->U[iPoint*nVar + k] = Sln0[k];
+    } // End for
+  } // End else (usePrevPoint)
 
   /* Reset closures */
 
@@ -244,8 +244,7 @@ void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double>
     bl->Us[iPoint] = turbParam[7];
   } // End else if
 
-  } // End for
-}
+} // End ResetSolution
 
 
 double TimeSolver::GetCFL0(){return CFL0;}
diff --git a/dart/src/wTimeSolver.h b/dart/src/wTimeSolver.h
index e4d914d8733ea671291b6640b1b2e14248414cff..dcf8baebf38a1481161c49548622792e7b9e34ec 100644
--- a/dart/src/wTimeSolver.h
+++ b/dart/src/wTimeSolver.h
@@ -41,8 +41,8 @@ public:
 private:
   double SetTimeStep(double CFL, double dx, double invVel);
   double ComputeSoundSpeed(double gradU2);
-  double ControlIntegration(size_t iPoint, BLRegion *bl, std::vector<double> temp_U, double CFLIn);
-  void ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> temp_U);
+  double ControlIntegration(size_t iPoint, BLRegion *bl, std::vector<double> Sln0, double CFLIn);
+  void ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> Sln0, bool usePrevPoint=false);
 
 };
 } // namespace dart
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index bce4235b212fcd20640fa96ed78ed263cf7a00e3..f4b16a9061054fa71f437cbe7cfbb9e365f1518f 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -102,16 +102,19 @@ ViscSolver::~ViscSolver()
  */
 void ViscSolver::InitMeshes(size_t region, std::vector<double> locM, std::vector<double> globM)
 {
+  bl[region]->coarseSize = locM.size();
   if (mapMesh)
   {
     std::vector<double> fLocMesh = meshConverter->CreateFineMesh(locM);
     std::vector<double> fGlobMesh = meshConverter->CreateFineMesh(globM);
+    bl[region]->fineSize = fLocMesh.size();
     bl[region]->SetMesh(fLocMesh, fGlobMesh);
     std::cout << "- " << bl[region]->name << ": Mesh mapped from " << locM.size() << " to " << bl[region]->GetnMarkers() << " elements." << std::endl;
   }
   else
   {
     bl[region]->SetMesh(locM, globM);
+    bl[region]->fineSize = locM.size();
     std::cout << "- " << bl[region]->name << ": Mesh elements " << bl[region]->GetnMarkers() << std::endl;
   }
 }
@@ -126,9 +129,9 @@ void ViscSolver::SendInvBC(size_t region, std::vector<double> _Ue,
 
   if (mapMesh)
   {
-    std::vector<double> UeFine = meshConverter->InterpolateToFineMesh(_Ue);
-    std::vector<double> MeFine = meshConverter->InterpolateToFineMesh(_Me);
-    std::vector<double> RhoeFine = meshConverter->InterpolateToFineMesh(_Rhoe);
+    std::vector<double> UeFine = meshConverter->InterpolateToFineMesh(_Ue, bl[region]->fineSize);
+    std::vector<double> MeFine = meshConverter->InterpolateToFineMesh(_Me, bl[region]->fineSize);
+    std::vector<double> RhoeFine = meshConverter->InterpolateToFineMesh(_Rhoe, bl[region]->fineSize);
 
     bl[region]->SetInvBC(UeFine, MeFine, RhoeFine);
     std::cout << "- " << bl[region]->name << ": Inviscid boundary interpolated." << std::endl;
@@ -147,8 +150,8 @@ void ViscSolver::SendExtVariables(size_t region, std::vector<double> _xxExt, std
 {
   if (mapMesh)
   {
-    std::vector<double> xxExtFine = meshConverter->InterpolateToFineMesh(_xxExt);
-    std::vector<double> deltaStarExtFine = meshConverter->InterpolateToFineMesh(_deltaStarExt);
+    std::vector<double> xxExtFine = meshConverter->InterpolateToFineMesh(_xxExt, bl[region]->fineSize);
+    std::vector<double> deltaStarExtFine = meshConverter->InterpolateToFineMesh(_deltaStarExt, bl[region]->fineSize);
     bl[region]->SetExtVariables(xxExtFine, deltaStarExtFine);
   }
   else
@@ -180,7 +183,7 @@ int ViscSolver::Run(unsigned int couplIter){
   }
 
   else{
-    tSolver->SetCFL0(0.5);
+    tSolver->SetCFL0(CFL0);
     tSolver->SetitFrozenJac(5);
     tSolver->SetinitSln(false);
     stagPtMvmt = false;
@@ -370,7 +373,16 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
  */
 std::vector<double> ViscSolver::ExtractBlowingVel(size_t iRegion)
 {
-  return bl[iRegion]->Blowing;
+  std::vector<double> cBlwVel;
+  if (mapMesh)
+  { 
+    cBlwVel = meshConverter->BlowingFineToCoarse(bl[iRegion]->Blowing, bl[iRegion]->coarseSize-1);
+  }
+  else
+  {
+    cBlwVel = bl[iRegion]->Blowing;
+  }
+  return cBlwVel;
 }
 
 /**
@@ -378,7 +390,17 @@ std::vector<double> ViscSolver::ExtractBlowingVel(size_t iRegion)
  */
 std::vector<double> ViscSolver::ExtractDeltaStar(size_t iRegion)
 {
-  return bl[iRegion]->DeltaStar;
+  std::vector<double> cDStar;
+  if (mapMesh)
+  { 
+    std::cout << "DeltaStar size" << bl[iRegion]->DeltaStar.size() << std::endl;
+    cDStar = meshConverter->VarsFineToCoarse(bl[iRegion]->DeltaStar, bl[iRegion]->coarseSize);
+  }
+  else
+  {
+    cDStar = bl[iRegion]->DeltaStar;
+  }
+  return cDStar;
 }
 
 /**
@@ -393,7 +415,16 @@ std::vector<double> ViscSolver::ExtractXx(size_t iRegion)
   for (size_t k=0; k<bl[iRegion]->GetnMarkers(); ++k){
     LocMarkersOut.push_back(bl[iRegion]->GetLocMarkers(k));
   }
-  return LocMarkersOut;
+  std::vector<double> cXxExt;
+  if (mapMesh)
+  { 
+    cXxExt = meshConverter->VarsFineToCoarse(LocMarkersOut, bl[iRegion]->coarseSize);
+  }
+  else
+  {
+    cXxExt = LocMarkersOut;
+  }
+  return cXxExt;
 }
 
 /**
@@ -446,7 +477,7 @@ void ViscSolver::ComputeDrag()
 std::vector<std::vector<double>> ViscSolver::ExtractSolution()
 {
   // [xx, x, deltaStar, theta, H, N, Ue, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta, VtExt, Me, Rhoe]
-  std::vector<std::vector<double>> blVectors(17);
+  std::vector<std::vector<double>> blVectors(18);
   unsigned int nVar = bl[0]->GetnVar();
 
   for(size_t iRegion=0; iRegion<bl.size(); ++iRegion)
@@ -470,6 +501,7 @@ std::vector<std::vector<double>> ViscSolver::ExtractSolution()
       blVectors[14].push_back(bl[iRegion]->GetVtExt(iPoint)); // VtExt
       blVectors[15].push_back(bl[iRegion]->GetMe(iPoint)); // Me
       blVectors[16].push_back(bl[iRegion]->GetRhoe(iPoint)); // Rhoe
+      blVectors[17].push_back(bl[iRegion]->Blowing[iPoint]); // Blowing velocity
     }
   }
   return blVectors;
diff --git a/dart/tests/bliNonLift.py b/dart/tests/bliNonLift.py
index 24c59b0b175398c2cf3b67b13faf676295bdc0ae..27b26d298646dc42ca10bf4d11fcfe9be557a34e 100755
--- a/dart/tests/bliNonLift.py
+++ b/dart/tests/bliNonLift.py
@@ -60,7 +60,7 @@ def main():
 
     # define flow variables
     Re = 1e7
-    alpha = 2.*math.pi/180
+    alpha = 0.*math.pi/180
     M_inf = 0.
 
     #user_xtr=[0.41,0.41]
@@ -68,7 +68,7 @@ def main():
     user_xtr=[None,None]
     user_Ti=None
 
-    mapMesh = 0
+    mapMesh = 1
     nFactor = 2
 
     # Time solver parameters
@@ -86,7 +86,7 @@ def main():
     # mesh the geometry
     print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
     tms['msh'].start()
-    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
     #pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.002, 'msLe' : 0.0001}
     msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
diff --git a/dart/tests/blicpp.py b/dart/tests/blicpp.py
index 0046b19aaae268d2fa76610b12c312b66da85145..3437e134e9f0719faaacb757a7f18c888c7bd3fd 100644
--- a/dart/tests/blicpp.py
+++ b/dart/tests/blicpp.py
@@ -56,10 +56,10 @@ def main():
 
     # define flow variables
     Re = 1e7
-    alpha = 2*math.pi/180
+    alpha = 15*math.pi/180
     M_inf = 0.2
 
-    meshFactor = 1
+    meshFactor = 2
     CFL0 = 1
 
     plotVar = 'Hk'
@@ -72,7 +72,7 @@ def main():
     # mesh the geometry
     print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
     tms['msh'].start()
-    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.0001}
     msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
diff --git a/dart/vCoupler.py b/dart/vCoupler.py
index 0884c4cda7679d426998c48f022d663160369d62..4d218f2b211e425cfd4de04d3024b10681f6844e 100644
--- a/dart/vCoupler.py
+++ b/dart/vCoupler.py
@@ -26,6 +26,7 @@ from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
 from fwk.coloring import ccolors
 
 import fwk
+import dart.viscUtils as viscU
 import numpy as np
 from matplotlib import pyplot as plt
 
@@ -110,6 +111,13 @@ class Coupler:
         print('Viscous ops terminated with exit code ', exitCode)
 
         error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
+        
+        """# Plot State.
+        plotVar = 'H'
+        blScalars, blVectors = viscU.ExtractVars(self.vSolver)
+        xtr=self.vSolver.Getxtr()
+        wData = viscU.Dictionary(blScalars, blVectors, xtr)
+        viscU.PlotVars(wData, plotVar)"""
 
         if error  <= self.couplTol and exitCode==0:
 
@@ -134,13 +142,41 @@ class Coupler:
           self.deltaStar[iRegion] = np.asarray(self.vSolver.ExtractDeltaStar(iRegion))
           self.xx[iRegion] = np.asarray(self.vSolver.ExtractXx(iRegion))
 
-        blwVelAF = np.concatenate((self.blwVel[0], self.blwVel[1][1:])) # Remove stagnation point doublon.
+        """plt.plot(self.blwVel[0], 'o')
+        plt.plot(self.blwVel[1])
+        plt.show()"""
+
+        blwVelAF = np.concatenate((self.blwVel[0], self.blwVel[1]))
         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]
 
+        """plotVar = 'H'
+        blScalars, blVectors = viscU.ExtractVars(self.vSolver)
+        xtr=self.vSolver.Getxtr()
+        wData = viscU.Dictionary(blScalars, blVectors, xtr)
+        for i in range(len(wData['x/c'])):
+          if wData['x/c'][i] == 1.0:
+            nUpper = i+1
+            break
+
+        plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], self.deltaStar[0], 'o')
+        plt.plot(wData['x/c'][:nUpper], wData['deltaStar'][:nUpper], 'x')
+        plt.legend(['Coarse', 'Fine'], frameon=False)
+        plt.show()
+
+        plt.plot(np.linspace(0,1, num=len(self.blwVel[0])), self.blwVel[0], 'o')
+        plt.plot(np.linspace(0,1,len(wData['Blowing'][:nUpper])), wData['Blowing'][:nUpper], 'x')
+        plt.legend(['Coarse', 'Fine'], frameon=False)
+        plt.show()
+
+        plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], self.xx[0], 'o')
+        plt.plot(wData['x/c'][:nUpper], wData['xx'][:nUpper], 'x')
+        plt.legend(['Coarse', 'Fine'], frameon=False)
+        plt.show()"""
+
         self.group[0].u = blwVelAF[self.group[0].connectListElems.argsort()]
         self.group[1].u = blwVelWK[self.group[1].connectListElems.argsort()]
 
diff --git a/dart/viscUtils.py b/dart/viscUtils.py
index 7cb75d80d32bcc20beaf49ca660da4ff3de1e227..45355a31391013e8e022f3f44a01bfc6d4896d22 100644
--- a/dart/viscUtils.py
+++ b/dart/viscUtils.py
@@ -24,6 +24,7 @@ def Dictionary(blScalars, blVectors, xtr):
   wData['VtExt']      = blVectors[14]
   wData['Me']         = blVectors[15]
   wData['Rhoe']       = blVectors[16]
+  wData['Blowing']    = blVectors[17]
   wData['xtrTop']     = xtr[0]
   wData['xtrBot']     = xtr[1]
   return wData
@@ -62,7 +63,7 @@ def ExtractVars(viscSolver):
 
   blScalars = [viscSolver.Cdt, viscSolver.Cdf, viscSolver.Cdp]
   
-  # [xx, x, deltaStar, theta, H, N, Ue, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta, VtExt, Me, Rhoe]
+  # [xx, x, deltaStar, theta, H, N, Ue, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta, VtExt, Me, Rhoe, Blowing]
 
   blVectors = viscSolver.ExtractSolution()
 
@@ -87,6 +88,6 @@ def PlotVars(wData, var='H'):
   plt.axvline(x=wData['xtrBot'], linestyle="--", color='#D95319')
   plt.xlabel('$x/c$')
   plt.ylabel(var)
-  plt.xlim([0,1])
+  plt.xlim([0,2])
   plt.legend(['Upper side', 'Lower side'], frameon=False)
   plt.show()
\ No newline at end of file