diff --git a/vii/pyVII/vCoupler.py b/vii/pyVII/vCoupler.py
index bab595e64faea57f833456b46f0f20ddc4566e11..5268b75b2838f3465bf4434511dc0e173fb11992 100644
--- a/vii/pyVII/vCoupler.py
+++ b/vii/pyVII/vCoupler.py
@@ -171,9 +171,6 @@ class Coupler:
     self.group[0].xx = xxAF[self.group[0].connectListNodes.argsort()]
     self.group[1].xx = xxWK[self.group[1].connectListNodes.argsort()]
 
-    plt.plot(np.linspace(0,2, len(self.group[0].u)),self.group[0].u)
-    plt.show()
-
     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])
diff --git a/vii/pyVII/vCoupler2.py b/vii/pyVII/vCoupler2.py
index 05090fbb883b566df433fcd65eb7a1988971dcc3..7e9aad362e8e80e6cdfb51da35384b456383d712 100644
--- a/vii/pyVII/vCoupler2.py
+++ b/vii/pyVII/vCoupler2.py
@@ -47,10 +47,11 @@ class Coupler:
     self.nElmUpperPrev = 0
     couplIter = 0
     cdPrev = 0
-    print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'
-    .format(self.iSolverAPI.GetAlpha()*180/M_PI, self.iSolverAPI.GetMinf(), self.vSolver.GetRe()/1e6))
-    print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
-    print(' --------------------------------------------------------')
+    if self.iSolverAPI.iSolver.verbose == 0 and self.vSolver.verbose == 0:
+      print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'
+      .format(self.iSolverAPI.GetAlpha()*180/M_PI, self.iSolverAPI.GetMinf(), self.vSolver.GetRe()/1e6))
+      print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
+      print(' --------------------------------------------------------')
     while couplIter < self.maxCouplIter:
 
       # Run inviscid solver
@@ -88,7 +89,6 @@ class Coupler:
       
       cdPrev = self.vSolver.Cdt
 
-
       print('- {:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.3f} {:>6.3f} {:>8.2f}'.format(couplIter,
       self.iSolverAPI.GetCl(), self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)))
 
diff --git a/vii/pyVII/vInterpolator.py b/vii/pyVII/vInterpolator.py
index 5c852d3c484886e9e9753865ba3a75dfc453c07a..c252a291ccff38040de5ca5073c668aba16ae9ae 100644
--- a/vii/pyVII/vInterpolator.py
+++ b/vii/pyVII/vInterpolator.py
@@ -57,8 +57,6 @@ class Interpolator:
                 dummy[:,1] += sym
                 self.idata[iReg] = np.row_stack((self.idata[iReg], dummy))
         self.vdata, self.vcg = self.GetMesh(_vmsh, visc=1)
-        plt.plot(self.vdata[0][:,0], self.vdata[0][:,1],'o')
-        plt.show()
         self.Xp, self.Yp, self.Zp, self.sortedRows, self.Xc, self.Yc, self.Zc, self.sortedElems = self.SortMesh(self.vdata, self.vcg)
         self.tms['Mesh sort'].stop()
         print(self.tms)
@@ -74,9 +72,9 @@ class Interpolator:
         ax.set_zlabel('z')
         #plt.gca().invert_yaxis()
         plt.axis('off')
-        plt.show()"""
+        plt.show()
 
-        """fig = plt.figure()
+        fig = plt.figure()
         ax = fig.add_subplot(projection='3d')
         ax.scatter(self.idata[0][:,0], self.idata[0][:,1], self.idata[0][:,2], s=0.3)
         ax.scatter(self.icg[0][:,0], self.icg[0][:,1], self.icg[0][:,2], s=0.3)
@@ -128,6 +126,8 @@ class Interpolator:
                 U[iReg][iNode,:] = [self.iSolver.U[row][0], self.iSolver.U[row][1], self.iSolver.U[row][2]]
                 M[iReg][iNode] = self.iSolver.M[row]
                 Rho[iReg][iNode] = self.iSolver.rho[row]
+        
+        #self.URani = U
 
         nD = self.cfg['nDim']
         Ux = [self.__RbfToSections(self.idata[0][:,:nD], U[0][:,0], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, 0:2], U[1][:,0], self.vdata[1][:,0:2])]
@@ -135,6 +135,8 @@ class Interpolator:
         Uz = [self.__RbfToSections(self.idata[0][:,:nD], U[0][:,2], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, 0:2], U[1][:,2], self.vdata[1][:,0:2])]
         Me = [self.__RbfToSections(self.idata[0][:,:nD], M[0], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, 0:2], M[1], self.vdata[1][:,0:2])]
         Rhoe = [self.__RbfToSections(self.idata[0][:,:nD], Rho[0], self.vdata[0][:,:nD]), self.__RbfToSections(self.idata[1][:, 0:2], Rho[1], self.vdata[1][:,0:2])]
+        
+        #self.URani = Ux[0]
         for iSec, y in enumerate(self.cfg['Sections']):
 
             for iReg in range(2):
@@ -236,7 +238,7 @@ class Interpolator:
                     blw[iSec][iReg] = np.concatenate((blw[iSec][iReg], blwUp[::-1], blwLw))
                 elif iReg == 1:
                     blw[iSec][iReg] = np.concatenate((blw[iSec][iReg], vSolver.ExtractBlowingVel(iSec, 2)))
-                #blw[iSec][iReg] = np.clip(blw[iSec][iReg], -0.01, 0.01)
+        
         for iReg in range(2):
 
             # Find index of all Sections in Y
