diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index 7de8cbdd2721cdca4f3389dc9774525c1564dfb8..d9c89a69c04e4a1d6fbc2e2e5f5010e618f5a974 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -27,7 +27,6 @@ void BLRegion::SetMesh(std::vector<double> locM, std::vector<double> globM)
   LocMarkers=locM;
   GlobMarkers=globM;
   nMarkers = locM.size();
-  std::cout << "nMarkers" << nMarkers << std::endl;
 
   /* Resize all variables */
 
@@ -84,11 +83,15 @@ void BLRegion::SetStagBC(double Re)
   
   Ret[0] = U[3]*U[0]*Re;
 
-  if (Ret[0] < 1){
-    Ret[0] = 0;
+  if (Ret[0] < 1)
+  {
+    Ret[0] = 1;
     U[3] = Ret[0]/(U[0]*Re);
   } // End if
 
+  //std::cout << "name " << name << " dv0 " << dv0 << " U3 " << U[3] << " U[0] " << U[0] << " \n" << std::scientific;
+
+
   DeltaStar[0] = U[0] * U[1];
 
   std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], Me[0], name);
@@ -101,6 +104,7 @@ void BLRegion::SetStagBC(double Re)
   Delta[0] = ClosParam[5];
   Cteq[0] = ClosParam[6];
   Us[0] = ClosParam[7];
+
 } // End SetStagBC
 
 void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<double> LowerCond)
@@ -133,9 +137,11 @@ void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<d
   Cteq[0] = ClosParam[6];
   Us[0] = ClosParam[7];
 
-  for (size_t k=0; k<nMarkers; ++k){
+  for (size_t k=0; k<nMarkers; ++k)
+  {
     Regime[k] = 1;
   }
+  PrintSolution(0);
 }
 
 void BLRegion::ComputeBlwVel()
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index 5378e868958f8fc1efe9ed70de5c9e0180c5490a..77ea0594c9f323cdcec41204f0d0ef201a223559 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -72,11 +72,6 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
   for (size_t i=0; i<nVar; ++i){
     temp_U.push_back(bl->U[iPoint*nVar+i]);
   } // End for
-  /*if (bl->name=="upper" && iPoint==67){
-    for (size_t k=0; k<nVar; ++k){
-      std::cout << "k" << temp_U[k];
-    }
-  }*/
 
   // Initialize time integration variables
 
@@ -97,50 +92,77 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
   Matrix<double, 5, 5> JacMatrix;
   Vector<double, 5> slnIncr;
   
