diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index e314c7285521f7b3c0f05d72bff983d4c11e96d3..916764c256db144782179057c8a895b2026c57ca 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -64,27 +64,29 @@ void TimeSolver::InitialCondition(size_t iPoint, BLRegion *bl)
 int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
 {
 
-  // Save initial condition
-
   unsigned int nVar = bl->GetnVar();
+  
+  /* Save initial condition */
 
   std::vector<double> Sln0;
   for (size_t i=0; i<nVar; ++i){
     Sln0.push_back(bl->U[iPoint*nVar+i]);
   } // End for
 
-  // Initialize time integration variables
+  /* Initialize time integration variables */
 
   double dx = bl->GetLocMarkers(iPoint) - bl->GetLocMarkers(iPoint-1);
 
+  /* Initial time step */
+
   double CFL = CFL0;
   unsigned int itFrozenJac = itFrozenJac0;
-
   double timeStep = SetTimeStep(CFL, dx, bl->GetVtExt(iPoint));
-
   Matrix<double, 5, 5> IMatrix;
   IMatrix = Matrix<double, 5, 5>::Identity() / timeStep;
 
+  /* Initial system residuals */
+
   Vector<double, 5> SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
   double normSysRes0 = SysRes.norm();
   double normSysRes = normSysRes0;
@@ -97,13 +99,6 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
 
   while (innerIters < maxIt){
 
-    if (bl->name == "lower" && iPoint == 425)
-    {
-      std::cout << "iter " << innerIters << std::endl;
-      bl->PrintSolution(iPoint-1);
-      std::cout << (normSysRes/normSysRes0) << " CFL" << CFL << std::endl;
-    }
-
     if(innerIters % itFrozenJac == 0)
     {
       if (nErrors >= 5)
@@ -198,7 +193,6 @@ double TimeSolver::ControlIntegration(size_t iPoint, BLRegion *bl, std::vector<d
 
   if (bl->U[iPoint*nVar + 3] <= 0)
   {
-    std::cout << "Negative Ue detected at point " << iPoint << std::endl;
     ResetSolution(iPoint, bl, Sln0, true);
     nErrors+=1;
     return 0.2;
@@ -264,12 +258,6 @@ void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double>
     bl->Us[iPoint] = turbParam[7];
   } // End else if
 
-  if (bl->name=="lower" && iPoint == 425)
-  {
-    std::cout << "end Reset" << std::endl;
-    bl->PrintSolution(iPoint);
-  }
-
 } // End ResetSolution
 
 
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index a22bc6e2c101b643dec131a578ede7e460821267..1c68fee80b27ae0c509398a59d98aa6be947b225 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -71,12 +71,13 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
     dissipRatio = 1;
   } // End else
   
-  bl->Ret[iPoint] = std::max(u[3] * u[0] * Re, 1.0);
-  bl->DeltaStar[iPoint] = u[0] * u[1];
+  bl->Ret[iPoint] = std::max(u[3] * u[0] * Re, 1.0); /* Reynolds number based on theta */
+  bl->DeltaStar[iPoint] = u[0] * u[1];               /* Displacement thickness */
 
   /* Boundary layer closure relations */
 
   if (bl->Regime[iPoint] == 0){
+    /* Laminar closure relations */
     std::vector<double> lamParam = bl->closSolver->LaminarClosures(u[0], u[1], bl->Ret[iPoint], bl->GetMe(iPoint), bl->name);
     bl->Hstar[iPoint] = lamParam[0];
     bl->Hstar2[iPoint] = lamParam[1];
@@ -87,7 +88,9 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
     bl->Cteq[iPoint] = lamParam[6];
     bl->Us[iPoint] = lamParam[7];
     } // End if
+    
   else if (bl->Regime[iPoint] == 1){
+    /* Turbulent closure relations */
     std::vector<double> turbParam = bl->closSolver->TurbulentClosures(u[0], u[1], u[4], bl->Ret[iPoint], bl->GetMe(iPoint), bl->name);
     bl->Hstar[iPoint] = turbParam[0];
     bl->Hstar2[iPoint] = turbParam[1];
@@ -98,7 +101,6 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
     bl->Cteq[iPoint] = turbParam[6];
     bl->Us[iPoint] = turbParam[7];
   } // End else if
-
   else{
     std::cout << "Wrong regime\n" << std::endl;
   } // End else
@@ -143,32 +145,8 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
   double Cteqa = bl->Cteq[iPoint];
   double Usa = bl->Us[iPoint];
 
-  /*if (bl->name=="lower" && iPoint == 2){
-    std::cout << "Reta" << bl->Ret[iPoint] << std::endl;
-    std::cout << "Hstara" << Hstara << std::endl;
-    std::cout << "Hstar2a" << Hstar2a << std::endl;
-    std::cout << "Hka" << Hka << std::endl;
-    std::cout << "Cfa" << Cfa << std::endl;
-    std::cout << "Cda" << Cda << std::endl;
-    std::cout << "deltaa" << deltaa << std::endl;
-    std::cout << "Usa" << Usa << std::endl;
-  }*/
-
   /* Space part */
 
-  /*if(bl->name=="lower" && iPoint==2){
-    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;
@@ -221,25 +199,7 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
     timeMatrix(4,4) = 1.0;
   }
 
-  /*if(bl->name=="lower" && iPoint==1){
-    std::cout << "timeMatrix" << timeMatrix << std::endl;
-    std::cout << "spaceVector" << spaceVector << std::endl;
-
-    std::cout << "dT_dx" << dT_dx << std::endl;
-    std::cout << "dH_dx" << dH_dx << std::endl;
-    std::cout << "due_dx" << due_dx << std::endl;
-    std::cout << "dCt_dx" << dCt_dx << std::endl;
-    std::cout << "dueExt_dx" << dueExt_dx << std::endl;
-    std::cout << "ddStarExt_dx" << ddeltaStarExt_dx << std::endl;
-    std::cout << "dHstar_dx" << dHstar_dx << std::endl;
-  }*/
-  Vector<double, 5> sysRes;
-  sysRes = timeMatrix.partialPivLu().solve(spaceVector);
-  /*if(bl->name=="upper" && iPoint==1){
-    std::cout << "sysRes" << sysRes << std::endl;
-  }*/
-
-  return sysRes;
+  return timeMatrix.partialPivLu().solve(spaceVector);
 }
 
 double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta)