@@ -322,6 +324,8 @@ class Interpolator:
 
 
             self.tms['Interpolation'].stop()
+            if iReg == 0:
+                self.URani = mappedBlw
             if iReg==0 or iReg==1:
                 for iElm, e in enumerate(self.blw[iReg].tag.elems):
                     if np.size(mappedBlw[self.icg[iReg][:,3]==e.no])>1:
diff --git a/vii/pyVII/vUtils.py b/vii/pyVII/vUtils.py
index e34467191158f649b7a0bb21097c2ca9dbfb19f7..b10033feeda92008388a1fcfb1f58d18fa9255e4 100644
--- a/vii/pyVII/vUtils.py
+++ b/vii/pyVII/vUtils.py
@@ -2,7 +2,7 @@ import vii
 from fwk.wutils import parseargs
 from fwk.coloring import ccolors
 
-def initBL(Re, Minf, CFL0, nSections, span=0, verbose=1):
+def initBL(Re, Minf, CFL0, nSections, span=0, _verbose=1):
     if Re<=0.:
         print(ccolors.ANSI_RED, "dart::vUtils Error : Negative Reynolds number.", Re, ccolors.ANSI_RESET)
         raise RuntimeError("Invalid parameter")
@@ -14,12 +14,12 @@ def initBL(Re, Minf, CFL0, nSections, span=0, verbose=1):
     if nSections < 0:
         print(ccolors.ANSI_RED, "dart::vUtils Fatal error : Negative number of sections.", nSections, ccolors.ANSI_RESET)
         raise RuntimeError("Invalid parameter")
-    if not(0<=verbose<=3):
-        print(ccolors.ANSI_RED, "dart::vUtils Fatal error : verbose not in valid range.", verbose, ccolors.ANSI_RESET)
+    if not(0<=_verbose<=3):
+        print(ccolors.ANSI_RED, "dart::vUtils Fatal error : verbose not in valid range.", _verbose, ccolors.ANSI_RESET)
         raise RuntimeError("Invalid parameter")
 
     args = parseargs()
-    solver = vii.ViscSolver(Re, Minf, CFL0, nSections, span, verbose=args.verb+1)
+    solver = vii.ViscSolver(Re, Minf, CFL0, nSections, span, verbose=_verbose)
     return solver
 
 
diff --git a/vii/src/wViscSolver.cpp b/vii/src/wViscSolver.cpp
index a89e87a35b8de41d7492970d170af6bf89fc978f..adecfc113dabd8f0111243fb303ab328744bfc43 100644
--- a/vii/src/wViscSolver.cpp
+++ b/vii/src/wViscSolver.cpp
@@ -15,6 +15,11 @@
 #include "wTimeSolver.h"
 #include <iomanip>
 
+#define ANSI_COLOR_RED "\x1b[1;31m"
+#define ANSI_COLOR_GREEN "\x1b[1;32m"
+#define ANSI_COLOR_YELLOW "\x1b[1;33m"
+#define ANSI_COLOR_RESET "\x1b[0m"
+
 /* --- Namespaces -------------------------------------------------------------------- */
 
 using namespace vii;
@@ -29,15 +34,13 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec
 
     PrintBanner();
 
-    /* Freestream parameters */
-
+    // Freestream parameters.
     Re = _Re;
     CFL0 = _CFL0;
 
     verbose = _verbose;
 
