diff --git a/dart/pyVII/vCoupler.py b/dart/pyVII/vCoupler.py
index 59e5ba09cfec3cf80d564e2b85f0fe17c700c55a..9de6006f9bc060cabee94e063d0a2177b9316dc2 100644
--- a/dart/pyVII/vCoupler.py
+++ b/dart/pyVII/vCoupler.py
@@ -36,7 +36,7 @@ class Coupler:
     self.iSolver=iSolver
     self.vSolver=vSolver
 
-    self.maxCouplIter = 1
+    self.maxCouplIter = 100
     self.couplTol = 1e-4
 
     self.tms=fwk.Timers()
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index 03a01b104e79926bbb0f3136e4a7153b53cbe7e3..4e3e81e05bff7b2a0c0164f5170f1bc8d5dc63a3 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -10,7 +10,7 @@ using namespace dart;
 
 BLRegion::BLRegion(double _xtrF)
 {
-    if (_xtrF < 0.0 && _xtrF != -1.0)
+    if (_xtrF < 0. && _xtrF != -1.)
     {
         std::cout << "Unrecognized transition location " << _xtrF << std::endl;
         throw;
@@ -45,18 +45,18 @@ void BLRegion::SetMesh(std::vector<double> x,
 
     if (mesh->GetDiffnElm()!=0)
     {
-        U.resize(nVar * nMarkers, 0);
-        Ret.resize(nMarkers, 0);
-        Cd.resize(nMarkers);
-        Cf.resize(nMarkers,0);
-        Hstar.resize(nMarkers,0);
-        Hstar2.resize(nMarkers,0);
-        Hk.resize(nMarkers,0);
-        Cteq.resize(nMarkers,0);
-        Us.resize(nMarkers,0);
-        DeltaStar.resize(nMarkers,0);
-        Delta.resize(nMarkers,0);
-        Regime.resize(nMarkers,0);
+        U.resize(nVar * nMarkers, 0.);
+        Ret.resize(nMarkers, 0.);
+        Cd.resize(nMarkers, 0.);
+        Cf.resize(nMarkers,0.);
+        Hstar.resize(nMarkers,0.);
+        Hstar2.resize(nMarkers,0.);
+        Hk.resize(nMarkers,0.);
+        Cteq.resize(nMarkers,0.);
+        Us.resize(nMarkers,0.);
+        DeltaStar.resize(nMarkers,0.);
+        Delta.resize(nMarkers,0.);
+        Regime.resize(nMarkers,0.);
     }
 } // End SetMesh
 
diff --git a/dart/src/wClosures.cpp b/dart/src/wClosures.cpp
index 6497f22565e06567e18014af70b1c18ba7e07d76..669c38ea985f4f925e5d716fb0a9337a56ceaf06 100644
--- a/dart/src/wClosures.cpp
+++ b/dart/src/wClosures.cpp
@@ -39,7 +39,7 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
 
     double delta = std::min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta);
 
-    double Hstar = 0.0;
+    double Hstar = 0.;
     double Hks = Hk - 4.35;
     if (Hk <= 4.35)
     {
@@ -54,8 +54,8 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
 
     /* Friction coefficient */
 
-    double tmp = 0.0;
-    double Cf = 0.0;
+    double tmp = 0.;
+    double Cf = 0.;
     if (Hk < 5.5)
     {
         tmp = (5.5-Hk)*(5.5-Hk)*(5.5-Hk) / (Hk+1.0);
@@ -69,7 +69,7 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
 
     /* Dissipation coefficient */
 
-    double Cd = 0.0;
+    double Cd = 0.;
     if (Hk < 4)
     {
         Cd = ( 0.00205*std::pow(4.0-Hk, 5.5) + 0.207 ) * (Hstar/(2*Ret));
@@ -89,10 +89,10 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
         Cf = 0;
     }
 
-    double Us = 0.0;
-    double Cteq = 0.0;
+    double Us = 0.;
+    double Cteq = 0.;
 
-    std::vector<double> ClosureVars(8, 0.0);
+    std::vector<double> ClosureVars(8, 0.);
 
     ClosureVars[0] = Hstar;
     ClosureVars[1] = Hstar2;
@@ -135,7 +135,7 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
     double gamma = 1.4;
     double Fc = std::sqrt(1+0.5*gamma*Me*Me);
 
-    double H0 = 0.0;
+    double H0 = 0.;
     if (Ret <= 400)
     {
         H0 = 4.0;
@@ -149,7 +149,7 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
         Ret = 200;
     }
 
-    double Hstar = 0.0;
+    double Hstar = 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;
@@ -179,17 +179,17 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
     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 Hkc = 0.;
+    double Cdw = 0.;
+    double Cdd = 0.;
+    double Cdl = 0.;
 
-    double Hmin = 0.0;
-    double Fl = 0.0;
-    double Dfac = 0.0;
+    double Hmin = 0.;
+    double Fl = 0.;
+    double Dfac = 0.;
 
-    double Cf = 0.0;
-    double Cd = 0.0;
+    double Cf = 0.;
+    double Cd = 0.;
 
     double n = 1.0;
     
@@ -232,7 +232,7 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
     Cd = n*(Cdw+ Cdd + Cdl);
     // Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk))**2)
 
-    std::vector<double> ClosureVars(8, 0.0);
+    std::vector<double> ClosureVars(8, 0.);
 
     ClosureVars[0] = Hstar;
     ClosureVars[1] = Hstar2;
diff --git a/dart/src/wInvBnd.h b/dart/src/wInvBnd.h
index e1aa68d25770072a0a13aeeaa9040de99ef6a890..c14f2e19e797babd547d75fc5f952e2dd875aaf9 100644
--- a/dart/src/wInvBnd.h
+++ b/dart/src/wInvBnd.h
@@ -31,11 +31,11 @@ class DART_API InvBnd
     double GetDeltaStar(size_t iPoint) const {return DeltaStar[iPoint];};
     double GetVtOld(size_t iPoint) const {return VtOld[iPoint];};
 
-    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;};
-    void SetVtOld(std::vector<double> _VtOld) {VtOld.resize(_VtOld.size(), 0.0); VtOld = _VtOld;};
+    void SetMe(std::vector<double> _Me) {Me.resize(_Me.size(), 0.); Me = _Me;};
+    void SetVt(std::vector<double> _Vt) {Vt.resize(_Vt.size(), 0.); Vt = _Vt;};
+    void SetRhoe(std::vector<double> _Rhoe) {Rhoe.resize(_Rhoe.size(), 0.); Rhoe = _Rhoe;};
+    void SetDeltaStar(std::vector<double> _DeltaStar) {DeltaStar.resize(_DeltaStar.size(), 0.); DeltaStar = _DeltaStar;};
+    void SetVtOld(std::vector<double> _VtOld) {VtOld.resize(_VtOld.size(), 0.); VtOld = _VtOld;};
 
 };
 } // namespace dart
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index 265f5c0961c52f3d4a2868f53014d798e2eb5a55..ac2f4cd9214938f805ee1069b07509d8562c56c5 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -67,7 +67,7 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
 
     /* Save initial condition */
 
