diff --git a/blast/blCoupler.py b/blast/blCoupler.py
index cc57b542a263e84764a677ceef884c1993945b1f..a52dc6c2311c909dd4267f30191207104c378253 100644
--- a/blast/blCoupler.py
+++ b/blast/blCoupler.py
@@ -25,15 +25,21 @@ import numpy as np
 import blast.blUtils as vutils
 
 class Coupler:
-    def __init__(self, iSolverAPI, vSolver, _maxCouplIter=150, _couplTol=1e-4, _iterPrint=1, _resetInv=False, sfx=''):
+    def __init__(self, iSolverAPI, vSolver, vconfig):
+
         self.isol = iSolverAPI
         self.vsol = vSolver
-        self.maxIter = _maxCouplIter
-        self.tol = _couplTol
-        self.resetInviscid = _resetInv
-        self.iterPrint = _iterPrint if self.isol.getVerbose() == 0 and self.vsol.verbose == 0 else 1
+        self.maxIter = vconfig.get('couplIter', 50)
+        self.tol = vconfig.get('couplTol', 1e-4)
+        self.resetInviscid = vconfig.get('resetInv', True)
+        self.iterPrint = vconfig.get('iterPrint', 1)
+        self.iterSaveRestart = vconfig.get('iterSaveRestart', float('inf'))
+        self.saveCouplingIters = vconfig.get('saveCouplingIters', False)
+        self.iterPrint = self.iterPrint if self.isol.getVerbose() == 0 and self.vsol.verbose == 0 else 1
         self.tms = fwk.Timers()
-        self.filesfx = sfx
+        self.filesfx = vconfig.get('filesfx', '')
+        self.restart = vconfig.get('restart_solution', False)
+
         print('')
         print(' ######################################################## ')
         print('|            ___  __   ___   _________________           |')
@@ -59,18 +65,20 @@ class Coupler:
         print(f'{"Verbosity level:":<30s} {self.isol.getVerbose():<2.0f}')
         print('')
 
-    def run(self, write=True):
+    def run(self, write=True, doit=False):
         # Aerodynamic coefficients.
         aeroCoeffs = {'Cl':[], 'Cd': [], 'Cdwake': []}
 
         # Convergence parameters.
         couplIter = 0
-        cdPrev = 0.0
+        if self.restart:
+            cdPrev = self.isol.cdPrev
+        else:
+            cdPrev = 0.0
 
+        self.isol.getBlowingBoundary(couplIter)
         while couplIter < self.maxIter:
-            # Impose blowing boundary condition in the inviscid solver.
             self.tms['processing'].start()
-            self.isol.getBlowingBoundary(couplIter)
             self.isol.interpolate('v2i')
             self.isol.setBlowingVelocity()
             self.tms['processing'].stop()
@@ -102,6 +110,8 @@ class Coupler:
             vEc = self.vsol.run()
             self.tms['viscous'].stop()
 
+            self.isol.getBlowingBoundary(couplIter+1)
+
             aeroCoeffs['Cl'].append(self.isol.getCl())
             aeroCoeffs['Cd'].append(self.isol.getCd() + self.vsol.Cdf)
             aeroCoeffs['Cdwake'].append(self.vsol.Cdt)
@@ -111,6 +121,11 @@ class Coupler:
             #cd = self.vsol.Cdt if self.vsol.Cdt != 0 else self.vsol.Cdf + self.isol.getCd()
             error = abs((cd - cdPrev) / cd) if cd != 0 else 1
 
+            # Save restart file
+            if couplIter % self.iterSaveRestart == 0:
+                _sfx = '_'+str(couplIter) if self.saveCouplingIters else ''
+                self.isol.save_restart(cd, sfx=_sfx)
+
             if error <= self.tol:
                 print(ccolors.ANSI_GREEN, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>6.4f} | {:>5.0f} {:>5.0f} | {:>6.3f}\n'.format(couplIter, self.isol.getCl(), self.isol.getTotalDrag(), self.vsol.Cdt, self.vsol.bodies[0].getAvgxtr(0), self.vsol.bodies[0].getAvgxtr(1), iEc, vEc, np.log10(error)), ccolors.ANSI_RESET)
                 if iEc != 0 or vEc != 0:
@@ -156,12 +171,10 @@ class Coupler:
         # Blowing velocity
         for ibody in range(self.isol.getnBodies()):
             for b in self.isol.iBnd[ibody]:
-                for i in range(len(b.blowingVel)):
-                    b.blowingVel[i] = 0.
+                b.setBlowingVelocity(np.zeros(b.getnElms()))
             for sec in self.isol.vBnd[ibody]:
                 for side in sec:
-                    for i in range(len(side.blowingVel)):
-                        side.blowingVel[i] = 0.
+                    side.setBlowingVelocity(np.zeros(side.getnElms()))
             self.isol.setBlowingVelocity()
 
         for ibody, body in enumerate(self.vsol.bodies):
diff --git a/blast/blUtils.py b/blast/blUtils.py
index b7edca099856014880da624d5402128223663948..2b4e24858dd22d491f3423f43bc7c183374470d5 100644
--- a/blast/blUtils.py
+++ b/blast/blUtils.py
@@ -112,11 +112,9 @@ def initBlast(iconfig, vconfig, iSolver='DART', task='analysis'):
     vsol = initBL(vconfig['Re'], isol.getMinf(), vconfig['CFL0'], vconfig['sections'], types, vconfig['spans'], vconfig['xtrF'], isol.getVerbose())
     isol.addViscousSolver(vsol)
 
-
-
     # Coupler
     import blast.blCoupler as blastCoupler
-    coupler = blastCoupler.Coupler(isol, vsol, _maxCouplIter=vconfig['couplIter'], _couplTol=vconfig['couplTol'], _iterPrint=vconfig['iterPrint'], _resetInv=vconfig['resetInv'], sfx=vconfig['sfx'])
+    coupler = blastCoupler.Coupler(isol, vsol, vconfig)
     return coupler, isol, vsol
 
 def mesh(file, pars):