-    /* Initialzing boundary layer */
-
+    // Initialzing boundary layer
     Sections.resize(nSections);
     for (size_t iSec=0; iSec < nSections; ++iSec)
     {
@@ -57,26 +60,21 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec
     Cdt = 0.; Cdf = 0.; Cdp = 0.;
     Span = _Span;
 
-    /* Pre/post processing flags */
-
+    // Pre/post processing flags
     solverExitCode = 0;
     stagPtMvmt.resize(Sections.size(), false);
 
-    /* Temporal solver initialization */
-
+    // Temporal solver initialization.
     tSolver = new TimeSolver(_CFL0, _Minf, Re);
 
-    /* User input numerical parameters */
-
+    // User input numerical parameters.
     std::cout << "--- Viscous solver setup ---\n"
               << "Reynolds number: " << Re << " \n"
               << "Freestream Mach number: " << _Minf << " \n"
               << "Initial CFL number: " << CFL0 << " \n"
               << std::endl;
-
     std::cout << "Boundary layer initialized" << std::endl;
-
-} // end ViscSolver
+}
 
 /**
  * @brief Viscous solver and subsequent dependancies destruction
@@ -90,7 +88,6 @@ ViscSolver::~ViscSolver()
             delete Sections[iSec][i];
         }
     }
-
     delete tSolver;
     std::cout << "~ViscSolver()\n";
 } // end ~ViscSolver
@@ -101,30 +98,28 @@ ViscSolver::~ViscSolver()
  */
 int ViscSolver::Run(unsigned int const couplIter)
 {   
-    bool lockTrans=false; /* Flagging activation of transition routines */
-    int pointExitCode = 0; /* Output of pseudo time integration (0: converged; -1: unsuccessful, failed; >0: unsuccessful nb iter) */
+    bool lockTrans=false;   /* Flagging activation of transition routines */
+    int pointExitCode = 0;  /* Output of pseudo time integration (0: converged; -1: unsuccessful, failed; >0: unsuccessful nb iter) */
     solverExitCode = 0;
 
+    double maxMach = 0.;
+
     for (size_t iSec=0; iSec<Sections.size(); ++iSec)
     {
+        tms["0-CheckIn"].start();
 
-        if (verbose>0){std::cout << "Sec " << iSec << " ";}
-        if (verbose>1){std::cout << "Mmax " << std::max(Sections[iSec][0]->invBnd->GetMaxMach(), Sections[iSec][1]->invBnd->GetMaxMach());}
-
-        /* Check for stagnation point movement */
+        if (verbose>1){std::cout << "Sec " << iSec << " ";}
 
+        // Check for stagnation point movement.
         if (couplIter != 0 && (Sections[iSec][0]->mesh->GetDiffnElm()))
         {
             stagPtMvmt[iSec] = true;
-            if (verbose>1){std::cout << "Stagnation point moved by " << Sections[iSec][0]->mesh->GetDiffnElm() << " elements." << std::endl;}
+            if (verbose>2) std::cout << "Stagnation point moved by " << Sections[iSec][0]->mesh->GetDiffnElm() << " elements." << std::endl;
         }
         else
-        {
             stagPtMvmt[iSec] = false;
-        }
-
-        /* Set boundary conditions */
 
+        // Set boundary conditions.
         Sections[iSec][0]->xtr = 1.; /* Upper side initial xtr */
         Sections[iSec][1]->xtr = 1.; /* Lower side initial xtr */
         Sections[iSec][2]->xtr = 0.; /* Wake initial xtr (fully turbulent) */
@@ -132,14 +127,17 @@ int ViscSolver::Run(unsigned int const couplIter)
         Sections[iSec][0]->SetStagBC(Re);
         Sections[iSec][1]->SetStagBC(Re);
 
-        /* Loop over the regions (upper, lower and wake) */
+        tms["0-CheckIn"].stop();
+
+        // Loop over the regions (upper, lower and wake).
         for (size_t iRegion = 0; iRegion < Sections[iSec].size(); ++iRegion)
         {
+            if (Sections[iSec][iRegion]->invBnd->GetMaxMach() > maxMach)
+                maxMach = Sections[iSec][iRegion]->invBnd->GetMaxMach();
 
-            if (verbose>0){std::cout << Sections[iSec][iRegion]->name << " ";}
-
-            /* Check if safe mode is required */
+            if (verbose>1) std::cout << Sections[iSec][iRegion]->name << " ";
 
+            // Check if safe mode is required.
             if (couplIter == 0 || Sections[iSec][iRegion]->mesh->GetDiffnElm() != 0 || stagPtMvmt[iSec] == true || solverExitCode != 0)
             {
                 tSolver->SetCFL0(1);
@@ -150,9 +148,7 @@ int ViscSolver::Run(unsigned int const couplIter)
                 {
                     std::vector<double> xxExtCurr(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.);
                     for (size_t iPoint=0; iPoint<xxExtCurr.size(); ++iPoint)
-                    {
                         xxExtCurr[iPoint] = Sections[iSec][iRegion]->mesh->GetLoc(iPoint);
-                    }
                     Sections[iSec][iRegion]->mesh->SetExt(xxExtCurr);
 
                     std::vector<double> DeltaStarZeros(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.0);
@@ -160,14 +156,12 @@ int ViscSolver::Run(unsigned int const couplIter)
 
                     std::vector<double> ueOld(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.0);
                     for (size_t iPoint=0; iPoint<ueOld.size(); ++iPoint)
-                    {
                         ueOld[iPoint] = Sections[iSec][iRegion]->invBnd->GetVt(iPoint);
-                    }
                     Sections[iSec][iRegion]->invBnd->SetVtOld(ueOld);
                 }
 
-                /* Keyword annoncing iteration safety level */
-                if (verbose>1){std::cout << "restart. ";}
+                // Keyword annoncing iteration safety level.
+                if (verbose>1) std::cout << "restart. ";
 
             }
             else
@@ -176,30 +170,25 @@ int ViscSolver::Run(unsigned int const couplIter)
                 tSolver->SetitFrozenJac(5);
                 tSolver->SetinitSln(false);
 
-                /* Keyword annoncing iteration safety level */
-                if (verbose>1){std::cout << "update. ";}
+                // Keyword annoncing iteration safety level.
+                if (verbose>1) std::cout << "update. ";
             }
 
-            /* Save number of points in each regions */
-
+            // Save number of points in each regions.
             Sections[iSec][iRegion]->mesh->prevnMarkers = Sections[iSec][iRegion]->mesh->GetnMarkers();
 
-            /* Reset regime */
-
+            // Reset regime.
             lockTrans = false;
             if (iRegion < 2) // Airfoil
-            {
                 for (size_t k = 0; k < Sections[iSec][iRegion]->mesh->GetnMarkers(); k++)
-                {
                     Sections[iSec][iRegion]->Regime[k] = 0;
-                }
-            }
+                    
             else // Wake
             {
                 for (size_t k = 0; k < Sections[iSec][iRegion]->mesh->GetnMarkers(); k++)
                     Sections[iSec][iRegion]->Regime[k] = 1;
 
-                /* Impose wake boundary condition */
+                // Impose wake boundary condition.
                 std::vector<double> UpperCond(Sections[iSec][iRegion]->GetnVar(), 0.);
                 std::vector<double> LowerCond(Sections[iSec][iRegion]->GetnVar(), 0.);
                 for (size_t k=0; k<UpperCond.size(); ++k)
@@ -209,7 +198,7 @@ int ViscSolver::Run(unsigned int const couplIter)
                 }
                 Sections[iSec][iRegion]->SetWakeBC(Re, UpperCond, LowerCond);
                 lockTrans = true;
-                if (verbose>2){std::cout << "Wake BC imposed." <<std::endl;}
+                if (verbose>2) std::cout << "Wake BC imposed." <<std::endl;
             }
 
             // Loop over points
@@ -219,51 +208,55 @@ int ViscSolver::Run(unsigned int const couplIter)
                 // Initialize solution at point
                 tSolver->InitialCondition(iPoint, Sections[iSec][iRegion]);
                 // Solve equations 
+                tms["1-TimeSolver"].start();
                 pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
+                tms["1-TimeSolver"].stop();
 
                 // Unsucessfull convergence output
-                if (pointExitCode != 0)
+                if (pointExitCode < 0)
                 {
                     if (verbose>1){std::cout << "Point " << iPoint << ": Convergence terminated with exitcode " << pointExitCode << std::endl;}
                     solverExitCode = -1; // Convergence failed
                 }
+                else if (pointExitCode > 0 and solverExitCode >=0)
+                    solverExitCode = pointExitCode; // Not fully converged
                 else
-                {
-                    if (verbose>2){std::cout << "Point " << iPoint << ": Successful convergence.\n";}
-                }
-
-                /* Transition */
+                    if (verbose>2) std::cout << "Point " << iPoint << ": Successful convergence.\n";
 
+                // Transition
+                tms["2-Transition"].start();
                 if (!lockTrans)
                 {
-                    /* Check if perturbation waves are growing or not. */
-
+                    // Check if perturbation waves are growing or not.
                     CheckWaves(iPoint, Sections[iSec][iRegion]);
 
-                    /* Amplification factor is compared to critical amplification factor. */
-
+                    // Amplification factor is compared to critical amplification factor.
                     if (Sections[iSec][iRegion]->U[iPoint * Sections[iSec][iRegion]->GetnVar() + 2] >= Sections[iSec][iRegion]->GetnCrit())
                     {
                         AverageTransition(iPoint, Sections[iSec][iRegion], 0);
-                        if (verbose>0)
+                        if (verbose>1)
                         {   
                             std::cout << std::fixed;
                             std::cout << std::setprecision(2);
                             std::cout << Sections[iSec][iRegion]->xtr
                                       << " (" << Sections[iSec][iRegion]->transMarker << ") ";
                         }
-                        if (verbose>2){std::cout << "" << std::endl;}
+                        if (verbose>2)
+                            std::cout << "" << std::endl;
                         lockTrans = true;
-                    } // End if
-                }     // End if (!lockTrans).
+                    }
+                }
+                tms["2-Transition"].stop();
 
             } // End for iPoints
 