diff --git a/dart/tests/bliLowRe.py b/dart/tests/bliLowRe.py
deleted file mode 100755
index 5e67e4e77c461ae34053f05b6c9f2deeb7aaeff2..0000000000000000000000000000000000000000
--- a/dart/tests/bliLowRe.py
+++ /dev/null
@@ -1,155 +0,0 @@
-#!/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) viscous flow around a NACA 0012
-# Amaury Bilocq
-#
-# Test the viscous-inviscid interaction scheme
-# Reference to the master's thesis: http://hdl.handle.net/2268/252195
-# Reference test cases with Naca0012 (different from master's thesis):
-# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
-#    -> msLe = 0.001, xtrTop = 0.061
-# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
-# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
-#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import dart.utils as floU
-import dart.default as floD
-
-import dart.viscous.Master.DCoupler as floVC
-import dart.viscous.Drivers.DDriver as floVSMaster
-
-import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
-import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
-
-import dart.viscous.Physics.DValidation as validation
-
-import tbox
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-import math
-
-def main():
-    # timer
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    # define flow variables
-    Re = 5e5
-    alpha = 2.*math.pi/180
-    M_inf = 0.2
-
-    #user_xtr=[0.41,0.41]
-    #user_xtr=[0,None]
-    user_xtr=[None, None]
-    user_Ti=None
-
-    mapMesh = 1
-    nFactor = 2
-
-    # Time solver parameters
-    Time_Domain = True # True for unsteady solver, False for steady solver
-    CFL0 = 1
-
-    # ========================================================================================
-
-    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
-    c_ref = 1
-    dim = 2
-    tol = 1e-4 # tolerance for the VII (usually one drag count)
-    case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
-
-    # 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.0001}
-    #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()
-
-    # set the problem and add medium, initial/boundary conditions
-    tms['pre'].start()
-    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
-    tms['pre'].stop()
-
-    # solve viscous problem
-
-    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-    tms['solver'].start()
-    isolver = floD.newton(pbl)
-
-    # Choose between steady and unsteady solver
-    if Time_Domain is True: # Run unsteady solver
-        vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
-        coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
-
-    elif Time_Domain is False: # Run steady solver 
-        vsolver=floVS_Std.Solver(Re)
-        coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
-    coupler.run()
-    tms['solver'].stop()
-
-    # extract Cp
-    Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
-    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
-    # display results
-    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
-    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
-    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
-
-    # display timers
-    tms['total'].stop()
-    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-    print('CPU statistics')
-    print(tms)
-    
-    # check results
-    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
-    tests = CTests()
-    if Re == 5e5 and M_inf == 0.2 and alpha == 2*math.pi/180:
-        tests.add(CTest('Cl', isolver.Cl, 0.2206, 5e-2)) # XFOIL 0.2135
-        tests.add(CTest('Cd', vsolver.Cd, 0.007233, 1e-3, forceabs=True)) # XFOIL 0.00706
-        tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 1e-3, forceabs=True)) # 0.00238
-        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.58, 0.03, forceabs=True)) # XFOIL 0.5642
-        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.90, 0.03, forceabs=True)) # XFOIL 0.9290
-    else:
-        print(ccolors.ANSI_RED + 'Test not defined for this flow', ccolors.ANSI_RESET)
-
-    tests.run()
-
-    # visualize solution and plot results
-    floD.initViewer(pbl)
-    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
-
-    val=validation.Validation(case)
-    val.main(isolver.Cl,vsolver.wData)
-
-    # eof
-    print('')
-
-if __name__ == "__main__":
-    main()
diff --git a/dart/tests/bliPolar.py b/dart/tests/bliPolar.py
deleted file mode 100755
index c4d6fd560b10d6f31071377c703826d7c9925471..0000000000000000000000000000000000000000
--- a/dart/tests/bliPolar.py
+++ /dev/null
@@ -1,167 +0,0 @@
-#!/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) viscous flow around a NACA 0012
-# Amaury Bilocq
-#
-# Test the viscous-inviscid interaction scheme
-# Reference to the master's thesis: http://hdl.handle.net/2268/252195
-# Reference test cases with Naca0012 (different from master's thesis):
-# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
-#    -> msLe = 0.001, xtrTop = 0.061
-# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
-# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
-#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-from numpy.core.function_base import linspace
-import dart.utils as floU
-import dart.default as floD
-
-import dart.viscous.Master.DCoupler as floVC
-import dart.viscous.Drivers.DDriver as floVSMaster
-
-import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
-import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
-
-import dart.viscous.Physics.DValidation as validation
-
-import tbox
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-import numpy as np
-from matplotlib import pyplot as plt
-
-import math
-
-def main():
-    # timer
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    # define flow variables
-    Re = 6e6
-    #alpha = 0.*math.pi/180
-    M_inf = 0.15
-
-    #user_xtr=[0.41,0.41]
-    #user_xtr=[0,None]
-    user_xtr=[0, 0]
-    user_Ti=None
-
-    mapMesh = 0
-    nFactor = 2
-
-    # Time solver parameters
-    Time_Domain = True # True for unsteady solver, False for steady solver
-    CFL0 = 1
-
-    # ========================================================================================
-
-    #U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
-    c_ref = 1
-    dim = 2
-    tol = 1e-4 # tolerance for the VII (usually one drag count)
-    case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
-
-    # 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.0001}
-    #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()
-
-    alphaVect = [
-   15.2600*math.pi/180,
-   16.3000*math.pi/180,
-   16.5*math.pi/180,
-   17.1300*math.pi/180,
-   17.2*math.pi/180,
-   17.5*math.pi/180,
-   17.7*math.pi/180,
-   18.0200*math.pi/180]
-
-    alphaSuc = []
-    clSuc = []
-    cdSuc = []
-    alphaFailed = []
-    clFailed = []
-    cdFailed = []
-
-    for alpha in alphaVect:
-      # set the problem and add medium, initial/boundary conditions
-      tms['pre'].start()
-      pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
-      tms['pre'].stop()
-
-      # solve viscous problem
-
-      print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-      tms['solver'].start()
-      isolver = floD.newton(pbl)
-
-      # Choose between steady and unsteady solver
-      if Time_Domain is True: # Run unsteady solver
-          vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
-          coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
-
-      elif Time_Domain is False: # Run steady solver 
-          vsolver=floVS_Std.Solver(Re)
-          coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
-      try:
-        coupler.run()
-        tms['solver'].stop()
-
-        # extract Cp
-        Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
-        tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
-        # display results
-        print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
-        print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
-        print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
-
-        # display timers
-        tms['total'].stop()
-        print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-        print('CPU statistics')
-        print(tms)
-
-        alphaSuc.append(alpha*180/math.pi)
-        clSuc.append(isolver.Cl)
-        cdSuc.append(vsolver.Cd)
-      except KeyboardInterrupt:
-        quit()
-      except:
-        alphaFailed.append(alpha*180/math.pi)
-        clFailed.append(isolver.Cl)
-        cdFailed.append(vsolver.Cd)
-
-
-    print(alphaSuc, clSuc, cdSuc)
-    print(alphaFailed, clFailed, cdFailed)
-
-if __name__ == "__main__":
-    main()
diff --git a/dart/tests/bliSeparated.py b/dart/tests/bliSeparated.py
deleted file mode 100755
index 49822dec4b639b2d9676cbc7bcccc5075654de51..0000000000000000000000000000000000000000
--- a/dart/tests/bliSeparated.py
+++ /dev/null
@@ -1,157 +0,0 @@
-#!/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) viscous flow around a NACA 0012
-# Amaury Bilocq
-#
-# Test the viscous-inviscid interaction scheme
-# Reference to the master's thesis: http://hdl.handle.net/2268/252195
-# Reference test cases with Naca0012 (different from master's thesis):
-# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
-#    -> msLe = 0.001, xtrTop = 0.061
-# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
-# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
-#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import dart.utils as floU
-import dart.default as floD
-
-import dart.viscous.Master.DCoupler as floVC
-import dart.viscous.Drivers.DDriver as floVSMaster
-
-import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
-import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
-
-import dart.viscous.Physics.DValidation as validation
-
-import tbox
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-import math
-
-def main():
-    # timer
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    # define flow variables
-    Re = 1e7
-    alpha = 15*math.pi/180
-    M_inf = 0.2
-    user_xtr=[None, None]
-    user_Ti=None
-
-    mapMesh = 1
-    nFactor = 2
-
-    # Time solver parameters
-    Time_Domain=True # True for unsteady solver, False for steady solver
-    CFL0 = 1
-
-    # ========================================================================================
-
-    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
-    c_ref = 1
-    dim = 2
-    tol = 1e-4 # tolerance for the VII (usually one drag count)
-    case='case15Xfoil.dat' # File containing results from XFOIL of the considered case
-
-
-    # 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.0001}
-    msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
-    tms['msh'].stop()
-
-    # set the problem and add medium, initial/boundary conditions
-    tms['pre'].start()
-    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
-    tms['pre'].stop()
-
-    # solve viscous problem
-    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-    tms['solver'].start()
-    isolver = floD.newton(pbl)
-
-  # Choose between steady and unsteady solver
-    if Time_Domain is True: # Run unsteady solver
-        vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
-        coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
-
-    elif Time_Domain is False: # Run steady solver 
-        vsolver=floVS_Std.Solver(Re)
-        coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
-    coupler.run()
-    tms['solver'].stop()
-
-    # extract Cp
-    Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
-    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
-    # display results
-    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
-    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
-    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
-
-    # display timers
-    tms['total'].stop()
-    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-    print('CPU statistics')
-    print(tms)
-
-    # check results
-    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
-    tests = CTests()
-    if Re == 1e7 and M_inf == 0.2 and alpha == 15*math.pi/180:
-      tests.add(CTest('Cl', isolver.Cl, 1.615, 5e-2)) # Xfoil 1.6654
-      tests.add(CTest('Cd', vsolver.Cd, 0.0143, 1e-3, forceabs=True)) # XFOIL 0.01460
-      tests.add(CTest('Cdp', vsolver.Cdp, 0.0103, 1e-3, forceabs=True)) # XFOIL 0.01255
-      tests.add(CTest('xtr_top', vsolver.xtr[0], 0.00660, 0.003, forceabs=True)) # XFOIL 0.0067
-      tests.add(CTest('xtr_bot', vsolver.xtr[1], 1, 0.03, forceabs=True)) # XFOIL 1.0000
-    elif Re == 5e6 and M_inf == 0.2 and alpha == 15*math.pi/180:
-      tests.add(CTest('Cl', isolver.Cl, 1.58, 5e-2)) # Xfoil 1.6109
-      tests.add(CTest('Cd', vsolver.Cd, 0.0201, 1e-3, forceabs=True)) # XFOIL 0.01944
-      tests.add(CTest('Cdp', vsolver.Cdp, 0.0130, 1e-3, forceabs=True)) # XFOIL 0.01544
-      tests.add(CTest('xtr_top', vsolver.xtr[0], 0.008, 0.003, forceabs=True)) # XFOIL 0.0081
-      tests.add(CTest('xtr_bot', vsolver.xtr[1], 1, 0.03, forceabs=True)) # XFOIL 1.0000
-    else:
-        raise Exception('Test not defined for this flow')
-
-    tests.run()
-
-    # visualize solution and plot results
-    floD.initViewer(pbl)
-    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
-
-    val=validation.Validation()
-    val.main(isolver.Cl,vsolver.wData)
-
-    # eof
-    print('')
-
-if __name__ == "__main__":
-    main()
diff --git a/dart/tests/bliTransonic.py b/dart/tests/bliTransonic.py
deleted file mode 100755
index 18ad37fa44b4f2b213b467697b67e4454710d297..0000000000000000000000000000000000000000
--- a/dart/tests/bliTransonic.py
+++ /dev/null
@@ -1,171 +0,0 @@
-#!/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) viscous flow around a NACA 0012
-# Amaury Bilocq
-#
-# Test the viscous-inviscid interaction scheme
-# Reference to the master's thesis: http://hdl.handle.net/2268/252195
-# Reference test cases with Naca0012 (different from master's thesis):
-# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
-#    -> msLe = 0.001, xtrTop = 0.061
-# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
-#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
-# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
-#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-import dart.utils as floU
-import dart.default as floD
-
-import dart.viscous.Master.DCoupler as floVC
-import dart.viscous.Drivers.DDriver as floVSMaster
-
-import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
-import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
-
-import dart.viscous.Physics.DValidation as validation
-
-import tbox
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-import numpy as np
-
-import math
-
-def main():
-    # timer
-    tms = fwk.Timers()
-    tms['total'].start()
-
-    # define flow variables
-    Re = 1e7
-    alpha = 2*math.pi/180
-    M_inf = 0.715
-
-    user_xtr=[None, None]
-    user_Ti=None
-
-    mapMesh = 1
-    nFactor = 2
-
-    # Time solver parameters
-    Time_Domain = True # True for unsteady solver, False for steady solver
-    CFL0 = 1
-
-    # ========================================================================================
-
-    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
-    c_ref = 1
-    dim = 2
-    tol = 1e-4 # tolerance for the VII (usually one drag count)
-    case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
-
-    # 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.001, 'msLe' : 0.0005}
-    msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream'])
-    tms['msh'].stop()
-
-    # set the problem and add medium, initial/boundary conditions
-    tms['pre'].start()
-    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
-    tms['pre'].stop()
-
-    # solve viscous problem
-
-    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
-    tms['solver'].start()
-    isolver = floD.newton(pbl)
-
-    # Choose between steady and unsteady solver
-    if Time_Domain is True: # Run unsteady solver
-
-        vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
-        coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
-
-    elif Time_Domain is False: # Run steady solver 
-
-        vsolver=floVS_Std.Solver(Re)
-        coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
-        
-    coupler.run()
-    tms['solver'].stop()
-
-    # extract Cp
-    Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
-    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
-    # display results
-    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
-    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
-    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
-
-    # display timers
-    tms['total'].stop()
-    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
-    print('CPU statistics')
-    print(tms)
-    
-    # check results
-    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
-    try:
-      tests = CTests()
-      if Re == 1e7 and M_inf == 0.715 and alpha == 2*math.pi/180:
-          tests.add(CTest('Cl', isolver.Cl, 0.72, 5e-2))
-          tests.add(CTest('Cd', vsolver.Cd, 0.0061, 1e-3, forceabs=True))
-          tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True))
-          tests.add(CTest('xtr_top', vsolver.xtr[0], 0.34, 0.08, forceabs=True))
-          tests.add(CTest('xtr_bot', vsolver.xtr[1],  0.5, 0.03, forceabs=True))
-      elif Re == 1e7 and M_inf == 0.72 and alpha == 2*math.pi/180:
-          tests.add(CTest('Cl', isolver.Cl, 0.674, 5e-2))
-          tests.add(CTest('Cd', vsolver.Cd, 0.0061, 1e-3, forceabs=True))
-          tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True))
-          tests.add(CTest('xtr_top', vsolver.xtr[0], 0.36, 0.08, forceabs=True))
-          tests.add(CTest('xtr_bot', vsolver.xtr[1],  0.5, 0.03, forceabs=True))
-      elif Re == 1e7 and M_inf == 0.73 and alpha == 2*math.pi/180 and user_xtr==[0, 0]:
-          tests.add(CTest('Cl', isolver.Cl, 0.671, 5e-2))
-          tests.add(CTest('Cd', vsolver.Cd, 0.0089, 1e-3, forceabs=True))
-          tests.add(CTest('Cdp', vsolver.Cdp, 0.0024, 1e-3, forceabs=True))
-          tests.add(CTest('xtr_top', vsolver.xtr[0], 0, 0.08, forceabs=True))
-          tests.add(CTest('xtr_bot', vsolver.xtr[1],  0, 0.03, forceabs=True))
-      else:
-          raise Exception('Test not defined for this flow')
-      tests.run()
-    except Exception as e:
-      print(e)
-
-    # visualize solution and plot results
-    floD.initViewer(pbl)
-    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
-
-    val=validation.Validation()
-    val.main(isolver.Cl,vsolver.wData)
-
-    # eof
-    print('')
-
-if __name__ == "__main__":
-    main()
diff --git a/dart/tests/bli_Polar.py b/dart/tests/bli_Polar.py
new file mode 100644
index 0000000000000000000000000000000000000000..45330d4eb301b60cbafb132352e7a1570f1949a8
--- /dev/null
+++ b/dart/tests/bli_Polar.py
@@ -0,0 +1,124 @@
+#!/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) viscous flow around a NACA 0012
+# Amaury Bilocq
+#
+# Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012 (different from master's thesis):
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
+#    -> msLe = 0.001, xtrTop = 0.061
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
+#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import dart.utils as floU
+import dart.default as floD
+import dart.viscUtils as viscU
+
+import dart.vCoupler as floVC
+import numpy as np
+
+import tbox
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def main():
+    print('Start')
+    # timer
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    # define flow variables
+    Re = 1e7
+    alpha=0
+    alphaSeq = np.arange(-5,5.5,0.5)
+    M_inf = 0.
+
+    meshFactor = 2
+    CFL0 = 1
+
+    # ========================================================================================
+
+    c_ref = 1
+    dim = 2
+
+    # mesh the geometry
+    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+    tms['msh'].start()
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.005, 'msLe' : 0.001}
+    msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
+    tms['msh'].stop()
+
+    # set the problem and add medium, initial/boundary conditions
+    tms['pre'].start()
+    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+    tms['pre'].stop()
+
+    # solve viscous problem
+
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    
+    isolver = floD.newton(pbl)
+    vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor)
+
+    coupler = floVC.Coupler(isolver, vsolver)
+
+    coupler.CreateProblem(blw[0], blw[1])
+
+    polar = coupler.RunPolar(alphaSeq)
+    #coupler.Run()
+    tms['solver'].stop()
+
+    # display timers
+    tms['total'].stop()
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    print('SOLVERS statistics')
+    print(coupler.tms)
+
+    viscU.PlotPolar(polar)
+    
+    # check results
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    tests.add(CTest('Cl1', polar[1,0], -0.430901, 5e-3))
+    tests.add(CTest('Cd1', polar[2,0], 0.005728, 1e-3, forceabs=True))
+    tests.add(CTest('Cl2', polar[1,-1], 0.430797, 5e-3))
+    tests.add(CTest('Cd2', polar[2,-1], 0.005789, 1e-3, forceabs=True))
+
+    tests.run()
+
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()
diff --git a/dart/tests/blicppPolar.py b/dart/tests/bli_Sep.py
similarity index 82%
rename from dart/tests/blicppPolar.py
rename to dart/tests/bli_Sep.py
index 63a5ebbfaab8dac39a97c49692f33b8d7e2d64a6..4477a8b355316d49c9ffc7abdadf88f90a7b3b67 100644
--- a/dart/tests/blicppPolar.py
+++ b/dart/tests/bli_Sep.py
@@ -39,6 +39,7 @@ import dart.default as floD
 import dart.viscUtils as viscU
 
 import dart.vCoupler as floVC
