diff --git a/sdpm/api/core.py b/sdpm/api/core.py
index 96caf8a2cd75a569e3dc088f563cfbcc87947af5..d7d6b38764efba6dad38097667848bda0192dd1b 100644
--- a/sdpm/api/core.py
+++ b/sdpm/api/core.py
@@ -17,7 +17,7 @@
 ## Initialize SDPM
 # Adrien Crovato
 
-def initSdpm(cfg, use_ad=False):
+def init_sdpm(cfg, use_ad=False):
     """Initialize SDPM components
 
     Parameters
@@ -82,7 +82,8 @@ def initSdpm(cfg, use_ad=False):
     _pbl.addBody(_bdy)
 
     # Solver
-    _sol = sdpm.Solver(_pbl)
+    vrb = cfg['Verb_lvl'] if 'Verb_lvl' in cfg else 1
+    _sol = sdpm.Solver(_pbl, vrb)
     _adj = sdpm.Adjoint(_sol) if use_ad else None
 
     # Return
diff --git a/sdpm/api/om.py b/sdpm/api/om.py
index 267f8668eb24d2d64bb469222461bd7b65e8b1a7..e052755b954b14a03eac0835fdce327ab701265a 100644
--- a/sdpm/api/om.py
+++ b/sdpm/api/om.py
@@ -107,15 +107,15 @@ class SdpmBuilder(Builder):
         use_ad : bool, optional
             Whether to use AD within SDPM (default: False)
         """
-        from .core import initSdpm
-        self.__sdpm = initSdpm(cfg, use_ad)
+        from .core import init_sdpm
+        self.__sdpm = init_sdpm(cfg, use_ad)
 
     def get_mesh(self):
         """Return OpenMDAO component to get the initial surface mesh coordinates
         """
         return SdpmMesh(bdy=self.__sdpm['bdy'])
 
-    def get_solver(self):
+    def get_solver(self, scenario_name=''):
         """Return OpenMDAO component containing the solver
         """
         return SdpmSolver(nf=self.__sdpm['n_f'], nm=self.__sdpm['n_m'], bdy=self.__sdpm['bdy'], sol=self.__sdpm['sol'], adj=self.__sdpm['adj'])
diff --git a/sdpm/cases/agard.py b/sdpm/cases/agard.py
index 742431f723c372cd11ea55ee682b5ddbb5e2350e..715427ab696967a5cf9b2fe28142f708f6c14d96 100644
--- a/sdpm/cases/agard.py
+++ b/sdpm/cases/agard.py
@@ -34,7 +34,7 @@ def get():
         # Unsteady motion
         'Num_wake_div': wnc, # number of cells per chord length in the wake
         'Frequencies': [0.01, 0.1, 0.3, 0.5], # list of reduced frequencies
-        'Modes': build_fpath('models/agard445_modes', __file__, 1), # input file containing the modes (optional)
+        'Modes': build_fpath('models/agard445_modes.csv', __file__, 1), # input file containing the modes (optional)
         'Num_modes': 4, # number of modes
         # Geometry
         'Symmetry': True, # whether only half (symmetric) model is provided
@@ -48,9 +48,9 @@ def get():
 
 def main():
     # Init SDPM
-    from sdpm.api.core import initSdpm
+    from sdpm.api.core import init_sdpm
     import numpy as np
-    sdpm = initSdpm(get())
+    sdpm = init_sdpm(get())
     # Run the solver
     sdpm['sol'].update()
     sdpm['sol'].run()
diff --git a/sdpm/src/sdpmSolver.cpp b/sdpm/src/sdpmSolver.cpp
index 1ee540b4c8a336dacb2bc75eaf6652cae50a5bf0..17cd2ddf059a2a18aa46abf66d35a11677c16956 100644
--- a/sdpm/src/sdpmSolver.cpp
+++ b/sdpm/src/sdpmSolver.cpp
@@ -29,7 +29,7 @@
 
 using namespace sdpm;
 