+            tms["3-Blowing"].start();
             Sections[iSec][iRegion]->ComputeBlwVel();
+            tms["3-Blowing"].stop();
             if (verbose>2){"Blowing velocities of " + Sections[iSec][iRegion]->name + " side computed.";}
 
-        }     // End for iRegion
-        if (verbose>0){std::cout << "\n";}
+        } // End for iRegion
+        if (verbose>1) std::cout << "\n";
     } // End for iSections
 
     for (size_t i=0; i<Sections.size(); ++i)
@@ -272,9 +265,39 @@ int ViscSolver::Run(unsigned int const couplIter)
     }
     ComputeTotalDrag();
 
-    if (solverExitCode != 0 && verbose>1){std::cout << "some points did not converge. ";}
-    if (verbose>0){std::cout << "\n" << std::endl;}
+    // Compute mean transition locations
+    std::vector<double> meanTransitions(2, 0.);
+    for (size_t iSec=0; iSec<Sections.size(); ++iSec)
+        for (size_t iReg=0; iReg<2; ++iReg)
+            meanTransitions[iReg] += Sections[iSec][iReg]->xtr;
+    for (size_t iReg=0; iReg<meanTransitions.size(); ++iReg)
+        meanTransitions[iReg]/=Sections.size();
 
+    if (verbose > 0)
+    {
+        if (verbose > 1) std::cout << "\n" << std::endl;
+        std::cout << std::fixed << std::setprecision(2)
+                  << "Viscous solver for Re " << Re/1e6 << "e6, Mach max "
+                  << maxMach
+                  << std::endl;
+
+        std::cout << std::setw(8) << "xTr_Up"
+                  << std::setw(12) << "xTr_Lw"
+                  << std::setw(12) << "Cd" << std::endl;
+        std::cout << "-----------------------------------\n";
+        std::cout << std::fixed << std::setprecision(5);
+        std::cout << std::setw(8) << meanTransitions[0]
+                  << std::setw(12) << meanTransitions[1]
+                  << std::setw(12) << Cdt << "\n";
+
+        if (verbose > 2)
+        {
+            std::cout << "\n" << std::endl;
+            std::cout << "Viscous solver CPU" << std::endl
+                        << tms;
+        }
+        std::cout << "\n" << std::endl;
+    }
     return solverExitCode; // exit code
 }
 
@@ -283,16 +306,13 @@ int ViscSolver::Run(unsigned int const couplIter)
  */
 void ViscSolver::CheckWaves(size_t iPoint, BLRegion *bl)
 {
-    /* Determine if the amplification waves start growing or not. */
 
     double logRet_crit = 2.492 * std::pow(1 / (bl->Hk[iPoint] - 1), 0.43) + 0.7 * (std::tanh(14 * (1 / (bl->Hk[iPoint] - 1)) - 9.24) + 1.0);
 
     if (bl->Ret[iPoint] > 0)
     { // Avoid log of negative number.
         if (std::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.
-        }
     }
 }
 