+import numpy as np
 
 import tbox
 import tbox.utils as tboxU
@@ -46,8 +47,6 @@ import fwk
 from fwk.testing import *
 from fwk.coloring import ccolors
 
-import numpy as np
-
 import math
 
 def main():
@@ -58,13 +57,14 @@ def main():
 
     # define flow variables
     Re = 1e7
-    alpha = np.linspace(-5, 5, 1)*math.pi/180
-    M_inf = 0.
+    alpha = 15*math.pi/180
+    #alphaSeq = np.arange(-5, 5.5, 0.5)
+    M_inf = 0.2
 
-    meshFactor = 1
+    meshFactor = 2
     CFL0 = 1
 
-    plotVar = 'Cf'
+    plotVar = 'H'
 
     # ========================================================================================
 
@@ -79,9 +79,8 @@ def main():
     tms['msh'].stop()
 
     # set the problem and add medium, initial/boundary conditions
-
     tms['pre'].start()
-    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha[0], 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
     tms['pre'].stop()
 
     # solve viscous problem
@@ -92,10 +91,11 @@ def main():
     isolver = floD.newton(pbl)
     vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor)
 
-    coupler = floVC.Coupler(isolver, vsolver, alpha)
+    coupler = floVC.Coupler(isolver, vsolver)
 
     coupler.CreateProblem(blw[0], blw[1])
 
