diff --git a/dart/pyVII/vCoupler.py b/dart/pyVII/vCoupler.py
index 7f0b8c97109498a8534218903db521294ef4ff17..79f4459069fb17b129aba63c8c5a5af758a87e3e 100644
--- a/dart/pyVII/vCoupler.py
+++ b/dart/pyVII/vCoupler.py
@@ -46,6 +46,8 @@ class Coupler:
 
     # Initilize meshes in c++ objects
 
+    dragIter = []
+
     self.nElmUpperPrev = 0
     couplIter = 0
     cdPrev = 0
@@ -77,20 +79,21 @@ class Coupler:
       self.iSolverAPI.SetBlowingVel(self.vSolver)
       self.tms['processing'].stop()
 
+      dragIter.append(self.vSolver.Cdt)
       error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
       
       if error  <= self.couplTol:
 
-        print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}\n'.format(couplIter,
+        print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.3f} {:>6.3f} {:>8.2f}\n'.format(couplIter,
         self.iSolverAPI.GetCl(), self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)),ccolors.ANSI_RESET)
         self.tms['tot'].stop()
         
-        return 0
+        return dragIter
       
       cdPrev = self.vSolver.Cdt
 
 
-      print('- {:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}'.format(couplIter,
+      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)))
 
       couplIter += 1
@@ -101,7 +104,7 @@ class Coupler:
       print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
       self.tms['tot'].stop()
 
-      return couplIter
+      return dragIter
 
 
   def RunPolar(self, alphaSeq):
diff --git a/dart/pyVII/vInterpolator.py b/dart/pyVII/vInterpolator.py
index 9dc992089c4ef81e4151820d9a420f20dd23d637..b259a98e00d68d58d69a21b27ced9acd6ed9da36 100644
--- a/dart/pyVII/vInterpolator.py
+++ b/dart/pyVII/vInterpolator.py
@@ -1,7 +1,6 @@
 import fwk
 
 import numpy as np
-import math as m
 from matplotlib import pyplot as plt
 from scipy.interpolate import bisplrep
 from scipy.interpolate import bisplev
@@ -10,7 +9,7 @@ from scipy.interpolate import RBFInterpolator
 
 class Interpolator:
 
-    def __init__(self, _dartSolver, bnd_wing, bnd_wake, nSections, nDim, nPoints=[200, 25], eps=0.01, k=[5, 5]) -> None:
+    def __init__(self, _dartSolver, bnd_wing, bnd_wake, nSections, nDim, nPoints=[500, 25], eps=0.05, k=[5, 5]) -> None:
         self.iSolver = _dartSolver
 
         self.bndWing = bnd_wing
@@ -21,6 +20,7 @@ class Interpolator:
         self.Sections, self.Wing, self.Wake = self.__GetWing(nSections, nDim, nPoints, eps, k)
         self.xxPrev = [[np.zeros(len(self.Sections[iSec][iReg][:,0])) for iReg in range(len(self.Sections[iSec]))] for iSec in range(len(self.Sections))]
         self.DeltaStarPrev = [[np.zeros(len(self.Sections[iSec][iReg][:,0])) for iReg in range(len(self.Sections[iSec]))] for iSec in range(len(self.Sections))]
+        self.interpMode = 2 # 0 For Rbf, 1 for 1d interpolation.
         self.tms = fwk.Timers()
 
     def GetInvBC(self, vSolver, couplIter):
@@ -52,14 +52,31 @@ class Interpolator:
             M_wk[iNode] = self.iSolver.M[n.row]
             Rho_wk[iNode] = self.iSolver.rho[n.row]
         
-        self.tms['Rbf'].start()
-        Ux = self.__InterpToSections(U_wg[:,0], U_wk[:,0])
-        Uy = self.__InterpToSections(U_wg[:,1], U_wk[:,1])
-        Uz = self.__InterpToSections(U_wg[:,2], U_wk[:,2])
-
-        Me = self.__InterpToSections(M_wg, M_wk)
-        Rho = self.__InterpToSections(Rho_wg, Rho_wk)
-        self.tms['Rbf'].stop()
+        self.tms['interpolation'].start()
+        if self.interpMode == 0:
+            Ux = self.__RealToSections(U_wg[:,0], U_wk[:,0])
+            Uy = self.__RealToSections(U_wg[:,1], U_wk[:,1])
+            Uz = self.__RealToSections(U_wg[:,2], U_wk[:,2])
+
+            Me = self.__RealToSections(M_wg, M_wk)
+            Rho = self.__RealToSections(Rho_wg, Rho_wk)
+        elif self.interpMode == 1:
+            Ux = self.__RbfToSections(U_wg[:,0], U_wk[:,0])
+            Uy = self.__RbfToSections(U_wg[:,1], U_wk[:,1])
+            Uz = self.__RbfToSections(U_wg[:,2], U_wk[:,2])
+
+            Me = self.__RbfToSections(M_wg, M_wk)
+            Rho = self.__RbfToSections(Rho_wg, Rho_wk)
+        elif self.interpMode == 2:
+            Ux = self.__InterpToSections(U_wg[:,0], U_wk[:,0])
+            Uy = self.__InterpToSections(U_wg[:,1], U_wk[:,1])
+            Uz = self.__InterpToSections(U_wg[:,2], U_wk[:,2])
+
+            Me = self.__InterpToSections(M_wg, M_wk)
+            Rho = self.__InterpToSections(Rho_wg, Rho_wk)
+        else:
+            raise RuntimeError('Unknown interpolation mode ', self.interpMode)
+        self.tms['interpolation'].stop()
 
         # Find stagnation point
         for iSec in range(len(self.Sections)):
@@ -67,13 +84,13 @@ class Interpolator:
             idxStgUp = np.argmin(np.sqrt(Ux[iSec][0]**2+Uy[iSec][0]**2))
             idxStgLw = np.argmin(np.sqrt(Ux[iSec][1]**2+Uy[iSec][1]**2))
 
-            if np.sqrt(Ux[iSec][0][idxStgUp]**2+Uy[iSec][0][idxStgUp]**2) >= np.sqrt(Ux[iSec][1][idxStgLw]**2+Uy[iSec][1][idxStgLw]**2):
-                reg = 1
-                idx = idxStgLw
-            else:
+            if np.sqrt(Ux[iSec][0][idxStgUp]**2+Uy[iSec][0][idxStgUp]**2) <= np.sqrt(Ux[iSec][1][idxStgLw]**2+Uy[iSec][1][idxStgLw]**2):
                 reg = 0
                 idx = idxStgUp
-            
+            else:
+                reg = 1
+                idx = idxStgLw
+
             xUp, xLw = self.__Remesh(self.Sections[iSec][0][:,0], self.Sections[iSec][1][:,0], idx, reg)
             yUp, yLw = self.__Remesh(self.Sections[iSec][0][:,1], self.Sections[iSec][1][:,1], idx, reg)
             zUp, zLw = self.__Remesh(self.Sections[iSec][0][:,2], self.Sections[iSec][1][:,2], idx, reg)
@@ -117,7 +134,7 @@ class Interpolator:
                 BlowingVel[iSec][iRegion] = np.asarray(vSolver.ExtractBlowingVel(iSec, iRegion))
                 self.DeltaStarPrev[iSec][iRegion] = vSolver.ExtractDeltaStar(iSec, iRegion)
                 self.xxPrev[iSec][iRegion] = vSolver.Extractxx(iSec, iRegion)
-
+    
         if self.nDim == 2:
 
             # Compute cell centroid x positions (in the viscous mesh)
@@ -137,8 +154,8 @@ class Interpolator:
                 cumDx += (xCoord[0][1][iCell] - xCoord[0][1][iCell-1]) # Cumulative space between cells
 
             xCellsWk = np.zeros(len(xCoord[0][2])-1)
-            xCellsWk[0] = (xCoord[0][2][1] - xCoord[0][2][0])/2
-            cumDx = xCoord[0][2][1] - xCoord[0][2][0]
+            xCellsWk[0] = self.Sections[iSec][0][-1, 0] + (xCoord[0][2][1] - xCoord[0][2][0])/2
+            cumDx = self.Sections[iSec][0][-1, 0] + xCoord[0][2][1] - xCoord[0][2][0]
             for iCell in range(1, len(xCoord[0][2])-1):
                 xCellsWk[iCell] = cumDx + (xCoord[0][2][iCell+1] - xCoord[0][2][iCell])/2
                 cumDx += (xCoord[0][2][iCell] - xCoord[0][2][iCell-1]) # Cumulative space between cells
@@ -150,15 +167,32 @@ class Interpolator:
             blwSplnWk = interp1d(xCellsWk, BlowingVel[0][2], bounds_error=False, fill_value="extrapolate")
             
             # Airfoil (Cell centroids are now computed in the inviscid mesh)
-            
-            for iNode in range(len(self.bndWing.nodes)-1):
-                if self.bndWing.nodes[iNode].pos[1]>=0:
-                    self.bndWing.setU(iNode, float(blwSplnUp(self.bndWing.nodes[iNode].pos[0])))
+            for iElm, e in enumerate(self.bndWing.tag.elems):
+
+                # Compute cg of element e
+                posCg = np.zeros(3)
+                for n in e.nodes:
+                    posCg[0] += n.pos[0]
+                    posCg[1] += n.pos[1]
+                    posCg[2] += n.pos[2]
+                posCg/=len(e.nodes)
+                
+                if posCg[1]>=0:
+                    self.bndWing.setU(iElm, float(blwSplnUp(posCg[0])))
                 else:
-                    self.bndWing.setU(iNode, float(blwSplnLw(self.bndWing.nodes[iNode].pos[0])))
+                    self.bndWing.setU(iElm, float(blwSplnLw(posCg[0])))
             