-  unsigned int nErrors = 0; // Number of errors encountered
+  nErrors = 0; // Number of errors encountered
   unsigned int innerIters = 0; // Inner (non-linear) iterations
 
   while (innerIters < maxIt){
 
-
-    /*if (bl->name=="upper" && iPoint==67){
-      bl->PrintSolution(iPoint);
-      std::cout << "CFL " << CFL << std::endl;
-    }*/
-    if(innerIters % itFrozenJac == 0){
+    if(innerIters % itFrozenJac == 0)
+    {
+      if (nErrors > 4)
+      {
+        return -1; // Convergence failed
+      }
+      CFL = ControlIntegration(iPoint, bl, temp_U, CFL);
+      itFrozenJac = itFrozenJac0;
       JacMatrix = SysMatrix->ComputeJacobian(iPoint, bl, SysRes, 1e-8);
     }
 
-    // Ici on peut faire autrement : .asDiagonal()
     slnIncr = (JacMatrix + IMatrix).partialPivLu().solve(-SysRes);
 
-    for(size_t k=0; k<nVar; ++k){
+    for(size_t k=0; k<nVar; ++k)
+    {
       bl->U[iPoint*nVar+k] += slnIncr(k);
-      //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);
     }
+    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();
+    if (bl->name == "wake" && iPoint == 20)
+    {
+      std::cout << normSysRes/normSysRes0 << std::scientific;
+      bl->PrintSolution(iPoint);
+    }
 
-    if (normSysRes <= tol * normSysRes0){
-      /* Successfull exit */
-      return 0;
+    if (normSysRes <= tol * normSysRes0)
+    {
+      return 0; /* Successfull exit */
     }
 
     // CFL adaptation.
 
     CFL = CFL0* pow(normSysRes0/normSysRes, 0.7);
-
     CFL = std::max(CFL, 1.0);
     timeStep = SetTimeStep(CFL, dx, bl->GetVtExt(iPoint));
     IMatrix = Matrix<double, 5, 5>::Identity() / timeStep;
 
     innerIters++;
   } // End while (innerIters < maxIt)
-  
+
+  for (size_t k=0; k<nVar; ++k)
+  {
+    if (std::isnan(bl->U[iPoint*nVar + k]))
+    {
+      ResetSolution(iPoint, bl, temp_U);
+      return innerIters; // New CFL.
+    }
+  }
+  if (normSysRes <= 1e-4 * normSysRes0)
+  {
+    return 0;
+  }
+  else
+  {
+    return innerIters;
+  }
+   
+
   ResetSolution(iPoint, bl, temp_U);
-  bl->PrintSolution(iPoint);
   return innerIters;
 } // End Integration
 
@@ -162,20 +184,43 @@ 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){
+
+  /* Check if there are NaN in the solution */
+  unsigned int nVar = bl->GetnVar();
+  for (size_t k=0; k<nVar; ++k)
+  {
+    if (std::isnan(bl->U[iPoint*nVar + k]))
+    {
+      //std::cout << "Point " << iPoint << " : Solution reset; NaN identified at position " << k << std::endl;
+      ResetSolution(iPoint, bl, temp_U);
+      itFrozenJac0 = 1; // Impose continuous Jacobian evaluation.
+      nErrors+=1;
+      return 0.3; // New CFL.
+    }
+  }
+  return CFLIn; // Keep CFL.
+}
+
 void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> temp_U){
 
-  std::cout << "Reseting point" << iPoint << std::endl;
+  //std::cout << "Reseting point" << iPoint << std::endl;
   unsigned int nVar = bl->GetnVar();
 
   /* Reset solution */
 
-  for(size_t k=0; k<nVar; ++k){
-    if (std::isnan(bl->U[iPoint*nVar+k])){
-      for (size_t j = 0; k<bl->GetnVar(); ++j){
-        if(!std::isnan(temp_U[j])){
+  for(size_t k=0; k<nVar; ++k)
+  {
+    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{
+        else
+        {
            bl->U[iPoint*nVar+j] = bl->U[(iPoint-1)*nVar+j];
         }
       break;
diff --git a/dart/src/wTimeSolver.h b/dart/src/wTimeSolver.h
index 07c160dc85184c9da2f4aa6e56ca2ee4a58ef729..e4d914d8733ea671291b6640b1b2e14248414cff 100644
--- a/dart/src/wTimeSolver.h
+++ b/dart/src/wTimeSolver.h
@@ -21,6 +21,7 @@ private:
   bool initSln;
   double gamma = 1.4;
   double Minf;
+  unsigned int nErrors;
 
   ViscFluxes *SysMatrix;
 
@@ -40,6 +41,7 @@ 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);
 
 };
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index b99a85ece9388de047c40bacfef4d4748b859134..83cd1a78b97b418758552ba0c14c7cebf45e19d1 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -134,7 +134,6 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
 
   double Mea = bl->GetMe(iPoint);
 
-  double Reta = bl->Ret[iPoint];
   double Hstara = bl->Hstar[iPoint];
   double Hstar2a = bl->Hstar2[iPoint];
   double Hka = bl->Hk[iPoint];
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index dc2ec99107de7d8c4d8808f73947bffefba34ce5..9d0e577eeda60896b732c8c2c9ed2151b28cc9a9 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -82,8 +82,6 @@ int ViscSolver::Run(unsigned int couplIter){
 
   /* Prepare time integration */
 
-  //std::cout << "couplIter" << couplIter << "nbElmUpper: " << nbElmUpper << "bl[0]->GetnMarkers()" << bl[0]->GetnMarkers() << "stagPtMvmt" << stagPtMvmt << std::endl;
-
   if ((couplIter == 0) || (nbElmUpper != bl[0]->GetnMarkers()) || (stagPtMvmt)){
     tSolver->SetCFL0(0.5);
     tSolver->SetitFrozenJac(1);
@@ -105,6 +103,7 @@ int ViscSolver::Run(unsigned int couplIter){
     stagPtMvmt = false;
     std::cout << "Updating solution" << std::endl;
   }
+
   nbElmUpper = bl[0]->GetnMarkers();
 
   /* Set boundary condition */
@@ -115,6 +114,8 @@ int ViscSolver::Run(unsigned int couplIter){
   /* Loop over the regions */
 
   bool lockTrans;
+  int solverExitCode = 0;
+  int pointExitCode;
 
   for (size_t iRegion = 0; iRegion < bl.size(); ++iRegion){
 
@@ -155,11 +156,17 @@ int ViscSolver::Run(unsigned int couplIter){
 
       /* Time integration */
 
-      tSolver->Integration(iPoint, bl[iRegion]);
+      pointExitCode = tSolver->Integration(iPoint, bl[iRegion]);
+      if (pointExitCode!=0)
+      {
+        std::cout << "Point " << iPoint << ": Convergence terminated with exitcode " << pointExitCode << std::endl;
+        solverExitCode = -1; // Convergence failed
+      }
 
       /* Transition */
 
-      if (!lockTrans){
+      if (!lockTrans)
+      {
         CheckWaves(iPoint, bl[iRegion]);
 
         // Free transition.
@@ -175,12 +182,13 @@ int ViscSolver::Run(unsigned int couplIter){
 
   ComputeDrag();
 
-  for (size_t iRegion=0; iRegion<bl.size(); ++iRegion){
+  for (size_t iRegion=0; iRegion<bl.size(); ++iRegion)
+  {
     bl[iRegion]->ComputeBlwVel();
   } // End for iRegion
 
-  return 0; // Successfull exit
-}
+  return solverExitCode; // exit code
+  }
 
 void ViscSolver::CheckWaves(size_t iPoint, BLRegion *bl)
 {
@@ -189,7 +197,8 @@ void ViscSolver::CheckWaves(size_t iPoint, BLRegion *bl)
   double logRet_crit = 2.492*pow(1/(bl->Hk[iPoint]-1), 0.43)
   + 0.7*(tanh(14*(1/(bl->Hk[iPoint]-1))-9.24)+1.0);
 
-  if (bl->Ret[iPoint] > 0){ // Avoid log of negative number.
+  if (bl->Ret[iPoint] > 0)
+  { // Avoid log of negative number.
     if (log10(bl->Ret[iPoint]) <= logRet_crit - 0.08){
       bl->U[iPoint*bl->GetnVar() + 2] = 0; // Set N to 0 at that point if waves do not grow.
     }
@@ -204,7 +213,8 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
 
   // Save laminar solution.
   std::vector<double> lamSol;
-  for (size_t k=0; k<nVar; ++k){
+  for (size_t k=0; k<nVar; ++k)
+  {
     lamSol.push_back(bl->U[iPoint*nVar + k]);
   }
 
@@ -212,7 +222,8 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
   bl->transMarker = iPoint; /* Save transition marker */
 
   /* Loop over remaining points in the region */
-  for (size_t k=0; k<bl->GetnMarkers() - iPoint; ++k){
+  for (size_t k=0; k<bl->GetnMarkers() - iPoint; ++k)
+  {
     bl->Regime[iPoint+k] = 1; // Set Turbulent regime. 
   }
 
@@ -288,31 +299,31 @@ void ViscSolver::ComputeDrag()
 
   unsigned int nVar = bl[0]->GetnVar();
   unsigned int lastWkPt = (bl[2]->GetnMarkers() - 1)*nVar;
+
   /* Normalize friction and dissipation coefficients. */
 
-  /* std::vector<double> frictioncoeff_upper;
+  std::vector<double> frictioncoeff_upper;
   for (size_t k=0; k<bl[0]->GetnMarkers(); ++k){
-    frictioncoeff_upper.push_back(bl[0]->U[k*nVar + 3] * bl[0]->Cf[k]);
+    frictioncoeff_upper.push_back(bl[0]->U[k*nVar + 3]*bl[0]->U[k*nVar + 3] * bl[0]->Cf[k]);
   }
   std::vector<double> frictioncoeff_lower;
   for (size_t k=0; k<bl[1]->GetnMarkers(); ++k){
-    frictioncoeff_upper.push_back(bl[1]->U[k*nVar + 3] * bl[1]->Cf[k]);
+    frictioncoeff_lower.push_back(bl[1]->U[k*nVar + 3]*bl[1]->U[k*nVar + 3] * bl[1]->Cf[k]);
   }
 
   /* Compute friction drag coefficient. */
 
-  /* double Cdf_upper = 0.0;
+  double Cdf_upper = 0.0;
   for (size_t k=1; k<bl[0]->GetnMarkers(); ++k){
-    Cdf_upper += (bl[0]->GetGlobMarkers(k) - bl[0]->GetGlobMarkers(k-1))*frictioncoeff_upper[k-1];
+    Cdf_upper += (frictioncoeff_upper[k]+frictioncoeff_upper[k-1])*(bl[0]->GetGlobMarkers(k) - bl[0]->GetGlobMarkers(k-1))/2;
   }
   double Cdf_lower = 0.0;
   for (size_t k=1; k<bl[1]->GetnMarkers(); ++k){
-    Cdf_upper += (bl[1]->GetGlobMarkers(k) - bl[1]->GetGlobMarkers(k-1))*frictioncoeff_lower[k-1];
-  }
+    Cdf_lower += (frictioncoeff_lower[k]+frictioncoeff_lower[k-1])*(bl[1]->GetGlobMarkers(k) - bl[1]->GetGlobMarkers(k-1))/2;
 
-  */
+  }
 
-  // Cdf = Cdf_upper + Cdf_lower; // No friction in the wake
+  Cdf = Cdf_upper + Cdf_lower; // No friction in the wake
   
   /* Compute total drag coefficient. */
 
@@ -320,5 +331,5 @@ void ViscSolver::ComputeDrag()
   
   /* Compute pressure drag coefficient. */
 
-  //Cdp = Cdt - Cdf;
+  Cdp = Cdt - Cdf;
 }
\ No newline at end of file
diff --git a/dart/tests/blicpp.py b/dart/tests/blicpp.py
index 76f57dd1b14c2613015a6f733d0c983a45f618e2..e5c52372fbcf8fbf97a7ea24fe2e81f28de53b35 100644
--- a/dart/tests/blicpp.py
+++ b/dart/tests/blicpp.py
@@ -56,7 +56,7 @@ def main():
     # define flow variables
     Re = 1e7
     alpha = 2*math.pi/180
-    M_inf = 0.715
+    M_inf = 0.
 
     #user_xtr=[0.41,0.41]
     #user_xtr=[0,None]
@@ -81,9 +81,9 @@ 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.005}
+    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/rae2822.geo', pars, ['field', 'airfoil', 'downstream'])
+    msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
     # set the problem and add medium, initial/boundary conditions
diff --git a/dart/vCoupler.py b/dart/vCoupler.py
index d5e0fa8914cdf897b3b9fc452374f66260912046..21f29243ac35f7ffafea0dc08b832d61d6899613 100644
--- a/dart/vCoupler.py
+++ b/dart/vCoupler.py
@@ -48,7 +48,8 @@ class Coupler:
     def Run(self):
 
       # Initilize meshes in c++ objects
-      
+
+      nElmUpperPrev = 0
       couplIter = 0
       cdPrev = 0
       while couplIter < self.maxCouplIter:
@@ -73,20 +74,25 @@ 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)
-        if couplIter == 0:
+
+        if ((len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]])) != (nElmUpperPrev or couplIter==0)):
+          
           # Initialize mesh          
           self.vSolver.InitMeshes(0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]], fMeshDict['globCoord'][:fMeshDict['bNodes'][1]])
           self.vSolver.InitMeshes(1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
           self.vSolver.InitMeshes(2, fMeshDict['locCoord'][fMeshDict['bNodes'][2]:], fMeshDict['globCoord'][fMeshDict['bNodes'][2]:])
 
           self.vSolver.SendExtVariables(0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]], fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]])
-          self.vSolver.SendExtVariables(1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
+          self.vSolver.SendExtVariables(1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])    
           self.vSolver.SendExtVariables(2, fMeshDict['locCoord'][fMeshDict['bNodes'][2]:], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:])
+       
         else:
+
           self.vSolver.SendExtVariables(0, fMeshDict['xxExt'][:fMeshDict['bNodes'][1]], fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]])
           self.vSolver.SendExtVariables(1, fMeshDict['xxExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
           self.vSolver.SendExtVariables(2, fMeshDict['xxExt'][fMeshDict['bNodes'][2]:], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:])
 
+        nElmUpperPrev = len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]])
         self.vSolver.SendInvBC(0, fMeshDict['vtExt'][:fMeshDict['bNodes'][1]], fMeshDict['Me'][:fMeshDict['bNodes'][1]], fMeshDict['rho'][:fMeshDict['bNodes'][1]])
         self.vSolver.SendInvBC(1, fMeshDict['vtExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['Me'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['rho'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
         self.vSolver.SendInvBC(2, fMeshDict['vtExt'][fMeshDict['bNodes'][2]:], fMeshDict['Me'][fMeshDict['bNodes'][2]:], fMeshDict['rho'][fMeshDict['bNodes'][2]:])
@@ -98,7 +104,7 @@ class Coupler:
 
         error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
 
-        if error  <= self.couplTol:
+        if error  <= self.couplTol and exitCode==0:
           print('')
           print(ccolors.ANSI_GREEN, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} < {:<2.3f}'
           .format(couplIter, self.vSolver.Cdt, np.log10(error), np.log10(self.couplTol)), ccolors.ANSI_RESET)