+    #coupler.RunPolar(alphaSeq)
     coupler.Run()
     tms['solver'].stop()
 
@@ -113,27 +113,28 @@ def main():
     print('CPU statistics')
     print(tms)
     
+    xtr=vsolver.Getxtr()
     
     
-    # visualize solution and plot results
+    """# visualize solution and plot results
     #floD.initViewer(pbl)
     tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True)
 
     blScalars, blVectors = viscU.ExtractVars(vsolver)
     xtr=vsolver.Getxtr()
     wData = viscU.Dictionary(blScalars, blVectors, xtr)
-    viscU.PlotVars(wData, plotVar)
+    viscU.PlotVars(wData, plotVar)"""
 
     
     # check results
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
-    if Re == 1e7 and M_inf == 0.2 and alpha == 5*math.pi/180:
-        tests.add(CTest('Cl', isolver.Cl, 0.56, 5e-2)) # XFOIL 0.58
-        tests.add(CTest('Cd', vsolver.Cdt, 0.0063, 1e-3, forceabs=True)) # XFOIL 0.00617
-        tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True)) # XFOIL 0.00176
-        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.059, 0.03, forceabs=True)) # XFOIL 0.0510
-        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.738, 0.03, forceabs=True)) # XFOIL 0.7473
+    if Re == 1e7 and M_inf == 0.2 and alpha == 15*math.pi/180:
+        tests.add(CTest('Cl', isolver.Cl, 1.61416, 5e-2))
+        tests.add(CTest('Cd', vsolver.Cdt, 0.01452, 1e-3, forceabs=True))
+        tests.add(CTest('Cdp', vsolver.Cdp, 0.0105, 1e-3, forceabs=True))
+        tests.add(CTest('xtr_top', xtr[0], 0.01, 0.03, forceabs=True))
+        tests.add(CTest('xtr_bot', xtr[1], 0.99, 0.03, forceabs=True))
     else:
         raise Exception('Test not defined for this flow')
 