@@ -301,11 +321,11 @@ void ViscSolver::CheckWaves(size_t iPoint, BLRegion *bl)
  */
 void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
 {
-    /* Averages solution on transition marker */
+    // Averages solution on transition marker.
 
     unsigned int nVar = bl->GetnVar();
 
-    /* Save laminar solution. */
+    // Save laminar solution.
 
     std::vector<double> lamSol(nVar, 0.0);
     for (size_t k = 0; k < lamSol.size(); ++k)
@@ -313,31 +333,29 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
         lamSol[k] = bl->U[iPoint * nVar + k];
     }
 
-    /* Set up turbulent portion boundary condition. */
+    // Set up turbulent portion boundary condition.
 
-    bl->transMarker = iPoint; /* Save transition marker */
+    bl->transMarker = iPoint; // Save transition marker.
 
-    /* Loop over remaining points in the region */
+    // Loop over remaining points in the region.
     for (size_t k = iPoint; k < bl->mesh->GetnMarkers(); ++k)
-    {
         bl->Regime[k] = 1; // Set Turbulent regime.
-    }
 
-    /* Compute transition location */
+    // Compute transition location.
     bl->xtr = (bl->GetnCrit() - bl->U[(iPoint - 1) * nVar + 2]) * ((bl->mesh->Getx(iPoint) - bl->mesh->Getx(iPoint - 1)) / (bl->U[iPoint * nVar + 2] - bl->U[(iPoint - 1) * nVar + 2])) + bl->mesh->Getx(iPoint - 1);
 
-    /*  Percentage of the subsection corresponding to a turbulent flow */
+    // Percentage of the subsection corresponding to a turbulent flow.
     double avTurb = (bl->mesh->Getx(iPoint) - bl->xtr) / (bl->mesh->Getx(iPoint) - bl->mesh->Getx(iPoint - 1));
     double avLam = 1. - avTurb;
 
-    /* Impose boundary condition */
+    // Impose boundary condition.
     double Cteq_trans;
     bl->closSolver->TurbulentClosures(Cteq_trans, bl->U[(iPoint - 1) * nVar + 0], bl->U[(iPoint - 1) * nVar + 1], 0.03, bl->Ret[iPoint - 1], bl->invBnd->GetMe(iPoint - 1), bl->name);
     bl->Cteq[iPoint - 1] = Cteq_trans;
     bl->U[(iPoint - 1) * nVar + 4] = 0.7 * bl->Cteq[iPoint - 1];
 
-    /* Avoid starting with ill conditioned IC for turbulent computation. (Regime was switched above) */
-    /* These initial guess values do not influence the converged solution */
+    // Avoid starting with ill conditioned IC for turbulent computation. (Regime was switched above).
+    // These initial guess values do not influence the converged solution.
     bl->U[iPoint * nVar + 4] = bl->U[(iPoint - 1) * nVar + 4]; // IC of transition point = turbulent BC (imposed @ previous point).
     bl->U[iPoint * nVar + 1] = 1.515;                          // Because we expect a significant drop of the shape factor.
 
@@ -345,19 +363,18 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
     int exitCode = tSolver->Integration(iPoint, bl);
 
     if (exitCode != 0 && verbose>1)
-    {
         std::cout << "Warning: Transition point turbulent computation terminated with exit code " << exitCode << std::endl;
-    }
-    else if(verbose>2){std::cout << "Sucessfull transition point turbulent computation.";}
+    else if(verbose>2)
+        std::cout << "Sucessfull transition point turbulent computation.";
 
-    /* Average both solutions (Now solution @ iPoint is the turbulent solution) */
-    bl->U[iPoint * nVar + 0] = avLam * lamSol[0] + avTurb * bl->U[(iPoint)*nVar + 0]; /* Theta */
-    bl->U[iPoint * nVar + 1] = avLam * lamSol[1] + avTurb * bl->U[(iPoint)*nVar + 1]; /* H */
-    bl->U[iPoint * nVar + 2] = 9.;                                                     /* N */
-    bl->U[iPoint * nVar + 3] = avLam * lamSol[3] + avTurb * bl->U[(iPoint)*nVar + 3]; /* ue */
-    bl->U[iPoint * nVar + 4] = avTurb * bl->U[(iPoint)*nVar + 4];                     /* Ct */
+    // Average both solutions (Now solution @ iPoint is the turbulent solution).
+    bl->U[iPoint * nVar + 0] = avLam * lamSol[0] + avTurb * bl->U[(iPoint)*nVar + 0];   /* Theta */
+    bl->U[iPoint * nVar + 1] = avLam * lamSol[1] + avTurb * bl->U[(iPoint)*nVar + 1];   /* H */
+    bl->U[iPoint * nVar + 2] = 9.;                                                      /* N */
+    bl->U[iPoint * nVar + 3] = avLam * lamSol[3] + avTurb * bl->U[(iPoint)*nVar + 3];   /* ue */
+    bl->U[iPoint * nVar + 4] = avTurb * bl->U[(iPoint)*nVar + 4];                       /* Ct */
 
-    /* Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations). */
+    // Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations).
     std::vector<double> lamParam(8, 0.);
     bl->closSolver->LaminarClosures(lamParam, bl->U[(iPoint - 1) * nVar + 0], bl->U[(iPoint - 1) * nVar + 1], bl->Ret[iPoint - 1], bl->invBnd->GetMe(iPoint - 1), bl->name);
     bl->Hstar[iPoint - 1] = lamParam[0];
@@ -394,40 +411,29 @@ void ViscSolver::ComputeSectionalDrag(std::vector<BLRegion *> const &bl, int i)
     unsigned int nVar = bl[0]->GetnVar();
     unsigned int lastWkPt = (bl[2]->mesh->GetnMarkers() - 1) * nVar;
 