-    std::vector<double> Sln0(nVar, 0.0);
+    std::vector<double> Sln0(nVar, 0.);
     for (size_t i = 0; i < nVar; ++i)
     {
         Sln0[i] = bl->U[iPoint * nVar + i];
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index 8982b4f96361a786ee07826c4907161870210b9f..5d0e0a19e3eca0f44194849127734cde70e7cb52 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -26,7 +26,7 @@ ViscFluxes::~ViscFluxes()
 Vector<double, 5> ViscFluxes::ComputeFluxes(size_t iPoint, BLRegion *bl)
 {
     unsigned int nVar = bl->GetnVar();
-    std::vector<double> u(nVar, 0.0);
+    std::vector<double> u(nVar, 0.);
     for (size_t k = 0; k < nVar; ++k)
     {
         u[k] = bl->U[iPoint * nVar + k];
@@ -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(nVar, 0.0);
+    std::vector<double> uUp(nVar, 0.);
     for (size_t k = 0; k < nVar; ++k)
     {
         uUp[k] = bl->U[iPoint * nVar + k];
     } // End for
 
-    double varSave = 0.0;
+    double varSave = 0.;
     for (size_t iVar = 0; iVar < nVar; ++iVar)
     {
         varSave = uUp[iVar];
@@ -188,31 +188,31 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
 
     timeMatrix(0, 0) = u[1] / u[3];
     timeMatrix(0, 1) = u[0] / u[3];
-    timeMatrix(0, 2) = 0.0;
+    timeMatrix(0, 2) = 0.;
     timeMatrix(0, 3) = u[0] * u[1] / (u[3] * u[3]);
-    timeMatrix(0, 4) = 0.0;
+    timeMatrix(0, 4) = 0.;
 
     timeMatrix(1, 0) = (1 + u[1] * (1 - Hstara)) / u[3];
     timeMatrix(1, 1) = (1 - Hstara) * u[0] / u[3];
-    timeMatrix(1, 2) = 0.0;
+    timeMatrix(1, 2) = 0.;
     timeMatrix(1, 3) = (2 - Hstara * u[1]) * u[0] / (u[3] * u[3]);
-    timeMatrix(1, 4) = 0.0;
+    timeMatrix(1, 4) = 0.;
 
-    timeMatrix(2, 0) = 0.0;
-    timeMatrix(2, 1) = 0.0;
-    timeMatrix(2, 2) = 1.0;
-    timeMatrix(2, 3) = 0.0;
-    timeMatrix(2, 4) = 0.0;
+    timeMatrix(2, 0) = 0.;
+    timeMatrix(2, 1) = 0.;
+    timeMatrix(2, 2) = 1.;
+    timeMatrix(2, 3) = 0.;
+    timeMatrix(2, 4) = 0.;
 
     timeMatrix(3, 0) = -c * u[1];
     timeMatrix(3, 1) = -c * u[0];
-    timeMatrix(3, 2) = 0.0;
-    timeMatrix(3, 3) = 1.0;
-    timeMatrix(3, 4) = 0.0;
+    timeMatrix(3, 2) = 0.;
+    timeMatrix(3, 3) = 1.;
+    timeMatrix(3, 4) = 0.;
 
-    timeMatrix(4, 0) = 0.0;
-    timeMatrix(4, 1) = 0.0;
-    timeMatrix(4, 2) = 0.0;
+    timeMatrix(4, 0) = 0.;
+    timeMatrix(4, 1) = 0.;
+    timeMatrix(4, 2) = 0.;
 
     if (bl->Regime[iPoint] == 1)
     {
@@ -221,7 +221,7 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
     }
     else
     {
-        timeMatrix(4, 3) = 0.0;
+        timeMatrix(4, 3) = 0.;
         timeMatrix(4, 4) = 1.0;
     }
 
@@ -241,7 +241,7 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
 
 double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta) const
 {
-    double Ax = 0.0;
+    double Ax = 0.;
 
     double Dgr = 0.08;
     double Hk1 = 3.5;
@@ -259,7 +259,7 @@ double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta) con
     else
     {
         double Rnorm = (log10(Ret) - (logRet_crit - Dgr)) / (2 * Dgr);
-        double Rfac = 0.0;
+        double Rfac = 0.;
         if (Rnorm <= 0)
         {
             Rfac = 0;
@@ -281,7 +281,7 @@ double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta) con
         if (Hk > Hk1)
         {
             double Hnorm = (Hk - Hk1) / (Hk2 - Hk1);
-            double Hfac;
+            double Hfac = 0.;
             if (Hnorm >= 1)
             {
                 Hfac = 1;
@@ -299,7 +299,7 @@ double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta) con
                 Ax2 = 0;
             }
             Ax = Hfac * Ax2 + (1 - Hfac) * Ax1;
-            Ax = std::max(Ax, 0.0);
+            Ax = std::max(Ax, 0.);
         }
     }
     return Ax;
diff --git a/dart/src/wViscMesh.cpp b/dart/src/wViscMesh.cpp
index eced1fede5ec0ca11003e46f1a891315e0ad050f..0a36a373895dc3037d2505baea171e5c4089cc32 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.0);
-    std::vector<double> dx(nElm,0.0);
-    alpha.resize(nElm, 0.0);
+    Loc.resize(nMarkers, 0.);
+    std::vector<double> dx(nElm,0.);
+    alpha.resize(nElm, 0.);
 
     for (size_t iPoint=0; iPoint<nElm; ++iPoint)
     {
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 24846a42f97e595c21b3f7629ce714a620e692e4..7297dfcc8be0f5f254fa51929ad3c11c5f36ed50 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -38,9 +38,9 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec
     Re = _Re;     /* Freestream Reynolds number */
 
     CFL0 = _CFL0;
-    Cdt = 0.0;
-    Cdf = 0.0;
-    Cdp = 0.0;
+    Cdt = 0.;
+    Cdf = 0.;
+    Cdp = 0.;
 
     solverExitCode = 0;
 
@@ -122,9 +122,9 @@ int ViscSolver::Run(unsigned int couplIter)
 
         /* 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]->xtr = 1.; /* Upper side initial xtr */
+        Sections[iSec][1]->xtr = 1.; /* Lower side initial xtr */
+        Sections[iSec][2]->xtr = 0.; /* Wake initial xtr (full turbulent) */
 
         Sections[iSec][0]->SetStagBC(Re);
         Sections[iSec][1]->SetStagBC(Re);
@@ -153,8 +153,8 @@ int ViscSolver::Run(unsigned int couplIter)
                 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)
+                std::vector<double> xxExtCurr(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.);
+                for (size_t iPoint=0; iPoint<xxExtCurr.size(); ++iPoint)
                 {
                     xxExtCurr[iPoint] = Sections[iSec][iRegion]->mesh->GetLoc(iPoint);
                 }
@@ -495,7 +495,7 @@ void ViscSolver::SetDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double
 
 std::vector<double> ViscSolver::ExtractxCoord(size_t iSec, size_t iRegion) const
 {
-    std::vector<double> xCoord(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.0);
+    std::vector<double> xCoord(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.);
     for (size_t iPoint=0; iPoint<xCoord.size(); ++iPoint)
     {
         xCoord[iPoint] = Sections[iSec][iRegion]->mesh->Getx(iPoint);
@@ -509,7 +509,7 @@ std::vector<double> ViscSolver::ExtractBlowingVel(size_t iSec, size_t iRegion) c
 }
 std::vector<double> ViscSolver::Extractxx(size_t iSec, size_t iRegion) const
 {   
-    std::vector<double> xx(Sections[iSec][iRegion]->mesh->GetnMarkers(),0.0);
+    std::vector<double> xx(Sections[iSec][iRegion]->mesh->GetnMarkers(),0.);
     for (size_t iPoint=0; iPoint<xx.size(); ++iPoint)
     {
         xx[iPoint] = Sections[iSec][iRegion]->mesh->GetLoc(iPoint);