diff --git a/dart/tests/blicpp.py b/dart/tests/bli_Trans.py
similarity index 92%
rename from dart/tests/blicpp.py
rename to dart/tests/bli_Trans.py
index a2eefe0be72b588f15ecf2004dab1b75f08ed966..21ceb6fd9340b770b7c80f2c02d8b6248b502c8a 100644
--- a/dart/tests/blicpp.py
+++ b/dart/tests/bli_Trans.py
@@ -57,14 +57,14 @@ def main():
 
     # define flow variables
     Re = 1e7
-    alpha = 2*math.pi/180
-    alphaSeq = np.linspace(0, 10, 10)
-    M_inf = 0.
+    alpha = 2.*math.pi/180
+    #alphaSeq = np.arange(-5, 5.5, 0.5)
+    M_inf = 0.715
 
-    meshFactor = 1
+    meshFactor = 2
     CFL0 = 1
 
-    plotVar = 'Cf'
+    plotVar = 'H'
 
     # ========================================================================================
 
@@ -74,8 +74,8 @@ 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}
-    msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.005}
+    msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
     # set the problem and add medium, initial/boundary conditions
@@ -95,7 +95,9 @@ def main():
 
     coupler.CreateProblem(blw[0], blw[1])
 
-    coupler.RunPolar(alphaSeq)
+    #coupler.RunPolar(alphaSeq)
+    coupler.resetInv = True
+    coupler.Run()
     tms['solver'].stop()
 
     # extract Cp
@@ -114,14 +116,14 @@ def main():
     
     
     
-    # visualize solution and plot results
+    """# visualize solution and plot results
     #floD.initViewer(pbl)
     tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True)
 
     blScalars, blVectors = viscU.ExtractVars(vsolver)
     xtr=vsolver.Getxtr()
     wData = viscU.Dictionary(blScalars, blVectors, xtr)
-    viscU.PlotVars(wData, plotVar)
+    viscU.PlotVars(wData, plotVar)"""
 
     
     # check results
diff --git a/dart/tests/bliHighLift.py b/dart/tests/bli_dualMesh.py
old mode 100755
new mode 100644
similarity index 65%
rename from dart/tests/bliHighLift.py
rename to dart/tests/bli_dualMesh.py
index 022ee4b8991cabf951e915eff74b94b73be394d5..20f135802802efe5e4e0ea3983ffcf80214ad65b
--- a/dart/tests/bliHighLift.py
+++ b/dart/tests/bli_dualMesh.py
@@ -36,14 +36,10 @@
 
 import dart.utils as floU
 import dart.default as floD
+import dart.viscUtils as viscU
 
-import dart.viscous.Master.DCoupler as floVC
-import dart.viscous.Drivers.DDriver as floVSMaster
-
-import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
-import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
-
-import dart.viscous.Physics.DValidation as validation
+import dart.vCoupler as floVC
+import numpy as np
 
 import tbox
 import tbox.utils as tboxU
@@ -54,32 +50,26 @@ from fwk.coloring import ccolors
 import math
 
 def main():
+    print('Start')
     # timer
     tms = fwk.Timers()
     tms['total'].start()
 
     # define flow variables
-    Re = 3e6
-    alpha = 10.*math.pi/180
-    M_inf = 0.15
-    user_xtr=[0, 0]
-    user_Ti=None
-
-    mapMesh = 1
-    nFactor = 3
+    Re = 1e7
+    alpha = 2*math.pi/180
+    #alphaSeq = np.arange(-5, 5.5, 0.5)
+    M_inf = 0.
 
-    # Time solver parameters
-    Time_Domain=True # True for unsteady solver, False for steady solver
+    meshFactor = 2
     CFL0 = 1
 
+    plotVar = 'H'
+
     # ========================================================================================
 
-    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
     c_ref = 1
     dim = 2
-    tol = 1e-4 # tolerance for the VII (usually one drag count)
-    case='case15Xfoil.dat' # File containing results from XFOIL of the considered case
-
 
     # mesh the geometry
     print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
@@ -94,19 +84,19 @@ def main():
     tms['pre'].stop()
 
     # solve viscous problem
+
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
     tms['solver'].start()
+    
     isolver = floD.newton(pbl)
+    vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor)
+
+    coupler = floVC.Coupler(isolver, vsolver)
 
-  # Choose between steady and unsteady solver
-    if Time_Domain is True: # Run unsteady solver
-        vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
-        coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
+    coupler.CreateProblem(blw[0], blw[1])
 
-    elif Time_Domain is False: # Run steady solver 
-        vsolver=floVS_Std.Solver(Re)
-        coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
-    coupler.run()
+    #coupler.RunPolar(alphaSeq)
+    coupler.Run()
     tms['solver'].stop()
 
     # extract Cp
@@ -115,34 +105,41 @@ def main():
     # display results
     print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
     print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
-    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
 
     # display timers
     tms['total'].stop()
     print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
     print('CPU statistics')
     print(tms)
-
+    
+    xtr=vsolver.Getxtr()
+    
+    
+    """# visualize solution and plot results
+    #floD.initViewer(pbl)
+    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True)
+
+    blScalars, blVectors = viscU.ExtractVars(vsolver)
+    xtr=vsolver.Getxtr()
+    wData = viscU.Dictionary(blScalars, blVectors, xtr)
+    viscU.PlotVars(wData, plotVar)"""
+
+    
     # check results
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
-    if Re==3e6 and M_inf==0.15 and alpha==10*math.pi/180 and user_xtr==[0., 0.]:
-        tests.add(CTest('Cl', isolver.Cl, 1.08, 5e-2)) # RANS 1.1 Xfoil 1.077
-        tests.add(CTest('Cd', vsolver.Cd, 0.0134, 1e-3, forceabs=True)) # RANS 0.0132 XFOIL 0.0128
-        tests.add(CTest('Cdp', vsolver.Cdp, 0.0058, 1e-3, forceabs=True))
-        tests.add(CTest('xtr_top', vsolver.xtr[0], 0, 1e-3, forceabs=True))
-        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0, 1e-3, forceabs=True))
+    if Re == 1e7 and M_inf == 0. and alpha == 2.*math.pi/180:
+        tests.add(CTest('Cl', isolver.Cl, 0.22113, 5e-3))
+        tests.add(CTest('Cd', vsolver.Cdt, 0.00534, 1e-3, forceabs=True))
+        tests.add(CTest('Cdp', vsolver.Cdp, 0.0009, 1e-3, forceabs=True))
+        tests.add(CTest('xtr_top', xtr[0], 0.20, 0.03, forceabs=True))
+        tests.add(CTest('xtr_bot', xtr[1], 0.50, 0.03, forceabs=True))
     else:
         raise Exception('Test not defined for this flow')
 
     tests.run()
 