-    /* Normalize friction and dissipation coefficients. */
-
+    // Normalize friction and dissipation coefficients.
     std::vector<double> frictioncoeff_upper(bl[0]->mesh->GetnMarkers(), 0.0);
     for (size_t k = 0; k < frictioncoeff_upper.size(); ++k)
-    {
         frictioncoeff_upper[k] = bl[0]->invBnd->GetVt(k) * bl[0]->invBnd->GetVt(k) * bl[0]->Cf[k];
-    }
+
     std::vector<double> frictioncoeff_lower(bl[1]->mesh->GetnMarkers(), 0.0);
     for (size_t k = 0; k < frictioncoeff_lower.size(); ++k)
-    {
         frictioncoeff_lower[k] = bl[1]->invBnd->GetVt(k) * bl[1]->invBnd->GetVt(k) * bl[1]->Cf[k];
-    }
-
-    /* Compute friction drag coefficient. */
 
+    // Compute friction drag coefficient.
     double Cdf_upper = 0.0;
     for (size_t k = 1; k < bl[0]->mesh->GetnMarkers(); ++k)
-    {
         Cdf_upper += (frictioncoeff_upper[k] + frictioncoeff_upper[k - 1]) * (bl[0]->mesh->Getx(k) - bl[0]->mesh->Getx(k - 1)) / 2;
-    }
     double Cdf_lower = 0.0;
     for (size_t k = 1; k < bl[1]->mesh->GetnMarkers(); ++k)
-    {
         Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * (bl[1]->mesh->Getx(k) - bl[1]->mesh->Getx(k - 1)) / 2;
-    }
 
     Cdf_sec[i] = Cdf_upper + Cdf_lower; // No friction in the wake
 
-    /* Compute total drag coefficient. */
-
+    // Compute total drag coefficient.
     Cdt_sec[i] = 2 * bl[2]->U[lastWkPt + 0] * pow(bl[2]->U[lastWkPt + 3], (bl[2]->U[lastWkPt + 1] + 5) / 2);
 
-    /* Compute pressure drag coefficient. */
-
+    // Compute pressure drag coefficient.
     Cdp_sec[i] = Cdt_sec[i] - Cdf_sec[i];
 }
 