-Solver::Solver(Problem &pbl) : _pbl(pbl), _cl(0), _cd(0), _cs(0), _cm(0)
+Solver::Solver(Problem &pbl, int vrb) : _pbl(pbl), _vrb(vrb), _cl(0), _cd(0), _cs(0), _cm(0)
 {
     // Say hi
     std::cout << std::endl;
@@ -247,29 +247,34 @@ void Solver::update()
  */
 void Solver::run()
 {
-    std::cout << "Running SDPM solver...\n"
-              << "-- Mean flow" << std::endl;
-    std::cout << "computing sources... " << std::flush;
+    std::cout << "Running SDPM solver..." << std::flush;
+    if (_vrb > 0)
+        std::cout << "\n-- Mean flow\n"
+                  << "computing sources... " << std::flush;
     sdpmVector3d ui = _pbl.getDragDir(); // freestream velocity
     sdpmVectorXd etau(_nP);
     for (auto body : _pbl.getBodies())
         for (auto e : body->getElements())
             etau(_rows[e]) = ui.dot(e->getCompressibleNormal().cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.)));
-    std::cout << "done." << std::endl;
-    std::cout << "solving for doublets... " << std::flush;
+    if (_vrb > 0)
+        std::cout << "done.\n"
+                  << "solving for doublets... " << std::flush;
     sdpmVectorXd rhs = _B0 * etau;
     sdpmVectorXd emu = solve(_A0, rhs);
-    std::cout << "done." << std::endl;
-    std::cout << "computing flow and loads on bodies... " << std::flush;
+    if (_vrb > 0)
+        std::cout << "done.\n"
+                  << "computing flow and loads on bodies... " << std::flush;
     post(etau, emu);
-    std::cout << "done." << std::endl;
+    if (_vrb > 0)
+        std::cout << "done." << std::endl;
 
     for (size_t imd = 0; imd < _pbl.getNumModes(); ++imd)
     {
         for (size_t ifq = 0; ifq < _pbl.getFrequencies().size(); ++ifq)
         {
-            std::cout << "-- Mode #" << imd << ", frequency #" << ifq << std::endl;
-            std::cout << "computing sources... " << std::flush;
+            if (_vrb > 0)
+                std::cout << "-- Mode #" << imd << ", frequency #" << ifq << "\n"
+                          << "computing sources... " << std::flush;
             sdpmVectorXcd etau1(_nP);
             for (auto body : _pbl.getBodies())
             {
@@ -279,9 +284,10 @@ void Solver::run()
                 for (size_t i = 0; i < elems.size(); ++i)
                     etau1(_rows[elems[i]]) = (uc0[i] + uc1[i] * sdpmComplex(0., 1.) * _pbl.getFrequencies()[ifq]).conjugate().dot(elems[i]->getCompressibleNormal().cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.)));
             }
-            std::cout << "done." << std::endl;
+            if (_vrb > 0)
+                std::cout << "done.\n"
+                          << "solving for doublets... " << std::flush;
             // Split complex set of equations into real set (twice the size)
-            std::cout << "solving for doublets... " << std::flush;
             sdpmMatrixXd A(2 * _nP, 2 * _nP);
             A.block(0, 0, _nP, _nP) = _A[ifq].real();
             A.block(_nP, _nP, _nP, _nP) = _A[ifq].real();
@@ -292,19 +298,26 @@ void Solver::run()
             rhs1.block(_nP, 0, _nP, 1) = (_B[ifq] * etau1).imag();
             sdpmVectorXd mu = solve(A, rhs1);
             sdpmVectorXcd emu1 = mu.block(0, 0, _nP, 1) + sdpmComplex(0, 1) * mu.block(_nP, 0, _nP, 1);
-            std::cout << "done." << std::endl;
-            std::cout << "computing flow and loads on bodies... " << std::flush;
+            if (_vrb > 0)
+                std::cout << "done.\n"
+                          << "computing flow and loads on bodies... " << std::flush;
             post(imd, ifq, etau1, emu1);
-            std::cout << "done." << std::endl;
+            if (_vrb > 0)
+                std::cout << "done." << std::endl;
         }
     }