-    # visualize solution and plot results
-    floD.initViewer(pbl)
-    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
-
-    val=validation.Validation()
-    val.main(isolver.Cl,vsolver.wData)
 
     # eof
     print('')
diff --git a/dart/tests/bli.py b/dart/tests/bli_lowLift.py
similarity index 64%
rename from dart/tests/bli.py
rename to dart/tests/bli_lowLift.py
index fd2f186f340c56719e664fcdca3dc8edb89fed6b..01084014b245d91a42d24910eeb3eaf51e6f48f5 100644
--- a/dart/tests/bli.py
+++ b/dart/tests/bli_lowLift.py
@@ -34,18 +34,23 @@
 # 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.utils as floU
 import dart.default as floD
-import dart.viscous.solver as floVS
-import dart.viscous.coupler as floVC
+import dart.viscUtils as viscU
+
+import dart.vCoupler as floVC
+import numpy as np
+
 import tbox
 import tbox.utils as tboxU
 import fwk
 from fwk.testing import *
 from fwk.coloring import ccolors
 
+import math
+
 def main():
+    print('Start')
     # timer
     tms = fwk.Timers()
     tms['total'].start()
@@ -53,31 +58,45 @@ def main():
     # define flow variables
     Re = 1e7
     alpha = 5*math.pi/180
-    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
-    M_inf = 0.5
+    #alphaSeq = np.arange(-5, 5.5, 0.5)
+    M_inf = 0.
+
+    meshFactor = 1
+    CFL0 = 5
+
+    plotVar = 'H'
+
+    # ========================================================================================
+
     c_ref = 1
     dim = 2
-    tol = 1e-4 # tolerance for the VII (usually one drag count)
 
     # 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.01}
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
     msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
-    # set the problem and add fluid, initial/boundary conditions
+    # set the problem and add medium, initial/boundary conditions
     tms['pre'].start()
     pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
     tms['pre'].stop()
 
     # solve viscous problem
+
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
     tms['solver'].start()
+    
     isolver = floD.newton(pbl)
-    vsolver = floVS.Solver(Re)
-    coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
-    coupler.run()
+    vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor)
+
+    coupler = floVC.Coupler(isolver, vsolver)
+
+    coupler.CreateProblem(blw[0], blw[1])
+
+    #coupler.RunPolar(alphaSeq)
+    coupler.Run()
     tms['solver'].stop()
 
     # extract Cp
@@ -86,44 +105,48 @@ def main():
     # display results
     print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
     print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
-    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
 
     # display timers
     tms['total'].stop()
     print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
     print('CPU statistics')
     print(tms)
+    print('SOLVERS statistics')
+    print(coupler.tms)
+
+    xtr=vsolver.Getxtr()
+    
+    """# visualize solution and plot results
+    #floD.initViewer(pbl)
+    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True)
 
-    # visualize solution and plot results
-    floD.initViewer(pbl)
-    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
+    blScalars, blVectors = viscU.ExtractVars(vsolver)
+    wData = viscU.Dictionary(blScalars, blVectors, xtr)
+    viscU.PlotVars(wData, plotVar)"""
 
+    
     # check results
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
-    if Re == 1e7 and M_inf == 0 and alpha == 5*math.pi/180:
-        tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.58
-        tests.add(CTest('Cd', vsolver.Cd, 0.0062, 1e-3, forceabs=True))
-        tests.add(CTest('Cdp', vsolver.Cdp, 0.0018, 1e-3, forceabs=True))
-        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.061, 0.03, forceabs=True)) # Xfoil 0.056
-        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.740, 0.03, forceabs=True))
-    elif Re == 1e7 and M_inf == 0.5 and alpha == 5*math.pi/180:
-        tests.add(CTest('Cl', isolver.Cl, 0.65, 5e-2)) # Xfoil 0.69
-        tests.add(CTest('Cd', vsolver.Cd, 0.0067, 1e-3, forceabs=True))
-        tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 1e-3, forceabs=True))
-        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.049, 0.03, forceabs=True)) # Xfoil 0.038
-        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.736, 0.03, forceabs=True))
-    elif Re == 1e7 and M_inf == 0 and alpha == 12*math.pi/180:
-        tests.add(CTest('Cl', isolver.Cl, 1.30, 5e-2)) # Xfoil 1.39
-        tests.add(CTest('Cd', vsolver.Cd, 0.011, 1e-3, forceabs=True))
-        tests.add(CTest('Cdp', vsolver.Cdp, 0.0064, 1e-3, forceabs=True))
-        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.010, 0.03, forceabs=True)) # Xfoil 0.008
-        tests.add(CTest('xtr_bot', vsolver.xtr[1], 1.000, 0.03, forceabs=True))
+    if Re == 1e7 and M_inf == 0. and alpha == 2.*math.pi/180:
+        tests.add(CTest('Cl', isolver.Cl, 0.2208, 5e-3))
+        tests.add(CTest('Cd', vsolver.Cdt, 0.00531, 1e-3, forceabs=True))
+        tests.add(CTest('Cdp', vsolver.Cdp, 0.0009, 1e-3, forceabs=True))
+        tests.add(CTest('xtr_top', xtr[0], 0.20, 0.03, forceabs=True))
+        tests.add(CTest('xtr_bot', xtr[1], 0.50, 0.03, forceabs=True))
+    if Re == 1e7 and M_inf == 0. and alpha == 5.*math.pi/180:
+        tests.add(CTest('Cl', isolver.Cl, 0.5488, 5e-3))
+        tests.add(CTest('Cd', vsolver.Cdt, 0.0062, 1e-3, forceabs=True))
+        tests.add(CTest('Cdp', vsolver.Cdp,0.0018, 1e-3, forceabs=True))
+        tests.add(CTest('xtr_top', xtr[0], 0.06, 0.03, forceabs=True))
+        tests.add(CTest('xtr_bot', xtr[1], 0.74, 0.03, forceabs=True))
     else:
         raise Exception('Test not defined for this flow')
 
     tests.run()
 
+
     # eof
     print('')
 
diff --git a/dart/tests/bliNonLift.py b/dart/tests/bli_python.py
similarity index 89%
rename from dart/tests/bliNonLift.py
rename to dart/tests/bli_python.py
index 27b26d298646dc42ca10bf4d11fcfe9be557a34e..ca729b16b449ff9e27ee27bac861781c2e107435 100755
--- a/dart/tests/bliNonLift.py
+++ b/dart/tests/bli_python.py
@@ -60,7 +60,7 @@ def main():
 
     # define flow variables
     Re = 1e7
-    alpha = 0.*math.pi/180
+    alpha = 5.*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 = 1
+    mapMesh = 0
     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.01}
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
     #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()
@@ -126,16 +126,18 @@ def main():
     print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
     print('CPU statistics')
     print(tms)