@@ -445,14 +451,12 @@ void ViscSolver::ComputeTotalDrag()
 
     else
     {
-        std::cout << "0 " << Cdt_sec[0] << std::endl;
         Cdt += (Sections[0][0]->mesh->Getz(0) - 0) * Cdt_sec[0];
         Cdp += (Sections[0][0]->mesh->Getz(0) - 0) * Cdp_sec[0];
         Cdf += (Sections[0][0]->mesh->Getz(0) - 0) * Cdf_sec[0];
         double dz = 0.;
         for(size_t k=1; k<Sections.size(); ++k)
         {
-            std::cout << k << " " << Cdt_sec[k] << std::endl;
             dz = (Sections[k][0]->mesh->Getz(0) - Sections[k-1][0]->mesh->Getz(0));
             Cdt += dz * Cdt_sec[k];
             Cdp += dz * Cdp_sec[k];
@@ -465,9 +469,9 @@ void ViscSolver::ComputeTotalDrag()
 
 void ViscSolver::SetGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z)
 {
-    if (verbose>1){std::cout << "Mesh of " << _x.size() << " elements on " << Sections[iSec][iRegion]->name << " side. ";}
+    if (verbose>2){std::cout << "Mesh of " << _x.size() << " elements on " << Sections[iSec][iRegion]->name << " side. ";}
     Sections[iSec][iRegion]->SetMesh(_x, _y, _z);
-    if (verbose>1){std::cout << "Mesh mapped to local frame of reference." << std::endl;}
+    if (verbose>3){std::cout << "Mesh mapped to local frame of reference." << std::endl;}
 }
 
 void ViscSolver::SetVelocities(size_t iSec, size_t iRegion, std::vector<double> vx, std::vector<double> vy, std::vector<double> vz)
@@ -495,27 +499,21 @@ void ViscSolver::SetVelocities(size_t iSec, size_t iRegion, std::vector<double>
 void ViscSolver::SetMe(size_t iSec, size_t iRegion, std::vector<double> _Me)
 {
     if (_Me.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
-    {
         throw std::runtime_error("dart::ViscSolver Mach number vector is not coherent with the mesh.");
-    }
     Sections[iSec][iRegion]->invBnd->SetMe(_Me);
 }
 
 void ViscSolver::SetRhoe(size_t iSec, size_t iRegion, std::vector<double> _Rhoe)
 {
     if (_Rhoe.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
-    {
         throw std::runtime_error("dart::ViscSolver Density vector is not coherent with the mesh.");
-    }
     Sections[iSec][iRegion]->invBnd->SetRhoe(_Rhoe);
 }
 
 void ViscSolver::SetxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt)
 {   
     if (_xxExt.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
-    {
         throw std::runtime_error("dart::ViscSolver External mesh vector is not coherent with the mesh.");
-    }
     Sections[iSec][iRegion]->mesh->SetExt(_xxExt);
 }
 
@@ -533,9 +531,7 @@ std::vector<double> ViscSolver::ExtractxCoord(size_t iSec, size_t iRegion) const
 {
     std::vector<double> xCoord(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.);
     for (size_t iPoint=0; iPoint<xCoord.size(); ++iPoint)
-    {
         xCoord[iPoint] = Sections[iSec][iRegion]->mesh->Getx(iPoint);
-    }
     return xCoord;
 }
 
@@ -547,9 +543,7 @@ std::vector<double> ViscSolver::Extractxx(size_t iSec, size_t iRegion) const
 {   
     std::vector<double> xx(Sections[iSec][iRegion]->mesh->GetnMarkers(),0.);
     for (size_t iPoint=0; iPoint<xx.size(); ++iPoint)
-    {
         xx[iPoint] = Sections[iSec][iRegion]->mesh->GetLoc(iPoint);
-    }
     return xx;
 }
 std::vector<double> ViscSolver::ExtractDeltaStar(size_t iSec, size_t iRegion) const
@@ -561,18 +555,14 @@ std::vector<double> ViscSolver::ExtractUe(size_t iSec, size_t iRegion) const
 {
     std::vector<double> ueCurr(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.0);
     for (size_t i=0; i<ueCurr.size(); ++i)
-    {
         ueCurr[i] = Sections[iSec][iRegion]->U[i * Sections[iSec][iRegion]->GetnVar() + 3];
-    }
     return ueCurr;
 }
 
 void ViscSolver::SetUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue)
 {
     if (_ue.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
-    {
         throw std::runtime_error("dart::ViscSolver External velocity vector is not coherent with the mesh.");
-    }
     Sections[iSec][iRegion]->invBnd->SetVtOld(_ue);
 }
 
@@ -584,9 +574,7 @@ std::vector<std::vector<double>> ViscSolver::GetSolution(size_t iSec)
 
     std::vector<std::vector<double>> Sln(16);
     for (size_t iVector=0; iVector<Sln.size(); ++iVector)
-    {
         Sln[iVector].resize(nMarkersUp + nMarkersLw + Sections[iSec][2]->mesh->GetnMarkers(), 0.);
-    }
     /* Upper side */
     size_t revPoint = nMarkersUp - 1;
     for (size_t iPoint=0; iPoint<nMarkersUp; ++iPoint)
@@ -655,10 +643,9 @@ std::vector<std::vector<double>> ViscSolver::GetSolution(size_t iSec)
         Sln[14][nMarkersUp + nMarkersLw + iPoint] = Sections[iSec][2]->Delta[iPoint];
 
         Sln[15][nMarkersUp + nMarkersLw + iPoint] = Sections[iSec][2]->mesh->Getx(iPoint);
-
     }
 
-    if (verbose>2){std::cout << "Solution structure sent." << std::endl;}
+    if (verbose>2) std::cout << "Solution structure sent." << std::endl;
 
     return Sln;
 
diff --git a/vii/src/wViscSolver.h b/vii/src/wViscSolver.h
index a981d88c0e33951c6f5581e8d2080ae4c716fda9..2a717f7822e5bc31e25f8f74a26aa639e7af96f8 100644
--- a/vii/src/wViscSolver.h
+++ b/vii/src/wViscSolver.h
@@ -3,6 +3,7 @@
 
 #include "vii.h"
 #include "wBLRegion.h"
+#include "wTimers.h"
 namespace vii
 {
 
@@ -33,6 +34,9 @@ private:
 
     double Span; /* Wing Span (0 if 2D case) */
 
+protected:
+    fwk::Timers tms;       ///< internal timers
+
 public:
     ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSections, double _Span=0., unsigned int verbose=1);
     ~ViscSolver();
diff --git a/vii/tests/bli25.py b/vii/tests/bli25.py
index c2129b908e6b58b7816ab2075a5088590fa907b8..b53fc2355df7e16355a3af3d253713ff3226f06b 100644
--- a/vii/tests/bli25.py
+++ b/vii/tests/bli25.py
@@ -64,7 +64,7 @@ def main():
     c_ref = 1.0 # reference length
     S_ref = spn # reference area
     fms = 1.0 # farfield mesh size
-    nms = 0.01 # nearfield mesh size
+    nms = 0.007 # nearfield mesh size
 
     # mesh the geometry
     print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
@@ -82,7 +82,13 @@ def main():
     pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0.25, 0., 0., 'wing', vsc = True)
     tms['pre'].stop()
 
-    config = {'nDim': dim, 'Sections':[-0.1, -0.05, 0.0, 0.05, 0.1], 'span':spn, 'Sym':[-0.1, 0.1], 'rbftype': 'linear', 'smoothing': 1e-6, 'degree': 0, 'neighbors': 10}
+    config = {'nDim': dim,
+              'Sections':[-0.1, -0.05, 0.0, 0.05, 0.1],
+              'span':spn, 'Sym':[-0.1, 0.1],
+              'rbftype': 'linear',
+              'smoothing': 1e-6,
+              'degree': 0,
+              'neighbors': 10}
     
     # solve problem
     iSolver = floD.newton(pbl)
@@ -99,7 +105,7 @@ def main():
 
     # VII
     iSolverAPI = vInterpol.Interpolator(iSolver, blw, imsh, vmsh, config)
-    vSolver = viscU.initBL(Re, M_inf, CFL0, len(config['Sections']), spn, 2)
+    vSolver = viscU.initBL(Re, M_inf, CFL0, len(config['Sections']), spn, 1)
     coupler = viscC.Coupler(iSolverAPI, vSolver, _maxIters=200)
 
     coupler.Run()
@@ -107,9 +113,6 @@ def main():
     vSolution = viscU.GetSolution(0, vSolver)
     viscU.WriteFile(vSolution, vSolver.GetRe())
 
-    plt.plot(vSolution['x/c'], vSolution['theta'])
-    plt.show()
-
     # Write corrected Cp
     iSolver.save(mshWriter)
     dartU.writeSlices(msh.name, config['Sections'], 2)
diff --git a/vii/tests/bli25_RandomTest.py b/vii/tests/bli25_RandomTest.py
new file mode 100644
index 0000000000000000000000000000000000000000..642d3760909eeedd6b19fc9ded6396a7bb511af5
--- /dev/null
+++ b/vii/tests/bli25_RandomTest.py
@@ -0,0 +1,161 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## Compute lifting (linear or nonlinear) potential flow around a NACA 0012 rectangular wing
+# Adrien Crovato
+#
+# Test the nonlinear shock-capturing capability and the Kutta condition
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import math
+import dart.default as floD
+import dart.utils as dartU # or from dartflo.dart import utils as dartU
+import tbox.utils as tboxU # or from dartflo.tbox import utils as tboxU
+import tboxVtk
+import numpy as np
+
+import vii.pyVII.vUtils as viscU
+import vii.pyVII.vCoupler2 as viscC
+import vii.pyVII.vInterpolator as vInterpol
+
+from matplotlib import pyplot as plt
+
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+def main():
+    # timer
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    # define flow variables
+    alpha = 2.*np.pi/180
+    M_inf = 0.715
+    dim = 3
+
+    CFL0 = 1
+    Re = 1e7
+    dim = 3
+
+    # define dimension and mesh size
+    spn = 0.2 # wing span
+    lgt = 6.0 # channel length
+    hgt = 6.0 # channel height
+    wdt = 3.0 # channel width
+    c_ref = 1.0 # reference length
+    S_ref = spn # reference area
+    fms = 1.0 # farfield mesh size
+    nms = 0.01 # nearfield mesh size
+    
+    U = []
+    for i in range(5):
+        # mesh the geometry
+        print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+        tms["msh"].start()
+        pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
+        imsh = viscU.mesh('/home/paul/lab/dartflo/dart/models/rae_25.msh', pars)
+        vmsh = viscU.mesh('/home/paul/lab/dartflo/vii/models/rae_25_visc.geo', pars)
+        msh, gmshWriter = floD.mesh(dim, 'models/rae_25.msh', pars, ['field', 'wing', 'symmetry', 'downstream'])
+        # morpher = floD.morpher(msh, dim, 'wing')
+        tms["msh"].stop()
+
+        # set the problem and add fluid, initial/boundary conditions
+        tms['pre'].start()
+        # Replce tp = 'teTip' by tp = 'wakeTip' for oneraM6.geo
+        pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0.25, 0., 0., 'wing', vsc = True)
+        tms['pre'].stop()
+
+        config = {'nDim': dim, 'Sections':[-0.1, -0.05, 0.0, 0.05, 0.1], 'span':spn, 'Sym':[-0.1, 0.1], 'rbftype': 'linear', 'smoothing': 1e-6, 'degree': 0, 'neighbors': 10}
+        
+        # solve problem
+        iSolver = floD.newton(pbl)
+
+        # VII
+        iSolverAPI = vInterpol.Interpolator(iSolver, blw, imsh, vmsh, config)
+        vSolver = viscU.initBL(Re, M_inf, CFL0, len(config['Sections']), spn, 2)
+        coupler = viscC.Coupler(iSolverAPI, vSolver, _maxIters=1)
+
+        coupler.Run()
+
+        U.append(iSolverAPI.URani)
+    
+    for i in range(len(U)):
+        plt.plot((U[i]-U[0])/U[0])
+    plt.show()
+
+    plt.show()
+    vSolution = viscU.GetSolution(0, vSolver)
+    viscU.WriteFile(vSolution, vSolver.GetRe())
+
+    plt.plot(vSolution['x/c'], vSolution['theta'])
+    plt.show()
+
+    # Write corrected Cp
+    iSolver.save(mshWriter)
+    dartU.writeSlices(msh.name, config['Sections'], 2)
+
+    cpCorr = [np.zeros((0,0)) for _ in range(len(config['Sections']))]
+    for i in range(len(config['Sections'])):
+        cpCorr[i] = tboxU.read('rae_25_slice_'+str(i)+'.dat')
+
+    # Plot cp
+
+    for i in range(len(config['Sections'])):
+        plt.plot(cpCorr[i][:,3], cpCorr[i][:,4], '-')
+        plt.plot(cpInv[i][:,3], cpInv[i][:,4], '--')
+        plt.xlabel('$x/c$')
+        plt.ylabel('$Cp$')
+        plt.title('$y/b$ = {:.3f}'.format(config['Sections'][i]))
+        plt.gca().invert_yaxis()
+        plt.legend(['Corrected', 'Inviscid'])
+        plt.show()
+
+    # display timers
+    tms['total'].stop()
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    print(coupler.tms)
+
+    # visualize solution
+    iSolverAPI.GetCp()
+    iSolverAPI.iSolver.save(gmshWriter)
+    floD.initViewer(pbl)
+
+    # check results
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    if alpha == 3*math.pi/180 and M_inf == 0.3 and spn == 1.0:
+        tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True))
+        tests.add(CTest('CL', solver.Cl, 0.135, 5e-2))
+        tests.add(CTest('CD', solver.Cd, 0.0062, 1e-2)) # Tranair (NF=0.0062, FF=0.0030), Panair 0.0035
+        tests.add(CTest('CS', solver.Cs, 0.0121, 5e-2))
+        tests.add(CTest('CM', solver.Cm, -0.0278, 1e-2))
+    else:
+        raise Exception('Test not defined for this flow')
+    tests.run()
+
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()