-    std::cout << "-- Generalized Aerodynamic Force matrices" << std::endl;
+    if (_vrb > 0)
+        std::cout << "-- Generalized Aerodynamic Force matrices" << std::endl;
     for (size_t ifq = 0; ifq < _pbl.getFrequencies().size(); ++ifq)
     {
-        std::cout << "computing GAF at frequency #" << ifq << "... " << std::flush;
+        if (_vrb > 0)
+            std::cout << "computing GAF at frequency #" << ifq << "... " << std::flush;
         computeGAF(ifq);
-        std::cout << "done." << std::endl;
+        if (_vrb > 0)
+            std::cout << "done." << std::endl;
     }
+    if (_vrb == 0)
+        std::cout << " done.";
     std::cout << std::endl;
 }
 
diff --git a/sdpm/src/sdpmSolver.h b/sdpm/src/sdpmSolver.h
index 313f019818401410c44ba5d07f8f623c99a8c277..ee6d5649538c1c48bd155049a4c1097d4cd7280b 100644
--- a/sdpm/src/sdpmSolver.h
+++ b/sdpm/src/sdpmSolver.h
@@ -36,6 +36,7 @@ class SDPM_API Solver : public sdpmObject
 protected:
     // Problem
     Problem &_pbl; ///< problem definition
+    int _vrb;      ///< verbosity level
     // AICs
     size_t _nP;                     ///< number of body panels
     std::map<Element *, int> _rows; ///< map linking an element to a matrix cell
@@ -72,7 +73,7 @@ private:
     std::vector<sdpmMatrixXd> _q1Im; ///< imaginary part of GAF matrices
 
 public:
-    Solver(Problem &pbl);
+    Solver(Problem &pbl, int vrb = 1);
     virtual ~Solver();
 
     std::vector<sdpmDouble> const &getPressure() const { return _cp; }
diff --git a/sdpm/tests/agard.py b/sdpm/tests/agard.py
index 42b924eb6e1c319eca6b0ba72d34fceb2cdefeb6..b70f59bbe328cdafc354db07f4c5bd0902fce324 100644
--- a/sdpm/tests/agard.py
+++ b/sdpm/tests/agard.py
@@ -50,7 +50,7 @@ def main():
     pbl.setFreestream(aoa, 0., minf)
     pbl.setUnsteady(n_modes, [kred])
     wing = sdpm.Body(msh, 'wing', 'te', crd, wnc, n_modes, 1)
-    x_modes, amp_modes = interpolate_modes(build_fpath('models/agard445_modes', __file__, 1), wing.getNodes())
+    x_modes, amp_modes = interpolate_modes(build_fpath('models/agard445_modes.csv', __file__, 1), wing.getNodes())
     for i in range(n_modes):
         wing.addMotion()
         wing.setFlexibleMotion(i, 1 / amp_modes[i], x_modes[:, 3*i:3*i+3])
diff --git a/sdpm/utils.py b/sdpm/utils.py
index c20891ab85ab776f82359c13105fe24b4af9a29a..169b9eafa45b24201a90232c1df6e3e6162cd237 100644
--- a/sdpm/utils.py
+++ b/sdpm/utils.py
@@ -118,7 +118,7 @@ def interpolate_modes(fname, surf_nodes, flat=True):
     Parameters
     ----------
     fname : string
-        name of the .csv file (without extension)
+        name of the CSV file containing the modes
     surf_nodes : array
         nodes of surface grid
     flat : bool
@@ -136,7 +136,7 @@ def interpolate_modes(fname, surf_nodes, flat=True):
     # Define number of columns containing coordinates to use
     ndim = 2 if flat else 3
     # Read modal data and remove duplicates
-    modes = np.loadtxt(fname+'.csv', delimiter=',', skiprows=1)
+    modes = np.loadtxt(fname, delimiter=',', skiprows=1)
     _, idx = np.unique(modes[:,:ndim], axis=0, return_index=True)
     modes = modes[idx,:]
     # Interpolate on surface mesh