+    print('SOLVERS statistics')
+    print(coupler.tms)
     
     # check results
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
-    if Re == 1e7 and M_inf == 0.2 and alpha == 5*math.pi/180:
-        tests.add(CTest('Cl', isolver.Cl, 0.56, 5e-2)) # XFOIL 0.58
-        tests.add(CTest('Cd', vsolver.Cd, 0.0063, 1e-3, forceabs=True)) # XFOIL 0.00617
-        tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True)) # XFOIL 0.00176
-        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.059, 0.03, forceabs=True)) # XFOIL 0.0510
-        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.738, 0.03, forceabs=True)) # XFOIL 0.7473
+    if Re == 1e7 and M_inf == 0. and alpha == 2*math.pi/180:
+        tests.add(CTest('Cl', isolver.Cl, 0.2197, 5e-2)) # XFOIL 0.58
+        tests.add(CTest('Cd', vsolver.Cd, 0.00529, 1e-3, forceabs=True)) # XFOIL 0.00617
+        tests.add(CTest('Cdp', vsolver.Cdp, 0.0008, 1e-3, forceabs=True)) # XFOIL 0.00176
+        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.20395, 0.03, forceabs=True)) # XFOIL 0.0510
+        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.49996, 0.03, forceabs=True)) # XFOIL 0.7473
     else:
         raise Exception('Test not defined for this flow')
 
diff --git a/dart/vCoupler.py b/dart/vCoupler.py
index 6224e6e078b8df1094e5f38468a934af5e31c1de..9e721714c9adad9b9173aa138c9d9e1f67912a8c 100644
--- a/dart/vCoupler.py
+++ b/dart/vCoupler.py
@@ -19,7 +19,6 @@
 ## VII Coupler
 # Paul Dechamps
 
-from multiprocessing import set_forkserver_preload
 from dart.viscous.Master.DAirfoil import Airfoil
 from dart.viscous.Master.DWake import Wake
 from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
@@ -33,208 +32,179 @@ 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 = 150
-        self.couplTol = 1e-4
-
-        self.tms=fwk.Timers()
-        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):
-
-      # Initilize meshes in c++ objects
-
-      nElmUpperPrev = 0
-      couplIter = 0
-      cdPrev = 0
-      print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'
-      .format(self.iSolver.pbl.alpha*180/math.pi, self.iSolver.pbl.M_inf, self.vSolver.Re/1e6))
-      print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
-      print(' --------------------------------------------------------')
-      while couplIter < self.maxCouplIter:
-
-        # Run inviscid solver
-        self.tms['inv'].start()
-        self.iSolver.run()
-        self.tms['inv'].stop()
-
-        # Extract inviscid boundary
-
-        for n in range(0, len(self.group)):
-
-          for i in range(0, len(self.group[n].boundary.nodes)):
-
-            self.group[n].v[i,0] = self.iSolver.U[self.group[n].boundary.nodes[i].row][0]
-            self.group[n].v[i,1] = self.iSolver.U[self.group[n].boundary.nodes[i].row][1]
-            self.group[n].v[i,2] = self.iSolver.U[self.group[n].boundary.nodes[i].row][2]
-            self.group[n].Me[i] = self.iSolver.M[self.group[n].boundary.nodes[i].row]
-            self.group[n].rhoe[i] = self.iSolver.rho[self.group[n].boundary.nodes[i].row]
-        
-        dataIn=[None,None]
-        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 ((len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]]) != nElmUpperPrev) or couplIter==0):
-
-          # print('STAG POINT MOVED')
-          
-          # 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(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]:])
-        # print('Inviscid boundary imposed')
-        # Run viscous solver
-
-        self.tms['visc'].start()
-        exitCode = self.vSolver.Run(couplIter)
-        self.tms['visc'].stop()
-        # 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.
-        if couplIter >= 12:
-          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:
-
-          #print(self.tms)
-          print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}\n'.format(couplIter, self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)),ccolors.ANSI_RESET)
-
-          """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)
-          print('')"""
-          return 0
-        
-        cdPrev = self.vSolver.Cdt
+  def __init__(self, iSolver, vSolver) -> None:
+    self.iSolver=iSolver
+    self.vSolver=vSolver
+
+    self.maxCouplIter = 150
+    self.couplTol = 1e-4
+
+    self.tms=fwk.Timers()
+    self.resetInv = False
+    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):
+
+    # Initilize meshes in c++ objects
+
+    self.nElmUpperPrev = 0
+    couplIter = 0
+    cdPrev = 0
+    print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'
+    .format(self.iSolver.pbl.alpha*180/math.pi, self.iSolver.pbl.M_inf, self.vSolver.Re/1e6))
+    print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
+    print(' --------------------------------------------------------')
+    while couplIter < self.maxCouplIter:
+
+      # Run inviscid solver
+      self.tms['inviscid'].start()
+      if self.resetInv:
+        self.iSolver.reset()
+      self.iSolver.run()
+      self.tms['inviscid'].stop()
+
+      # Extract inviscid boundary
+      self.tms['processing'].start()
+      self.__ImposeInvBC(couplIter)
+      self.tms['processing'].stop()
+
+      # Run viscous solver
+      self.tms['viscous'].start()
+      exitCode = self.vSolver.Run(couplIter)
+      self.tms['viscous'].stop()
+      # print('Viscous ops terminated with exit code ', exitCode)
+
+      self.tms['processing'].start()
+      self.__ImposeBlwVel()
+      self.tms['processing'].stop()
+
+      error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
+      
+      if error  <= self.couplTol and exitCode==0:
+
+        print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}\n'.format(couplIter,
+        self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)),ccolors.ANSI_RESET)
+        return 0
+      
+      cdPrev = self.vSolver.Cdt
 
-        xtr=self.vSolver.Getxtr()
 
-      #print('{:>5s} {:>12s} {:>12s} {:>12s} {:>6s} {:>12s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
-        if couplIter%5==0:
-          print(' {:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}'.format(couplIter, self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)))
-        
-        """print('')
-        print(ccolors.ANSI_BLUE, '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)
-        print('')"""
-
-        # Retreive and set blowing velocities.
-
-        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))
-
-        """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)
+      if couplIter%5==0:
         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
+        print(' {:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}'.format(couplIter,
+        self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)))
 
-        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()
+      couplIter += 1
+      self.tms['processing'].stop()
+    else:
+      print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
+      return couplIter
 
-        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()"""
+  def RunPolar(self, alphaSeq):
+    
+    import math
 
-        self.group[0].u = blwVelAF[self.group[0].connectListElems.argsort()]
-        self.group[1].u = blwVelWK[self.group[1].connectListElems.argsort()]
+    self.iSolver.printOn = False
+    self.vSolver.printOn = False
 
-        self.group[0].deltaStar = deltaStarAF[self.group[0].connectListNodes.argsort()]
-        self.group[1].deltaStar = deltaStarWK[self.group[1].connectListNodes.argsort()]
+    coeffTable = np.empty((3,len(alphaSeq)))
+    coeffTable[0,:] = alphaSeq
 