-            for iNode in range(len(self.bndWake.nodes)-1):
-                self.bndWake.setU(iNode, float(blwSplnWk(self.bndWake.nodes[iNode].pos[0])))
+            for iElm, e in enumerate(self.bndWake.tag.elems):
+
+                # Compute cg of element e
+                posCg = np.zeros(3)
+                for n in e.nodes:
+                    posCg[0] += n.pos[0]
+                    posCg[1] += n.pos[1]
+                    posCg[2] += n.pos[2]
+                posCg/=len(e.nodes)
+
+                self.bndWake.setU(iElm, float(blwSplnWk(posCg[0])))
             
 
         elif self.nDim == 3:
@@ -190,6 +224,8 @@ class Interpolator:
     def GetCm(self):
         return self.iSolver.Cm
 
+    def __RealToSections(self, var, varWk):
+        pass
 
     def __InterpToSections(self, var, varWk):
 
@@ -201,6 +237,35 @@ class Interpolator:
         var = np.reshape(var, [len(var), 1])
         completeWing = np.hstack((self.Wing, var))
 
+        wingUp = completeWing[completeWing[:,2]>=0, :]
+        wingLw = completeWing[completeWing[:,2]<=0, :]
+
+        mappedVar = [[np.zeros(0) for _ in range(3)] for _ in range(len(self.Sections))]
+
+        for iSec in range(len(mappedVar)):
+            for iReg in range(len(mappedVar[iSec])):
+                if iReg==0:
+                    interpFunc = interp1d(wingUp[:,0], wingUp[:,3])
+                if iReg==1:
+                    interpFunc = interp1d(wingLw[:,0], wingLw[:,3])
+                if iReg==2:
+                    interpFunc = interp1d(completeWake[:,0], completeWake[:,3])
+                
+                mappedVar[iSec][iReg] = interpFunc(self.Sections[iSec][iReg][:,0])
+
+        return mappedVar
+
+
+    def __RbfToSections(self, var, varWk):
+
+        varWk = np.reshape(varWk, [len(varWk), 1])
+        completeWake = np.hstack((self.Wake, varWk))
+        _, idx = np.unique(completeWake[:,:3], axis=0, return_index=True)
+        completeWake = completeWake[np.sort(idx)]
+
+        var = np.reshape(var, [len(var), 1])
+        completeWing = np.hstack((self.Wing, var))
+
         wingUp = completeWing[completeWing[:,2]>=0]
         wingLw = completeWing[completeWing[:,2]<=0]
 
@@ -209,11 +274,13 @@ class Interpolator:
         for iSec in range(len(mappedVar)):
             for iReg in range(len(mappedVar[iSec])):
                 if iReg==0:
-                    mappedVar[iSec][iReg] = RBFInterpolator(wingUp[:,:3], wingUp[:,3], neighbors=200, smoothing=0.1, kernel='linear', degree=0)(self.Sections[iSec][iReg])
+                    mappedVar[iSec][iReg] = RBFInterpolator(wingUp[:,:3], wingUp[:,3], neighbors=10, smoothing=0.1, epsilon=0.6, kernel='gaussian', degree=1)(self.Sections[iSec][iReg])
+                    #mappedVar[iSec][iReg] = RBFInterpolator(wingUp[:,:3], wingUp[:,3])(self.Sections[iSec][iReg])
                 if iReg==1:
-                    mappedVar[iSec][iReg] = RBFInterpolator(wingLw[:,:3], wingLw[:,3], neighbors=200, smoothing=0.1, kernel='linear', degree=0)(self.Sections[iSec][iReg])
+                    mappedVar[iSec][iReg] = RBFInterpolator(wingLw[:,:3], wingLw[:,3], neighbors=10, smoothing=0.1, epsilon=0.6, kernel='gaussian', degree=1)(self.Sections[iSec][iReg])
+                    #mappedVar[iSec][iReg] = RBFInterpolator(wingLw[:,:3], wingLw[:,3])(self.Sections[iSec][iReg])
                 if iReg==2:
-                    mappedVar[iSec][iReg] = RBFInterpolator(completeWake[:,:2], completeWake[:,3], smoothing=0.1, kernel='linear', degree=0)(self.Sections[iSec][iReg][:,:2])
+                    mappedVar[iSec][iReg] = RBFInterpolator(completeWake[:,:2], completeWake[:,3], neighbors=10, smoothing=0.1, kernel='linear', degree=0)(self.Sections[iSec][iReg][:,:2])
         return mappedVar
 
 