-        self.group[0].xx = xxAF[self.group[0].connectListNodes.argsort()]
-        self.group[1].xx = xxWK[self.group[1].connectListNodes.argsort()]
+    for ind, alpha in enumerate(alphaSeq):
+      self.iSolver.pbl.update(alpha*math.pi/180)
 
-        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])
+      self.Run()
 
-        couplIter += 1
-      else:
-        print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
-        return couplIter
+      # Collect Coeff
 
+      coeffTable[1,ind] = self.iSolver.Cl
+      coeffTable[2,ind] = self.vSolver.Cdt
 
-    def RunPolar(self, alphaSeq):
-      
-      import math
+    return coeffTable
+
+  def __ImposeBlwVel(self):
+
+    # Retreive and set blowing velocities.
 
-      self.iSolver.printOn = False
-      self.vSolver.printOn = False
+    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))
 
-      coeffTable = np.empty((2,len(alphaSeq)))
+    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]
 
-      for ind, alpha in enumerate(alphaSeq):
-        self.iSolver.pbl.update(alpha*math.pi/180)
+    self.group[0].u = blwVelAF[self.group[0].connectListElems.argsort()]
+    self.group[1].u = blwVelWK[self.group[1].connectListElems.argsort()]
 
-        self.Run()
+    self.group[0].deltaStar = deltaStarAF[self.group[0].connectListNodes.argsort()]
+    self.group[1].deltaStar = deltaStarWK[self.group[1].connectListNodes.argsort()]
 
-        # Collect Coeff
+    self.group[0].xx = xxAF[self.group[0].connectListNodes.argsort()]
+    self.group[1].xx = xxWK[self.group[1].connectListNodes.argsort()]
 
-        coeffTable[0,ind] = self.iSolver.Cl
-        coeffTable[1,ind] = self.vSolver.Cdt
+    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])
 
-      plt.plot(coeffTable[0,:], coeffTable[1,:], 'x-')
-      plt.xlabel('$c_l$ (-)')
-      plt.ylabel('$c_d$ (-)')
-      plt.show()
+  def __ImposeInvBC(self, couplIter):
+
+    for n in range(0, len(self.group)):
+
+      for i in range(0, len(self.group[n].boundary.nodes)):
+
+        self.group[n].v[i,0] = self.iSolver.U[self.group[n].boundary.nodes[i].row][0]
+        self.group[n].v[i,1] = self.iSolver.U[self.group[n].boundary.nodes[i].row][1]
+        self.group[n].v[i,2] = self.iSolver.U[self.group[n].boundary.nodes[i].row][2]
+        self.group[n].Me[i] = self.iSolver.M[self.group[n].boundary.nodes[i].row]
+        self.group[n].rhoe[i] = self.iSolver.rho[self.group[n].boundary.nodes[i].row]
+    
+    dataIn=[None,None]
+    for n in range(0,len(self.group)):
+      self.group[n].connectListNodes, self.group[n].connectListElems, dataIn[n]  = self.group[n].connectList()
+    fMeshDict, _, _ = GroupConstructor().mergeVectors(dataIn, 0, 1)
+
+    if ((len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]]) != self.nElmUpperPrev) or couplIter==0):
+
+      # print('STAG POINT MOVED')
+      
+      # 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(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]:])
+
+    self.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]:])
+    # print('Inviscid boundary imposed')
+  
+
+
+  
 
 
 
diff --git a/dart/viscUtils.py b/dart/viscUtils.py
index 45355a31391013e8e022f3f44a01bfc6d4896d22..763036be821231b502f18cf9f64cdbfc25608207 100644
--- a/dart/viscUtils.py
+++ b/dart/viscUtils.py
@@ -90,4 +90,16 @@ def PlotVars(wData, var='H'):
   plt.ylabel(var)
   plt.xlim([0,2])
   plt.legend(['Upper side', 'Lower side'], frameon=False)
+  plt.show()
+
+def PlotPolar(coeffTable):
+  from matplotlib import pyplot as plt
+
+  plt.plot(coeffTable[0,:], coeffTable[1,:], 'x-')
+  plt.xlabel('$AoA$ (deg)')
+  plt.ylabel('$c_l$ (-)')
+  plt.show()
+  plt.plot(coeffTable[1,:], coeffTable[2,:], 'x-')
+  plt.xlabel('$c_l$ (-)')
+  plt.ylabel('$c_d$ (-)')
   plt.show()
\ No newline at end of file
diff --git a/dart/viscous/Drivers/DGroupConstructor.py b/dart/viscous/Drivers/DGroupConstructor.py
index 176a67873cd83849814e28c1c479679880c53b2e..fab1f757d6968f25339b20eba95e6f76dba6dcda 100755
--- a/dart/viscous/Drivers/DGroupConstructor.py
+++ b/dart/viscous/Drivers/DGroupConstructor.py
@@ -97,7 +97,7 @@ class GroupConstructor():
         else:
 
             Mesh         = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][:,0], dataIn[1][:,0]))
-            MeshbNodes   = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][:,0])]
+            MeshbNodes   = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][:,0]) + len(dataIn[1][:,0])]
 
             xx           = np.concatenate((xx_up, xx_lw, xx_wk))
             vtExt        = np.concatenate((vtExt_up, vtExt_lw, vtExt_wk))
diff --git a/dart/viscous/Master/DCoupler.py b/dart/viscous/Master/DCoupler.py
index 295cfc5be8c0fd36f36b82748101a09511ea06ce..eea34c41eb4e9dc6590bce8bede32035da54a897 100755
--- a/dart/viscous/Master/DCoupler.py
+++ b/dart/viscous/Master/DCoupler.py
@@ -66,6 +66,8 @@ class Coupler:
             self.isolver.run()  
             self.tms['inv'].stop()
 
+            self.tms['processing'].start()
+
             for n in range(0, len(self.group)):
 
                 for i in range(0, len(self.group[n].boundary.nodes)):
@@ -82,11 +84,15 @@ class Coupler:
                 Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)
                 tboxU.write(Cp, 'Cpinv_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
 
+            self.tms['processing'].stop()
+
             # Run viscous solver
             print('+------------------------------ Viscous solver --------------------------------+')
             self.tms['visc'].start()
             self.vsolver.run(self.group)
             self.tms['visc'].stop()
+              
+            self.tms['processing'].start()
 
             for n in range(0, len(self.group)):
                 for i in range(0, self.group[n].nE):
@@ -97,9 +103,7 @@ class Coupler:
 
             # Check convergence and output coupling iteration.
 
-            print('tms @ it', it, self.tms)
-            if it == 84:
-              converged = 1
+            #print('tms @ it', it, self.tms)
 
             if error <= self.tol:
               print(self.tms)
@@ -125,11 +129,13 @@ class Coupler:
             # Prepare next iteration.
             it += 1
             self.vsolver.it += 1
+            self.tms['processing'].stop()
+
 
         # Save results
 
-        self.isolver.save(self.writer)    # Inviscid solver files.
+        """self.isolver.save(self.writer)    # Inviscid solver files.
         self.vsolver.writeFile()          # Viscous solver files.
         Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)    # Pressure coefficient.
         tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
-        print('\n')
+        print('\n')"""