@@ -231,11 +298,19 @@ class Interpolator:
         """
 
         if reg==0:
-            vLw = np.insert(vectorLw, 0, np.flip(vectorUp[1:idx]))
-            vUp = np.delete(vectorUp, np.s_[:idx-1])
+            if idx!=0:
+                vLw = np.insert(vectorLw, 0, np.flip(vectorUp[1:idx]))
+                vUp = np.delete(vectorUp, np.s_[:idx-1])
+            else:
+                vLw = vectorLw
+                vUp = vectorUp
         elif reg==1:
-            vUp = np.insert(vectorUp, 0, np.flip(vectorLw[1:idx]))
-            vLw = np.delete(vectorLw, np.s_[:idx-1])
+            if idx!=0:
+                vUp = np.insert(vectorUp, 0, np.flip(vectorLw[1:idx]))
+                vLw = np.delete(vectorLw, np.s_[:idx-1])
+            else:
+                vUp = vectorUp
+                vLw = vectorLw
         else:
             raise RuntimeError('Incorrect region specified during remeshing.')
         
diff --git a/dart/pyVII/vInterpolator2.py b/dart/pyVII/vInterpolator2.py
new file mode 100644
index 0000000000000000000000000000000000000000..ab546c8ee6883727cc0e4281d37eda712b643e9f
--- /dev/null
+++ b/dart/pyVII/vInterpolator2.py
@@ -0,0 +1,396 @@
+import numpy as np
+from matplotlib import pyplot as plt
+from fwk.coloring import ccolors
+
+# 3D spline with Bézier curves
+from scipy.interpolate import bisplrep
+from scipy.interpolate import bisplev
+
+# 2D spline with Bézier curves
+from scipy.interpolate import splprep
+from scipy.interpolate import splev
+from scipy.integrate import cumtrapz
+
+from scipy.interpolate import interp1d
+
+from scipy.interpolate import RBFInterpolator
+
+import fwk
+
+
+
+class Interpolator:
+    def __init__(self, _dartSolver, bnd_wing, bnd_wake, _cfg):
+        self.iSolver = _dartSolver
+
+        self.bndWing = bnd_wing
+        self.bndWake = bnd_wake
+        self.invStruct, self.viscStruct = self.__GetWing(bnd_wing, bnd_wake, _cfg)
+        self.config = _cfg
+
+        self.xxPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
+        self.DeltaStarPrev = [[np.zeros(len(self.viscStruct[iSec][iReg][:,0])) for iReg in range(2)] for iSec in range(len(self.viscStruct))]
+        self.tms = fwk.Timers()
+
+    def GetInvBC(self, vSolver, couplIter):
+        self.tms['InvToVisc'].start()
+
+        # Extract parameters from the inviscid solver
+        # Particular care must be taken with space directions; In the inviscid, the airfoil is in the x,z plane.
+        # While setting up the viscous solver, y and z must be inverted. 
+
+        # Wing
+
+        nElm = len(self.bndWing.nodes)
+        U_wg, M_wg, Rho_wg = np.zeros((nElm, 3)), np.zeros(nElm), np.zeros(nElm)
+        for iNode, n in enumerate(self.bndWing.nodes):
+            U_wg[iNode,0] = self.iSolver.U[n.row][0]
+            U_wg[iNode,1] = self.iSolver.U[n.row][1]
+            U_wg[iNode,2] = self.iSolver.U[n.row][2]
+            M_wg[iNode] = self.iSolver.M[n.row]
+            Rho_wg[iNode] = self.iSolver.rho[n.row]
+        
+        # Wake
+
+        nElm = len(self.bndWake.nodes)
+        U_wk, M_wk, Rho_wk = np.zeros((nElm, 3)), np.zeros(nElm), np.zeros(nElm)
+        for iNode, n in enumerate(self.bndWake.nodes):
+            U_wk[iNode,0] = self.iSolver.U[n.row][0]
+            U_wk[iNode,1] = self.iSolver.U[n.row][2] # x and y inverted
+            U_wk[iNode,2] = self.iSolver.U[n.row][1]
+            M_wk[iNode] = self.iSolver.M[n.row]
+            Rho_wk[iNode] = self.iSolver.rho[n.row]
+
+        Ux = self.__RbfToSections(U_wg[:,0], U_wk[:,0])
+        Uy = self.__RbfToSections(U_wg[:,1], U_wk[:,1])
+        Uz = self.__RbfToSections(U_wg[:,2], U_wk[:,2])
+
+        Me = self.__RbfToSections(M_wg, M_wk)
+        Rho = self.__RbfToSections(Rho_wg, Rho_wk)
+
+        # Find stagnation point
+        for iSec in range(len(self.viscStruct)):
+
+            #if couplIter == 0:
+            self.stgPt = np.argmin(Ux[iSec][0]**2+Uy[iSec][0]**2)
+
+            xUp, xLw = self.__Remesh(self.viscStruct[iSec][0][:,0], self.stgPt)
+            yUp, yLw = self.__Remesh(self.viscStruct[iSec][0][:,1], self.stgPt)
+            zUp, zLw = self.__Remesh(self.viscStruct[iSec][0][:,2], self.stgPt)
+
+            vSolver.SetGlobMesh(iSec, 0, xUp, yUp, zUp)
+            vSolver.SetGlobMesh(iSec, 1, xLw, yLw, zLw)
+            vSolver.SetGlobMesh(iSec, 2, self.viscStruct[iSec][1][:,0],
+                                         self.viscStruct[iSec][1][:,1],
+                                         self.viscStruct[iSec][1][:,2])
+            
+            UxUp, UxLw = self.__Remesh(Ux[iSec][0], self.stgPt)
+            UyUp, UyLw = self.__Remesh(Uy[iSec][0], self.stgPt)
+            UzUp, UzLw = self.__Remesh(Uz[iSec][0], self.stgPt)
+
+            vSolver.SetVelocities(iSec, 0, UxUp, UyUp, UzUp)
+            vSolver.SetVelocities(iSec, 1, UxLw, UyLw, UzLw)
+            vSolver.SetVelocities(iSec, 2, Ux[iSec][1], Uz[iSec][1], Uy[iSec][1])
+
+            MeUp, MeLw = self.__Remesh(Me[iSec][0], self.stgPt)
+
+            vSolver.SetMe(iSec, 0, MeUp)
+            vSolver.SetMe(iSec, 1, MeLw)
+            vSolver.SetMe(iSec, 2, Me[iSec][1])
+
+            RhoUp, RhoLw = self.__Remesh(Rho[iSec][0], self.stgPt)
+
+            vSolver.SetRhoe(iSec, 0, RhoUp)
+            vSolver.SetRhoe(iSec, 1, RhoLw)
+            vSolver.SetRhoe(iSec, 2, Rho[iSec][1])
+
+            DeltaStarPrevUp, DeltaStarPrevLw = self.__Remesh(self.DeltaStarPrev[iSec][0], self.stgPt)
+            vSolver.SetDeltaStarExt(iSec, 0, DeltaStarPrevUp)
+            vSolver.SetDeltaStarExt(iSec, 1, DeltaStarPrevLw)
+            vSolver.SetDeltaStarExt(iSec, 2, self.DeltaStarPrev[iSec][1])
+
+            xxPrevUp, xxPrevLw = self.__Remesh(self.xxPrev[iSec][0], self.stgPt)
+            vSolver.SetxxExt(iSec, 0, xxPrevUp)
+            vSolver.SetxxExt(iSec, 1, xxPrevLw)
+            vSolver.SetxxExt(iSec, 2, self.xxPrev[iSec][1])
+            
+            
+
+    def SetBlowingVel(self, vSolver):
+        BlowingVel = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
+        for iSec in range(len(self.viscStruct)):
+            for iRegion in range(3):
+                BlowingVel[iSec][iRegion] = np.asarray(vSolver.ExtractBlowingVel(iSec, iRegion))
+
+
+        if self.config['nDim'] == 2:
+            BlwWing = np.concatenate((BlowingVel[0][0][::-1], BlowingVel[0][1]))
+            BlwWake = BlowingVel[0][2]
+            self.DeltaStarPrev[iSec][0] = np.concatenate((np.asarray(vSolver.ExtractDeltaStar(iSec, 0))[::-1], np.asarray(vSolver.ExtractDeltaStar(iSec, 1))[1:]), axis=None)
+            self.DeltaStarPrev[iSec][1] = np.asarray(vSolver.ExtractDeltaStar(iSec, 2))
+            self.xxPrev[iSec][0] = np.concatenate((np.asarray(vSolver.Extractxx(iSec, 0))[::-1], np.asarray(vSolver.Extractxx(iSec, 1))[1:]), axis=None)
+            self.xxPrev[iSec][1] = np.asarray(vSolver.Extractxx(iSec, 2))
+
+            # Compute cell centers
+            viscCgWing = np.zeros((len(self.viscStruct[0][0][:,0]) - 1, 3))
+            for iElm in range(len(self.viscStruct[0][0][:,0]) - 1):
+                viscCgWing[iElm, 0] = self.viscStruct[0][0][iElm, 0] + self.viscStruct[0][0][iElm + 1, 0]
+                viscCgWing[iElm, 1] = self.viscStruct[0][0][iElm, 1] + self.viscStruct[0][0][iElm + 1, 1]
+                viscCgWing[iElm, 2] = self.viscStruct[0][0][iElm, 2] + self.viscStruct[0][0][iElm + 1, 2]
+            viscCgWing/=2
+
+            viscCgWake = np.zeros((len(self.viscStruct[0][1][:,0]) - 1, 3))
+            for iElm in range(len(self.viscStruct[0][1][:,0]) - 1):
+                viscCgWake[iElm, 0] = self.viscStruct[0][1][iElm, 0] + self.viscStruct[0][1][iElm + 1, 0]
+                viscCgWake[iElm, 1] = self.viscStruct[0][1][iElm, 1] + self.viscStruct[0][1][iElm + 1, 1]
+                viscCgWake[iElm, 2] = self.viscStruct[0][1][iElm, 2] + self.viscStruct[0][1][iElm + 1, 2]
+            viscCgWake/=2
+
+            invCgWing = np.zeros((len(self.bndWing.tag.elems), 3))
+            for iElm, e in enumerate(self.bndWing.tag.elems):
+                for n in e.nodes:
+                    invCgWing[iElm, 0] +=n.pos[0]
+                    invCgWing[iElm, 1] +=n.pos[1]
+                    invCgWing[iElm, 2] +=n.pos[2]
+                invCgWing[iElm, :] /= e.nodes.size()
+            invCgWake = np.zeros((len(self.bndWake.tag.elems), 3))
+            for iElm, e in enumerate(self.bndWake.tag.elems):
+                for n in e.nodes:
+                    invCgWake[iElm, 0] +=n.pos[0]
+                    invCgWake[iElm, 1] +=n.pos[1]
+                    invCgWake[iElm, 2] +=n.pos[2]
+                invCgWake[iElm, :] /= e.nodes.size()
+            
+            """plt.plot(self.invStruct[0][:,0], self.invStruct[0][:,1], 'x')
+            plt.plot(invCgWing[:, 0], invCgWing[:,1], 'o')
+            plt.show()"""
+
+        self.tms['Interpolation'].start()
+        mappedBlwWing = RBFInterpolator(viscCgWing, BlwWing, neighbors=5, kernel='linear', smoothing=1e-6)(invCgWing)
+        mappedBlwWake = RBFInterpolator(viscCgWake, BlwWake, neighbors=5, kernel='linear', smoothing=1e-6)(invCgWake)
+        self.tms['Interpolation'].stop()
+        for iElm in range(len(self.bndWing.tag.elems)):
+            self.bndWing.setU(iElm, mappedBlwWing[iElm])
+        for iElm in range(len(self.bndWake.tag.elems)):
+            self.bndWake.setU(iElm, mappedBlwWake[iElm])
+
+
+    def RunSolver(self):
+        self.iSolver.run()
+
+    def ResetSolver(self):
+        self.iSolver.reset()
+
+    def GetCp(self,bnd):
+        import dart.utils as floU
+
+        Cp = floU.extract(bnd, self.iSolver.Cp)
+        return Cp
+
+    def GetAlpha(self):
+        return self.iSolver.pbl.alpha
+
+    def GetMinf(self):
+        return self.iSolver.pbl.M_inf
+    
+    def GetCl(self):
+        return self.iSolver.Cl
+    
+    def GetCm(self):
+        return self.iSolver.Cm
+
+    def __Remesh(self, vec, idx_stag):
+        xUp = vec[:idx_stag+1]
+        xUp = xUp[::-1]
+
+        xLw = vec[idx_stag:]
+        return xUp, xLw
+
+    def __RbfToSections(self, var, varWk):
+        varWk = np.reshape(varWk, [len(varWk), 1])
+        completeWake = np.hstack((self.invStruct[1], varWk))
+        _, idx = np.unique(completeWake[:,:3], axis=0, return_index=True)
+        completeWake = completeWake[np.sort(idx)]
+
+        var = np.reshape(var, [len(var), 1])
+        completeWing = np.hstack((self.invStruct[0], var))
+
+        mappedVar = [[np.zeros(0) for _ in range(3)] for _ in range(len(self.viscStruct))]
+
+        for iSec in range(len(mappedVar)):
+            for iReg in range(len(mappedVar[iSec])):
+                if iReg==0: # Wing
+                    mappedVar[iSec][iReg] = RBFInterpolator(completeWing[:,:3], completeWing[:,3], neighbors=5, kernel='linear', smoothing=1e-2)(self.viscStruct[iSec][0])
+                if iReg==1: # Wake
+                    mappedVar[iSec][iReg] = RBFInterpolator(completeWake[:,:3], completeWake[:,3], neighbors=5, kernel='linear', smoothing=1e-2)(self.viscStruct[iSec][1])
+        
+        """plt.plot(completeWing[:,0], completeWing[:,3], 'x')
+        plt.plot(completeWake[:,0], completeWake[:,3], 'x')
+        plt.plot(self.viscStruct[0][0][:,0], mappedVar[0][0])
+        plt.plot(self.viscStruct[0][1][:,0], mappedVar[0][1])
+        plt.show()"""
+        return mappedVar
+
+
+    def __GetWing(self, bnd_wing, bnd_wake, config):
+
+        if config['nDim'] == 2 and config['nSections'] != 1:
+            raise RuntimeError('Cannot create more than 1 section on 2D geometry. Specified:', config['nSections'])
+
+        viscStruct = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0], 3)) for i in range(2)] for _ in range(config['nSections'])]
+
+        invStruct = [np.zeros((len(bnd_wake.nodes) if iReg==1 else len(bnd_wing.nodes), 3)) for iReg in range(2)]
+
+        for iNode, n in enumerate(bnd_wing.nodes):
+            invStruct[0][iNode, 0] = n.pos[0]
+            invStruct[0][iNode, 1] = n.pos[1]
+            invStruct[0][iNode, 2] = n.pos[2]
+        for iNode, n in enumerate(bnd_wake.nodes):
+            invStruct[1][iNode, 0] = n.pos[0]
+            invStruct[1][iNode, 1] = n.pos[2]
+            invStruct[1][iNode, 2] = n.pos[1]
+
+        if config['nDim'] == 2:
+            spanDistri = np.zeros(1)
+            spanDistri[0] = invStruct[0][0,2] # Span direction is the z direction 
+
+        elif config['nDim'] == 3:
+            invStruct[0][:,[1,2]] = invStruct[0][:,[2,1]] # Inverse y and z so that spanwise direction is in z
+            invStruct[1][:,[1,2]] = invStruct[1][:,[2,1]] # Inverse y and z so that spanwise direction is in z
+            spanDistri = np.linspace(min(invStruct[0][:,2]), max(invStruct[0][:,2]), config['nSections'])
+            spanDistri[0] +=0.1
+            spanDistri[-1] -=0.1
+
+        for iSec, z in enumerate(spanDistri):
+
+            # Create x distribution at the corresponding span position
+            minX = min(invStruct[0][abs(invStruct[0][:,2]-z)<=config['eps'], 0])
+            chordLength = max(invStruct[0][abs(invStruct[0][:,2]-z)<=config['eps'], 0]) - minX
+            xDistri = minX + 1/2*(1 + np.cos(np.linspace(0, np.pi, config['nPoints'][0])))*chordLength
+            nElmSide = len(xDistri)
+            viscStruct[iSec][0][:,0] = np.concatenate((xDistri, xDistri[::-1]), axis=None)
+            viscStruct[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
+
+            wakeLength = np.max(invStruct[1][:,0]) - np.max(viscStruct[iSec][0][:,0])
+            viscStruct[iSec][1][:,0] = np.max(viscStruct[iSec][0][:,0]) + (np.max(viscStruct[iSec][0][:,0]) + np.cos(np.linspace(np.pi, np.pi/2, config['nPoints'][1]))) * wakeLength
+            viscStruct[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
+
+            if config['nDim'] == 2:
+                upper, lower = self.__SeparateArf(invStruct[0][:,0], invStruct[0][:,1])
+                viscStruct[iSec][0][:nElmSide,1] = interp1d(upper[:,0], upper[:,1])(viscStruct[iSec][0][:nElmSide,0])
+                viscStruct[iSec][0][nElmSide:,1] = interp1d(lower[:,0], lower[:,1])(viscStruct[iSec][0][nElmSide:,0])
+                viscStruct[iSec][0] = np.delete(viscStruct[iSec][0], nElmSide, axis=0)
+
+                viscStruct[iSec][1][:,1] = interp1d(invStruct[1][:,0], invStruct[1][:,1])(viscStruct[iSec][1][:,0])
+
+                """plt.plot(viscStruct[iSec][0][:,0], viscStruct[iSec][0][:,1], '-')
+                plt.plot(viscStruct[iSec][0][:nElmSide,0], viscStruct[iSec][0][:nElmSide,1], 'x')
+                plt.plot(viscStruct[iSec][0][nElmSide:,0], viscStruct[iSec][0][nElmSide:,1], 'o')
+                plt.plot(viscStruct[iSec][1][:,0], np.zeros(len(viscStruct[iSec][1][:,0])), 'x-')
+                plt.xlim([-0.5,2])
+                plt.ylim([-0.5, 0.5])
+                plt.show()"""
+
+            elif config['nDim'] == 3:
+                nPointsMin = (config['k'][0]+1) * (config['k'][1]+1) # See Scipy bisplrep doc
+
+                wingSpan = max(invStruct[0][:,2]) - min(invStruct[0][:,2])
+                portion = invStruct[0][abs(invStruct[0][:,2] - z) < 0.01 * wingSpan,:]
+
+                portionUpper, portionLower = self.__SeparateArf(portion[:,0], portion[:,1], portion[:,2])
+                if len(portionUpper)<=nPointsMin or len(portionLower)<=nPointsMin:
+                    print(ccolors.ANSI_YELLOW,'dart::Interpolator Warning: Found {:.0f} points within 1% of spanwise position z={:.2f}. Trying 10%. '.format(
+                        len(portion), z),ccolors.ANSI_RESET)
+                    portion = invStruct[0][abs(invStruct[0][:,2] - z) < 0.1 * wingSpan,:]
+                    if len(portion)<=nPointsMin:
+                        print(ccolors.ANSI_RED, 'Failed. Continuing\n', ccolors.ANSI_RESET)
+                    else:
+                        print(ccolors.ANSI_GREEN, 'Sucess.\n', ccolors.ANSI_RESET)
+
+                # Upper
+                sWing = (len(portionUpper[:,0])-np.sqrt(2*len(portionUpper[:,0])) + len(portionUpper[:,0])+np.sqrt(2*len(portionUpper[:,0])))/2
+                tckWingUp = bisplrep(portionUpper[:,0], portionUpper[:,2], portionUpper[:,1], kx=config['k'][0], ky=config['k'][1], s=sWing)
+                for iPoint in range(nElmSide):
+                    viscStruct[iSec][0][iPoint,1] = bisplev(viscStruct[iSec][0][iPoint,0], viscStruct[iSec][0][iPoint,2], tckWingUp)
+                
+                # Lower
+                sWing = (len(portionLower[:,0])-np.sqrt(2*len(portionLower[:,0])) + len(portionLower[:,0])+np.sqrt(2*len(portionLower[:,0])))/2
+                tckWingUp = bisplrep(portionLower[:,0], portionLower[:,2], portionLower[:,1], kx=config['k'][0], ky=config['k'][1], s=sWing)
+                for iPoint in range(nElmSide, len(viscStruct[iSec][0][:,0])):
+                    viscStruct[iSec][0][iPoint,1] = bisplev(viscStruct[iSec][0][iPoint,0], viscStruct[iSec][0][iPoint,2], tckWingUp)
+
+                # Wake
+                sWake = (len(invStruct[1][:,0])-np.sqrt(2*len(invStruct[1][:,0])) + len(invStruct[1][:,0])+np.sqrt(2*len(invStruct[1][:,0])))/2
+                tckWingUp = bisplrep(invStruct[1][:,0], invStruct[1][:,2], invStruct[1][:,1], kx=config['k'][0], ky=config['k'][1], s=sWake)
+                for iPoint in range(len(viscStruct[iSec][1][:,0])):
+                    viscStruct[iSec][1][iPoint,1] = bisplev(viscStruct[iSec][1][iPoint,0], viscStruct[iSec][1][iPoint,2], tckWingUp)
+
+
+        if config['nDim'] ==3:
+            fig = plt.figure()
+            ax = fig.add_subplot(projection='3d')
+            ax.scatter(invStruct[0][:,0], invStruct[0][:,1], invStruct[0][:,2], s=0.2)
+            for iSec in range(len(viscStruct)):
+                ax.scatter(viscStruct[iSec][0][:,0], viscStruct[iSec][0][:,1], viscStruct[iSec][0][:,2], s=0.9, color='red')
+            ax.set_ylim([-0.5, 0.5])
+            ax.set_xlabel('x')
+            ax.set_ylabel('y')
+            ax.set_zlabel('z')
+            plt.show()
+        return invStruct, viscStruct
+
+    def __SeparateArf(self, xInp, yInp, zInp=None):
+        if zInp is None:
+            zInp = np.zeros(len(xInp))
+        if len(xInp) != len(yInp) or len(xInp) != len(zInp):
+            raise RuntimeError ('dart::Interpolator::__SeparateArf Missmatch between intput vector sizes.')
+
+        Q=np.zeros((len(xInp), 3))
+        Q[:,0] = xInp.copy()
+        Q[:,1] = yInp.copy()
+        Q[:,2] = zInp.copy()
+
+        Q = Q[Q[:, 0].argsort()]
+
+        le = Q[0,:]
+        te = Q[-1,:]
+        ymean = (le[1] + te[1]) / 2
+
+        upper, lower, save = le, le, le
+
+        flag1 = 0
+        iPoint = 0
+        for i in range(len(Q[:,0])):
+            if Q[iPoint,1] == 0.0:
+                if flag1 == 1:
+                    Q=np.delete(Q, iPoint, axis = 0)
+                    iPoint-=1
+                else:
+                    flag1 = 1
+            iPoint +=1
+
+        for iPoint in range(1, len(Q[:,0])):
+            if Q[iPoint,1] > ymean:
+                upper = np.vstack((upper, Q[iPoint,:]))
+            elif Q[iPoint, 1] < ymean:
+                lower = np.vstack((lower, Q[iPoint,:]))
+            else:
+                save = np.vstack((save, Q[iPoint, :]))
+        
+        _, idx = np.unique(save, return_index=True, axis=0)
+        save=save[np.sort(idx)]
+
+        #save = np.delete(save, 0, axis=0)
+
+        
+        try:
+            if len(save[:,0])>0:
+                upper = np.vstack((upper, save))
+                lower = np.vstack((lower, save))
+        except:
+            pass
+
+        upper = np.vstack((upper, te))
+        lower = np.vstack((lower, te))
+
+        return upper, lower
\ No newline at end of file
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index c1ed9ac5717d3cb311712ce3fb061a8eb8917baf..5cf782b8d8d2bae78bcc6409e7c5844903cbb631 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -36,6 +36,7 @@ void BLRegion::SetMesh(std::vector<double> x,
 
     if (mesh->GetDiffnElm()!=0)
     {
+        std::cout << "Resizing vectors" << std::endl;
         U.resize(nVar * nMarkers, 0);
         Ret.resize(nMarkers, 0);
         Cd.resize(nMarkers);
diff --git a/dart/src/wClosures.cpp b/dart/src/wClosures.cpp
index 4a91354875a4eee6dca9804a71d83bf94e5c7c03..6aa84a0c1b343387767ce1e278e13058c0454db9 100644
--- a/dart/src/wClosures.cpp
+++ b/dart/src/wClosures.cpp
@@ -17,120 +17,180 @@ Closures::~Closures()
 std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret, double Me, std::string name)
 {
 
-    double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me); // Kinematic shape factor
-    double Hstar2 = (0.064 / (Hk - 0.8) + 0.251) * Me * Me;   // Density shape parameter
-    double delta = (3.15 + H + (1.72 / (Hk - 1))) * theta;
+    H = std::max(H, 1.00005);
+
+    /* Kinematic shape factor H* */
+
+    double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
+    if (name == "wake")
+    {
+        Hk = std::max(Hk, 1.00005);
+    }
+    else
+    {
+        Hk = std::max(Hk, 1.05000);
+    }
+
+    /* Density shape parameter */
+    
+    double Hstar2 = (0.064/(Hk-0.8)+0.251)*Me*Me;
+
+    /* Boundary layer thickness */
+
+    double delta = std::min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta);
+
+    double Hstar = 0.0;
     double Hks = Hk - 4.35;
-    double Hstar;
     if (Hk <= 4.35)
     {
-        double dens = Hk + 1;
-        Hstar = 0.0111 * Hks * Hks / dens - 0.0278 * Hks * Hks * Hks / dens + 1.528 - 0.0002 * (Hks * Hk) * (Hks * Hk);
-        Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me); // Whitfield's minor additional compressibility correction
+        Hstar = 0.0111*Hks*Hks/(Hk+1) - 0.0278*Hks*Hks*Hks/(Hk+1) + 1.528 -0.0002*(Hks*Hk)*(Hks*Hk);
     }
-    else if (Hk > 4.35)
+    else
     {
-        Hstar = 1.528 + 0.015 * Hks * Hks / Hk;
-        Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me); // Whitfield's minor additional compressibility correction
+        Hstar = 1.528 + 0.015*Hks*Hks/Hk;
     }
+    // Whitfield's minor additional compressibility correction
+    Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me);
+
+    /* Friction coefficient */
 
-    double tmp;
-    double Cf;
-    double Cd;
+    double tmp = 0.0;
+    double Cf = 0.0;
     if (Hk < 5.5)
     {
-        tmp = (5.5 - Hk) * (5.5 - Hk) * (5.5 - Hk) / (Hk + 1);
-        Cf = (-0.07 + 0.0727 * tmp) / Ret;
+        tmp = (5.5-Hk)*(5.5-Hk)*(5.5-Hk) / (Hk+1.0);
+        Cf = (-0.07 + 0.0727*tmp)/Ret;
     }
     else if (Hk >= 5.5)
     {
-        tmp = 1 - 1 / (Hk - 4.5);
-        Cf = (-0.07 + 0.015 * tmp * tmp) / Ret;
+        tmp = 1.0 - 1.0/(Hk-4.5);
+        Cf = (-0.07 + 0.015*tmp*tmp) /Ret;
     }
+
+    /* Dissipation coefficient */
+
+    double Cd = 0.0;
     if (Hk < 4)
     {
-        Cd = (0.00205 * pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2 * Ret));
+        Cd = ( 0.00205*std::pow(4.0-Hk, 5.5) + 0.207 ) * (Hstar/(2*Ret));
     }
-    else if (Hk >= 4)
+    else
     {
-        double HkCd = (Hk - 4) * (Hk - 4);
-        double denCd = 1 + 0.02 * HkCd;
-        Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2 * Ret));
+        double HkCd = (Hk-4)*(Hk-4);
+        double denCd = 1+0.02*HkCd;
+        Cd = ( -0.0016*HkCd/denCd + 0.207) * (Hstar/(2*Ret));
     }
+
+    /* Wake relations */
+    
     if (name == "wake")
     {
-        Cd = 2 * (1.10 * (1.0 - 1.0 / Hk) * (1.0 - 1.0 / Hk) / Hk) * (Hstar / (2 * Ret));
+        Cd = 2*(1.10*(1.0-1.0/Hk)*(1.0-1.0/Hk)/Hk)* (Hstar/(2*Ret));
         Cf = 0;
     }
 
-    double Cteq = 0;
-    double Us = 0;
+    double Us = 0.0;
+    double Cteq = 0.0;
 
-    std::vector<double> ClosureVars;
+    std::vector<double> ClosureVars(8, 0.0);
 
-    ClosureVars.push_back(Hstar);
-    ClosureVars.push_back(Hstar2);
-    ClosureVars.push_back(Hk);
-    ClosureVars.push_back(Cf);
-    ClosureVars.push_back(Cd);
-    ClosureVars.push_back(delta);
-    ClosureVars.push_back(Cteq);
-    ClosureVars.push_back(Us);
+    ClosureVars[0] = Hstar;
+    ClosureVars[1] = Hstar2;
+    ClosureVars[2] = Hk;
+    ClosureVars[3] = Cf;
+    ClosureVars[4] = Cd;
+    ClosureVars[5] = delta;
+    ClosureVars[6] = Cteq;
+    ClosureVars[7] = Us;
     return ClosureVars;
 }
 
 std::vector<double> Closures::TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name)
 {
 
-    double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
-    Hk = std::max(Hk, 1.0005);
-    double Hstar2 = ((0.064 / (Hk - 0.8)) + 0.251) * Me * Me;
+    H = std::max(H, 1.00005);
+    Ct = std::min(Ct, 0.3);
+    //Ct = std::max(std::min(Ct, 0.30), 0.0000001);
 
-    double Fc = sqrt(1 + 0.2 * Me * Me);
+    double Hk = (H - 0.29*Me*Me)/(1+0.113*Me*Me);
+    if (name == "wake")
+    {
+        Hk = std::max(Hk, 1.00005);
+    }
+    else
+    {
+        Hk = std::max(Hk, 1.05000);
+    }
 
-    double H0;
+    double Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me*Me;
+    double gamma = 1.4;
+    double Fc = std::sqrt(1+0.5*gamma*Me*Me);
+
+    double H0 = 0.0;
     if (Ret <= 400)
     {
-        H0 = 4;
+        H0 = 4.0;
     }
-    else
+    else if (Ret > 400)
     {
-        H0 = 3 + (400 / Ret);
+        H0 = 3 + (400/Ret);
+    }
+    if (Ret <= 200)
+    {
+        Ret = 200;
     }
 
-    double Hstar;
+    double Hstar = 0.0;
     if (Hk <= H0)
     {
-        Hstar = ((0.5 - 4 / Ret) * ((H0 - Hk) / (H0 - 1)) * ((H0 - Hk) / (H0 - 1))) * (1.5 / (Hk + 0.5)) + 1.5 + 4 / Ret;
+        Hstar = ((0.5-4/Ret)*((H0-Hk)/(H0-1))*((H0-Hk)/(H0-1)))*(1.5/(Hk+0.5))+1.5+4/Ret;
     }
-    else
+    else if (Hk > H0)
     {
-        Hstar = (Hk - H0) * (Hk - H0) * (0.007 * log(Ret) / ((Hk - H0 + 4 / log(Ret)) * (Hk - H0 + 4 / log(Ret))) + 0.015 / Hk) + 1.5 + 4 / Ret;
+        Hstar = (Hk-H0)*(Hk-H0)*(0.007*std::log(Ret)/((Hk-H0+4/std::log(Ret))*(Hk-H0+4/std::log(Ret))) + 0.015/Hk) + 1.5 + 4/Ret;
     }
-    Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me); /* Whitfield's minor additional compressibility correction */
+    /* Whitfield's minor additional compressibility correction */ 
+    Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me);
+
+
+    double logRt = std::log(Ret/Fc);
+    logRt = std::max(logRt, 3.0);
+    double arg = std::max(-1.33*Hk, -20.0);
+
+    /* Equivalent normalized wall slip velocity */
+    double Us = (Hstar/2)*(1-4*(Hk-1)/(3*H));
 
-    double Us = (Hstar / 2) * (1 - 4 * (Hk - 1) / (3 * H)); /* Equivalent normalized wall slip velocity */
-    double delta = std::min((3.15 + H + (1.72 / (Hk - 1))) * theta, 12 * theta);
-    double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
+    /* Boundary layer thickness */
 
-    double Cf;
-    double Cd;
-    double Hkc;
-    double Cdw;
-    double Cdd;
-    double Cdl;
-    double Cteq;
-    int n;
+    double delta = std::min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta);
+
+    double Ctcon = 0.5/(6.7*6.7*0.75);
+    double Hkc = 0.0;
+    double Cdw = 0.0;
+    double Cdd = 0.0;
+    double Cdl = 0.0;
+
+    double Hmin = 0.0;
+    double Fl = 0.0;
+    double Dfac = 0.0;
+
+    double Cf = 0.0;
+    double Cd = 0.0;
+
+    double n = 1.0;
+    
     if (name == "wake")
     {
-        Us = std::min(Us, 0.99995);
-        n = 2;  // two wake halves
-        Cf = 0; // no friction inside the wake
+        if (Us > 0.99995)
+        {
+            Us = 0.99995;
+        }
+        n = 2.0; // two wake halves
+        Cf = 0.0; // no friction inside the wake
         Hkc = Hk - 1;
-        Cdw = 0;                                                                              // wall contribution
-        Cdd = (0.995 - Us) * Ct * Ct;                                                         // turbulent outer layer contribution
-        Cdl = 0.15 * (0.995 - Us) * (0.995 - Us) / Ret;                                       // laminar stress contribution to outer layer
-        Cteq = sqrt(4 * Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - Us) * (Hk * Hk) * H)); // Here it is Cteq^0.5
+        Cdw = 0.0;  // wall contribution
+        Cdd = (0.995-Us)*Ct*Ct; // turbulent outer layer contribution
+        Cdl = 0.15*(0.995-Us)*(0.995-Us)/Ret; // laminar stress contribution to outer layer
     }
     else
     {
@@ -138,30 +198,36 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
         {
             Us = 0.98;
         }
-        n = 1;
-        Cf = 1 / Fc * (0.3 * exp(-1.33 * Hk) * pow((log10(Ret / Fc)), -1.74 - 0.31 * Hk) + 0.00011 * (tanh(4 - Hk / 0.875) - 1));
-        Hkc = std::max(Hk - 1 - 18 / Ret, 0.01);
-
-        // Dissipation coefficient
-        double Hmin = 1 + 2.1 / log(Ret);
-        double Fl = (Hk - 1) / (Hmin - 1);
-        double Dfac = 0.5 + 0.5 * tanh(Fl);
-        Cdw = 0.5 * (Cf * Us) * Dfac;
-        Cdd = (0.995 - Us) * Ct * Ct;
-        Cdl = 0.15 * (0.995 - Us) * (0.995 - Us) / Ret;
-        Cteq = sqrt(Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - Us) * (Hk * Hk) * H)); // Here it is Cteq^0.5
-    }
-    Cd = n * (Cdw + Cdd + Cdl);
-
-    std::vector<double> ClosureVars;
-
-    ClosureVars.push_back(Hstar);
-    ClosureVars.push_back(Hstar2);
-    ClosureVars.push_back(Hk);
-    ClosureVars.push_back(Cf);
-    ClosureVars.push_back(Cd);
-    ClosureVars.push_back(delta);
-    ClosureVars.push_back(Cteq);
-    ClosureVars.push_back(Us);
+        n = 1.0;
+        Cf = (1/Fc)*(0.3*std::exp(arg)* std::pow(logRt/2.3026, (-1.74-0.31*Hk)) + 0.00011*(std::tanh(4-(Hk/0.875))-1));
+        Hkc = std::max(Hk-1-18/Ret, 0.01);
+        /* dissipation coefficient */
+        Hmin = 1 + 2.1/std::log(Ret);
+        Fl = (Hk-1)/(Hmin-1);
+        Dfac = 0.5+0.5*std::tanh(Fl);
+        Cdw = 0.5*(Cf*Us) * Dfac;
+        Cdd = (0.995-Us)*Ct*Ct;
+        Cdl = 0.15*(0.995-Us)*(0.995-Us)/Ret;
+    }
+    if (n*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)<0)
+    {
+        std::cout << "Negative sqrt encountered " << std::endl;
+        throw;
+    }
+    double Cteq = std::sqrt(n*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H));
+    Cd = n*(Cdw+ Cdd + Cdl);
+    // Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk))**2)
+
+    std::vector<double> ClosureVars(8, 0.0);
+
+    ClosureVars[0] = Hstar;
+    ClosureVars[1] = Hstar2;
+    ClosureVars[2] = Hk;
+    ClosureVars[3] = Cf;
+    ClosureVars[4] = Cd;
+    ClosureVars[5] = delta;
+    ClosureVars[6] = Cteq;
+    ClosureVars[7] = Us;
     return ClosureVars;
-}
+
+}
\ No newline at end of file
diff --git a/dart/src/wInvBnd.h b/dart/src/wInvBnd.h
index 77a88cbbe0ca5bebfe51c3abf7b309be82492399..0fa1fa85ef2ef7c3424fd06a456b6bd6e5ff0be2 100644
--- a/dart/src/wInvBnd.h
+++ b/dart/src/wInvBnd.h
@@ -29,10 +29,10 @@ class DART_API InvBnd
     double GetRhoe(size_t iPoint) const {return Rhoe[iPoint];};
     double GetDeltaStar(size_t iPoint) const {return DeltaStar[iPoint];};
 
-    void SetMe(std::vector<double> _Me) {Me = _Me;};
-    void SetVt(std::vector<double> _Vt) {Vt = _Vt;};
-    void SetRhoe(std::vector<double> _Rhoe) {Rhoe = _Rhoe;};
-    void SetDeltaStar(std::vector<double> _DeltaStar) {DeltaStar = _DeltaStar;};
+    void SetMe(std::vector<double> _Me) {Me.resize(_Me.size(), 0.0); Me = _Me;};
+    void SetVt(std::vector<double> _Vt) {Vt.resize(_Vt.size(), 0.0); Vt = _Vt;};
+    void SetRhoe(std::vector<double> _Rhoe) {Rhoe.resize(_Rhoe.size(), 0.0); Rhoe = _Rhoe;};
+    void SetDeltaStar(std::vector<double> _DeltaStar) {DeltaStar.resize(_DeltaStar.size(), 0.0); DeltaStar = _DeltaStar;};
 
 };
 } // namespace dart
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index 49a5d4f0044b695d3bea6aaba21926526db5e2f4..d5f16791fa70f2248607b445c133f9fb38456145 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -53,15 +53,15 @@ void TimeSolver::InitialCondition(size_t iPoint, BLRegion *bl)
         bl->Delta[iPoint] = bl->Delta[iPoint - 1];
         bl->DeltaStar[iPoint] = bl->DeltaStar[iPoint - 1];
     }
-    if (bl->Regime[iPoint] == 1 && bl->U[iPoint * nVar + 4] == 0)
+    if (bl->Regime[iPoint] == 1 && bl->U[iPoint * nVar + 4] <= 0)
     {
         bl->U[iPoint * nVar + 4] = 0.03;
     } // End if
 
-    if (bl->U[iPoint*nVar+3]<0)
+    /*if (bl->U[iPoint*nVar+3]<0)
     {
         std::cout << "dart::TimeSolver negative velocity encountered at point " << iPoint << std::endl;
-    }
+    }*/
 } // End InitialCondition
 
 int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
@@ -71,10 +71,10 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
 
     /* Save initial condition */
 
-    std::vector<double> Sln0;
+    std::vector<double> Sln0(nVar, 0.0);
     for (size_t i = 0; i < nVar; ++i)
     {
-        Sln0.push_back(bl->U[iPoint * nVar + i]);
+        Sln0[i] = bl->U[iPoint * nVar + i];
     } // End for
 
     /* Initialize time integration variables */
@@ -115,10 +115,23 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
         {
             bl->U[iPoint * nVar + k] += slnIncr(k);
         }
-
+        bl->U[iPoint * nVar + 0] = std::max(bl->U[iPoint * nVar + 0], 1e-6);
+        bl->U[iPoint * nVar + 1] = std::max(bl->U[iPoint * nVar + 1], 1.0005);
         SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
         normSysRes = (SysRes + slnIncr / timeStep).norm();
 
+        if (std::isnan(normSysRes))
+        {
+            ResetSolution(iPoint, bl, Sln0, true);
+            SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
+            CFL = 1.0;
+            timeStep = SetTimeStep(CFL, dx, bl->invBnd->GetVt(iPoint));
+            IMatrix = Matrix<double, 5, 5>::Identity() / timeStep;
+            normSysRes = (SysRes + slnIncr / timeStep).norm();
+            itFrozenJac = 1;
+
+        }
+
         if (normSysRes <= tol * normSysRes0)
         {
             return 0; /* Successfull exit */
@@ -132,8 +145,10 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
 
         innerIters++;
     } // End while (innerIters < maxIt)
-
-    ResetSolution(iPoint, bl, Sln0);
+    if (std::isnan(normSysRes) || normSysRes/normSysRes0 > 1e-3)
+    {
+        ResetSolution(iPoint, bl, Sln0);
+    }
     return innerIters;
 } // End Integration
 
diff --git a/dart/src/wTimeSolver.h b/dart/src/wTimeSolver.h
index ff0d6900e45609594c4d5704688a2cc04eba6c31..e45d68906cfc54f5f77a7c1a153e65fbd487e5a6 100644
--- a/dart/src/wTimeSolver.h
+++ b/dart/src/wTimeSolver.h
@@ -25,7 +25,7 @@ private:
     ViscFluxes *SysMatrix;
 
 public:
-    TimeSolver(double _CFL0, double _Minf, double Re, unsigned long _maxIt = 50, double _tol = 1e-8, unsigned int _itFrozenJac = 5);
+    TimeSolver(double _CFL0, double _Minf, double Re, unsigned long _maxIt = 100, double _tol = 1e-8, unsigned int _itFrozenJac = 5);
     ~TimeSolver();
     void InitialCondition(size_t iPoint, BLRegion *bl);
     int Integration(size_t iPoint, BLRegion *bl);
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index 709040486c8ce4ef75996f19bb684b3e3f9af11c..60854e3de990330597fcd7e597e2282fc16d1293 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -26,10 +26,10 @@ ViscFluxes::~ViscFluxes()
 Vector<double, 5> ViscFluxes::ComputeFluxes(size_t iPoint, BLRegion *bl)
 {
     unsigned int nVar = bl->GetnVar();
-    std::vector<double> u;
+    std::vector<double> u(nVar, 0.0);
     for (size_t k = 0; k < nVar; ++k)
     {
-        u.push_back(bl->U[iPoint * nVar + k]);
+        u[k] = bl->U[iPoint * nVar + k];
     } // End for
     return BLlaws(iPoint, bl, u);
 }
@@ -39,13 +39,13 @@ Matrix<double, 5, 5> ViscFluxes::ComputeJacobian(size_t iPoint, BLRegion *bl, Ve
 
     unsigned int nVar = bl->GetnVar();
     Matrix<double, 5, 5> JacMatrix;
-    std::vector<double> uUp;
+    std::vector<double> uUp(nVar, 0.0);
     for (size_t k = 0; k < nVar; ++k)
     {
-        uUp.push_back(bl->U[iPoint * nVar + k]);
+        uUp[k] = bl->U[iPoint * nVar + k];
     } // End for
 
-    double varSave;
+    double varSave = 0.0;
     for (size_t iVar = 0; iVar < nVar; ++iVar)
     {
         varSave = uUp[iVar];
@@ -75,7 +75,7 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
         dissipRatio = 1;
     } // End else
 
-    bl->Ret[iPoint] = u[3] * u[0] * Re;  /* Reynolds number based on theta */
+    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 */
@@ -129,8 +129,8 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
     double cExt = 4 / (M_PI * dxExt);
 
     /* Amplification routine */
-    double dN_dx;
-    double Ax;
+    double dN_dx = 0;
+    double Ax = 0;
 
     if (bl->xtrF != -1 && bl->Regime[iPoint] == 0)
     {
@@ -138,11 +138,6 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
         dN_dx = (u[2] - bl->U[(iPoint - 1) * nVar + 2]) / dx;
 
     } // End if
-    else
-    {
-        Ax = 0;
-        dN_dx = 0;
-    } // End else
 
     double Mea = bl->invBnd->GetMe(iPoint);
 
@@ -155,7 +150,7 @@ 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 == "upper" && iPoint == 168)
+    /*if (bl->name == "upper" && iPoint == 99)
   {
     std::cout << "Regime " << bl->Regime[iPoint] << " \n"
               << "Mea " << Mea << " \n"
@@ -168,7 +163,7 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
               << "delta " << deltaa << " \n"
               << "Cteq " << Cteqa << " \n"
               << "Us " << Usa << std::endl;
-  } */
+  }*/
 
     /* Space part */
 
@@ -227,12 +222,23 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
         timeMatrix(4, 4) = 1.0;
     }
 
+    /*Vector<double, 5> result = timeMatrix.partialPivLu().solve(spaceVector);
+    if (std::isnan(result.norm()))
+    {
+        std::cout << "Point " << iPoint << " det " << timeMatrix.determinant() << std::endl;
+        std::cout << timeMatrix << std::endl;
+        std::cout << spaceVector << std::endl;
+        std::cout << dT_dx << " " << (2 + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx << " " << " cf " << - Cfa / 2 << std::endl;
+        std::cout << " u[0] " << u[0] << " u[1] " << u[1] <<  " u[2] " << u[2] <<  " u[3] " << u[3] << " u[4] " << u[4] << " Usa " << Usa << std::endl;
+        throw;
+    }*/
+
     return timeMatrix.partialPivLu().solve(spaceVector);
 }
 
 double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta)
 {
-    double Ax;
+    double Ax = 0.0;
 
     double Dgr = 0.08;
     double Hk1 = 3.5;
@@ -250,7 +256,7 @@ double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta)
     else
     {
         double Rnorm = (log10(Ret) - (logRet_crit - Dgr)) / (2 * Dgr);
-        double Rfac;
+        double Rfac = 0.0;
         if (Rnorm <= 0)
         {
             Rfac = 0;
diff --git a/dart/src/wViscMesh.cpp b/dart/src/wViscMesh.cpp
index 69546187ae84ef7da3753f76de5761acb72cbc4e..eced1fede5ec0ca11003e46f1a891315e0ad050f 100644
--- a/dart/src/wViscMesh.cpp
+++ b/dart/src/wViscMesh.cpp
@@ -43,9 +43,9 @@ void ViscMesh::SetGlob(std::vector<double> _x, std::vector<double> _y, std::vect
 
 void ViscMesh::ComputeBLcoord()
 {
-    Loc.resize(nMarkers, 0);
-    std::vector<double> dx(nElm,0);
-    alpha.resize(nElm, 0);
+    Loc.resize(nMarkers, 0.0);
+    std::vector<double> dx(nElm,0.0);
+    alpha.resize(nElm, 0.0);
 
     for (size_t iPoint=0; iPoint<nElm; ++iPoint)
     {
diff --git a/dart/src/wViscMesh.h b/dart/src/wViscMesh.h
index decfbe845e2cc54877b9bfd4b5a982c1d5f14be3..21d481f9e6d9add20bd7103da5093bbb96cc9ce1 100644
--- a/dart/src/wViscMesh.h
+++ b/dart/src/wViscMesh.h
@@ -37,7 +37,7 @@ class DART_API ViscMesh
     double Getz(size_t iPoint) const {return z[iPoint];};
     double GetExt(size_t iPoint) const {return Ext[iPoint];};
 
-    unsigned int GetDiffnElm() const {return nMarkers-prevnMarkers;};
+    int GetDiffnElm() const {return nMarkers-prevnMarkers;};
 
     void SetGlob(std::vector<double> _x, std::vector<double> _y, std::vector<double> _z);
     void SetExt(std::vector<double> ExtM) {Ext = ExtM;};
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 5e911afa8923ee602c45fc1e11e379095cc766a4..6bc2ac7dfb3f26daac90a782af0b12784c5c4532 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -102,101 +102,75 @@ int ViscSolver::Run(unsigned int couplIter)
 
     for (size_t iSec=0; iSec<Sections.size(); ++iSec)
     {
-        std::cout << "\n - Sec " << iSec << " ";
+        std::cout << "\n - Sec " << iSec << " stagPtMvmt " << stagPtMvmt[iSec];
 
-        /* Prepare time integration */
+        /* Check for stagnation point movement */
 
-        bool resetSolution = false;
-        if ((couplIter == 0) ||
-             Sections[iSec][0]->mesh->GetDiffnElm()!=0 ||
-             Sections[iSec][1]->mesh->GetDiffnElm()!=0 ||
-             Sections[iSec][2]->mesh->GetDiffnElm()!=0 ||
-            (stagPtMvmt[iSec]))
+        if (couplIter != 0 && (Sections[iSec][0]->mesh->GetDiffnElm()))
         {
-            resetSolution = true;
+            stagPtMvmt[iSec] = true;
+            if (printOn){std::cout << "Stagnation point moved by " << Sections[iSec][0]->mesh->GetDiffnElm() << " elements." << std::endl;}
         }
         else
         {
-            resetSolution = false;
+            stagPtMvmt[iSec] = false;
         }
-        /*std::cout << "resetSolution " << resetSolution << std::endl;
 
-        std::cout << Sections[iSec][0]->mesh->GetDiffnElm()  <<
-             Sections[iSec][1]->mesh->GetDiffnElm()  <<
-             Sections[iSec][2]->mesh->GetDiffnElm() << std::endl;*/
-        resetSolution = true;
-        if (resetSolution)
+        /* Set boundary condition */
+
+        Sections[iSec][0]->xtr = 1.0; /* Upper side initial xtr */
+        Sections[iSec][1]->xtr = 1.0; /* Lower side initial xtr */
+        Sections[iSec][2]->xtr = 0.0; /* Wake initial xtr (full turbulent) */
+
+        Sections[iSec][0]->SetStagBC(Re);
+        Sections[iSec][1]->SetStagBC(Re);
+
+
+        /* Loop over the regions */
+        for (size_t iRegion = 0; iRegion < Sections[iSec].size(); ++iRegion)
         {
-            tSolver->SetCFL0(1);
-            tSolver->SetitFrozenJac(1);
-            tSolver->SetinitSln(true);
+            /* Output region name */
+
+            if (printOn){std::cout << Sections[iSec][iRegion]->name << " ";}
+
+            /* Check if we need to enter safe mode */
 
-            /* Special case where the external flow mesh cannot be used for the interaction law */
-            for (size_t iRegion=0; iRegion<Sections[iSec].size(); ++iRegion)
+            if (couplIter == 0 || Sections[iSec][iRegion]->mesh->GetDiffnElm() != 0 || stagPtMvmt[iSec] == true)
             {
-                std::vector<double> xxExtCurr;
+                tSolver->SetCFL0(1);
+                tSolver->SetitFrozenJac(1);
+                tSolver->SetinitSln(true);
+
+                std::vector<double> xxExtCurr(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.0);
                 for (size_t iPoint=0; iPoint<Sections[iSec][iRegion]->mesh->GetnMarkers(); ++iPoint)
                 {
-                    xxExtCurr.push_back(Sections[iSec][iRegion]->mesh->GetLoc(iPoint));
+                    xxExtCurr[iPoint] = Sections[iSec][iRegion]->mesh->GetLoc(iPoint);
                 }
                 Sections[iSec][iRegion]->mesh->SetExt(xxExtCurr);
-            }
 
-            if (couplIter != 0 && (Sections[iSec][0]->mesh->GetDiffnElm()))
-            {
-                stagPtMvmt[iSec] = true;
-                if (printOn)
+                if (couplIter == 0)
                 {
-                    std::cout << "Stagnation point moved" << std::endl;
+                    std::vector<double> DeltaStarZeros(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.0);
+                    Sections[iSec][iRegion]->invBnd->SetDeltaStar(DeltaStarZeros);
                 }
+
+                if (printOn){std::cout << "restart. ";}
+
             }
             else
             {
-                stagPtMvmt[iSec] = false;
-            }
-            if (printOn)
-            {
-                std::cout << "restart. ";
+                tSolver->SetCFL0(CFL0);
+                tSolver->SetitFrozenJac(5);
+                tSolver->SetinitSln(false);
             }
-        }
 
-        else
-        {
-            tSolver->SetCFL0(CFL0);
-            tSolver->SetitFrozenJac(5);
-            tSolver->SetinitSln(false);
-            stagPtMvmt[iSec] = false;
-            if (printOn)
-            {
-                std::cout << "update. ";
-            }
-        }
+            /* Save number of points in each regions */
 
-        for (size_t iRegion=0; iRegion<Sections[iSec].size(); ++iRegion)
-        { 
             Sections[iSec][iRegion]->mesh->prevnMarkers = Sections[iSec][iRegion]->mesh->GetnMarkers();
-        }
-
-        /* Set boundary condition */
-
-        Sections[iSec][0]->xtr = 1.0; /* Upper side initial xtr */
-        Sections[iSec][1]->xtr = 1.0; /* Lower side initial xtr */
-        Sections[iSec][2]->xtr = 0.0; /* Wake initial xtr (full turbulent) */
-
-        Sections[iSec][0]->SetStagBC(Re);
-        Sections[iSec][1]->SetStagBC(Re);
 
+            /* Reset regime */
 
-        /* Loop over the regions */
-        for (size_t iRegion = 0; iRegion < Sections[iSec].size(); ++iRegion)
-        {
-            if (printOn)
-            {
-                std::cout << Sections[iSec][iRegion]->name << " ";
-            }
             lockTrans = false;
-            
-            /* Reset regime */
             if (iRegion < 2)
             {
                 for (size_t k = 0; k < Sections[iSec][iRegion]->mesh->GetnMarkers(); k++)
@@ -212,10 +186,11 @@ int ViscSolver::Run(unsigned int couplIter)
                 }
 
                 /* Impose wake boundary condition */
-                std::vector<double> UpperCond;
-                std::vector<double> LowerCond;
+                std::vector<double> UpperCond(Sections[iSec][iRegion]->GetnVar(), 0);
+                std::vector<double> LowerCond(Sections[iSec][iRegion]->GetnVar(), 0);
                 UpperCond = std::vector<double>(Sections[iSec][0]->U.end() - Sections[iSec][0]->GetnVar(), Sections[iSec][0]->U.end()); /* Base std::vector constructor to slice */
                 LowerCond = std::vector<double>(Sections[iSec][1]->U.end() - Sections[iSec][1]->GetnVar(), Sections[iSec][1]->U.end()); /* Base std::vector constructor to slice */
+                std::cout << " Size uppercond " << UpperCond.size()  << " Size lowercond " << LowerCond.size() << std::endl;
                 Sections[iSec][iRegion]->SetWakeBC(Re, UpperCond, LowerCond);
                 lockTrans = true;
             }
@@ -229,8 +204,9 @@ int ViscSolver::Run(unsigned int couplIter)
                 tSolver->InitialCondition(iPoint, Sections[iSec][iRegion]);
 
                 /* Time integration */
-                /*if (iRegion == 0 && iPoint <4)
+                if (iRegion == 1 && iPoint == 1)
                 {
+                    std::scientific;
                 std::cout << " iSec " << iSec
                           << " iPoint " << iPoint
                           << " Loc " << Sections[iSec][iRegion]->mesh->GetLoc(iPoint)
@@ -240,17 +216,31 @@ int ViscSolver::Run(unsigned int couplIter)
                           << " DeltaStarExt " << Sections[iSec][iRegion]->invBnd->GetDeltaStar(iPoint)
                           << std::endl;
                 Sections[iSec][iRegion]->PrintSolution(iPoint);
-                }*/
+                }
 
                 pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
 
+                /* if (iRegion == 1 && iPoint == 265)
+                {
+                std::cout << " iSec " << iSec
+                          << " iPoint " << iPoint
+                          << " Loc " << Sections[iSec][iRegion]->mesh->GetLoc(iPoint)
+                          << " Vt " << Sections[iSec][iRegion]->invBnd->GetVt(iPoint)
+                          << " Me " << Sections[iSec][iRegion]->invBnd->GetMe(iPoint)
+                          << " xxExt " << Sections[iSec][iRegion]->mesh->GetExt(iPoint)
+                          << " DeltaStarExt " << Sections[iSec][iRegion]->invBnd->GetDeltaStar(iPoint)
+                          << std::scientific;
+                Sections[iSec][iRegion]->PrintSolution(iPoint);
+                } */
+
+                /* Unsucessfull convergence output */
+
                 if (pointExitCode != 0)
                 {
-                    if (printOn)
-                    {
-                        std::cout << "Point " << iPoint << ": Convergence terminated with exitcode " << pointExitCode << std::endl;
-                    }
+                    if (printOn){std::cout << "Point " << iPoint << ": Convergence terminated with exitcode " << pointExitCode << std::endl;}
+                    stagPtMvmt[iSec] = true;
                     solverExitCode = -1; // Convergence failed
+                    
                 }
 
                 /* Transition */
@@ -280,10 +270,11 @@ int ViscSolver::Run(unsigned int couplIter)
     } // End for iSections
 
     ComputeDrag(Sections[0]);
-    if (printOn)
+    if (solverExitCode != 0)
     {
-        std::cout << "\n" << std::endl;
+        std::cout << "some points did not converge. ";
     }
+    if (printOn){std::cout << "\n" << std::endl;}
 
     return solverExitCode; // exit code
 }
@@ -316,7 +307,7 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
     unsigned int nVar = bl->GetnVar();
 
     // Save laminar solution.
-    std::vector<double> lamSol(nVar);
+    std::vector<double> lamSol(nVar, 0.0);
     for (size_t k = 0; k < lamSol.size(); ++k)
     {
         lamSol[k] = bl->U[iPoint * nVar + k];
@@ -353,7 +344,7 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
 
     if (exitCode != 0)
     {
-        std::cout << "Warning level 3: Transition point turbulent computation terminated with exit code " << exitCode << std::endl;
+        std::cout << "Warning: Transition point turbulent computation terminated with exit code " << exitCode << std::endl;
     }
 
     /* Average both solutions (Now solution @ iPoint is the turbulent solution) */
@@ -515,7 +506,7 @@ void ViscSolver::PrintBanner()
     std::cout << "#                                                   ///,        ////           #" << std::endl;
     std::cout << "#       Paul Dechamps                               \\  /,      /  >.           #" << std::endl;
     std::cout << "#       ULiege 2022                                  \\  /,   _/  /.            #" << std::endl;
-    std::cout << "#       v0.1.0 - Memphis                              \\_  /_/   /.             #" << std::endl;
+    std::cout << "#                                                     \\_  /_/   /.             #" << std::endl;
     std::cout << "#                                                      \\__/_   <               #" << std::endl;
     std::cout << "#                                                      /<<< \\_\\_               #" << std::endl;
     std::cout << "#   ██╗  ██╗ ██████╗ ██████╗ ██╗   ██╗███████╗        /,)^>>_._ \\              #" << std::endl;
diff --git a/dart/tests/bli3.py b/dart/tests/bli3.py
index 5ab3cfe3b37eba033d85e4a003867d9be41d1f81..8083eee7f8e359e237658011d7a5e08e70e7a417 100644
--- a/dart/tests/bli3.py
+++ b/dart/tests/bli3.py
@@ -31,7 +31,7 @@ import dart.default as floD
 
 import dart.pyVII.vUtils as viscU
 import dart.pyVII.vCoupler as viscC
-import dart.pyVII.vInterpolator as vInterpol
+import dart.pyVII.vInterpolator2 as vInterpol
 
 import fwk
 from fwk.testing import *
@@ -74,7 +74,8 @@ def main():
     tms['pre'].stop()
 
     # solve problem
-    iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw[0], blw[1], nSections, dim)
+    config = {'nDim': dim, 'nSections':nSections, 'nPoints':[500, 25], 'eps':0.05, 'k':[5, 5]}
+    iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw[0], blw[1], config)
     vSolver = viscU.initBL(Re, M_inf, CFL0, nSections)
     #iSolverAPI = viscAPI.dartAPI(floD.newton(pbl), bnd, wk, nSections, vInterp)
     coupler = viscC.Coupler(iSolverAPI, vSolver)
diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py
index 7544e5eb81e5226c9b874296c9a72d6a74364d5f..d827250aafd2650a4a28dd2cce477b95c2811ba5 100644
--- a/dart/tests/bli_lowLift.py
+++ b/dart/tests/bli_lowLift.py
@@ -34,12 +34,14 @@
 # This test is provided to ensure that the solver works properly.
 # Mesh refinement may have to be performed to obtain physical results.
 
+from tkinter import NS
 import dart.default as floD
 
 import dart.pyVII.vUtils as viscU
 import dart.pyVII.vCoupler as viscC
-import dart.pyVII.vInterpolator as vInterpol
+import dart.pyVII.vInterpolator2 as vInterpol
 
+from matplotlib import pyplot as plt
 import numpy as np
 
 import tbox
@@ -58,7 +60,7 @@ def main():
 
     # define flow variables
     Re = 1e7
-    alpha = 5.*math.pi/180
+    alpha = 12.*math.pi/180
     #alphaSeq = np.arange(-5, 5.5, 0.5)
     M_inf = 0.
 
@@ -87,13 +89,22 @@ def main():
 
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
     tms['solver'].start()
-    
-    iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw[0], blw[1], nSections, dim)
-    vSolver = viscU.initBL(Re, M_inf, CFL0, nSections)
-    coupler = viscC.Coupler(iSolverAPI, vSolver)
+
+    config = {'nDim': dim, 'nSections':nSections, 'nPoints':[500, 50], 'eps':0.02}
+
+    #iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw[0], blw[1], nSections, dim)
+    dragIter = [np.zeros(0) for i in range(1)]
+    for i in range(len(dragIter)):
+        iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw[0], blw[1], config)
+        vSolver = viscU.initBL(Re, M_inf, CFL0, nSections)
+        coupler = viscC.Coupler(iSolverAPI, vSolver)
 
     #coupler.RunPolar(alphaSeq)
-    coupler.Run()
+        dragIter[i] = coupler.Run()
+        plt.plot(dragIter[i])
+    plt.title('{:.0f} real. {:.2f} deg AoA'.format(len(dragIter), alpha*180/np.pi))
+    plt.xlabel('Iters')
+    plt.ylabel('$c_d$')
     tms['solver'].stop()
 
     # extract Cp
@@ -112,6 +123,8 @@ def main():
     print('SOLVERS statistics')
     print(coupler.tms)
     print(iSolverAPI.tms)
+    plt.show()
+
 
     """# visualize solution and plot results
     #floD.initViewer(pbl)
@@ -126,7 +139,7 @@ def main():
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
     if Re == 1e7 and M_inf == 0. and alpha == 2.*math.pi/180:
-        tests.add(CTest('Cl', iSolverAPI.GetCl(), 0.2208, 5e-3))
+        tests.add(CTest('Cl', iSolverAPI.GetCl(), 0.2271, 5e-3)) # Xfoil 0.2271
         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', vSolver.Getxtr(0,0), 0.20, 0.03, forceabs=True))