diff --git a/flow/benchmark/onera.py b/flow/benchmark/onera.py
index b65789de12aa720ff7e0f9fd74bee1dc0034f590..962b830e8d1237cd3515d8319ebe109e42e73fff 100644
--- a/flow/benchmark/onera.py
+++ b/flow/benchmark/onera.py
@@ -16,8 +16,8 @@
 # limitations under the License.
 
 
-# @brief Compute flow around the Onera M6 wing at 3 degrees AOA and Mach 0.84 for benchmark
-# @authors Adrien Crovato
+## Compute flow around the Onera M6 wing at 3 degrees AOA and Mach 0.84 for benchmark
+# Adrien Crovato
 
 import numpy as np
 import flow.default as floD
@@ -88,7 +88,7 @@ def main():
     tms['solver'].start()
     solver = newton(pbl)
     solver.run()
-    solver.save(0, writer)
+    solver.save(writer)
     tms['solver'].stop()
 
     # display timers
diff --git a/flow/cases/coyote.py b/flow/cases/coyote.py
index dbda2d576cf39b90ae5bc78300f06fcca6e21fe8..d5c0b5641c1ddd3b9a3d4f147b8d5ff5df99bd23 100644
--- a/flow/cases/coyote.py
+++ b/flow/cases/coyote.py
@@ -16,7 +16,7 @@
 # limitations under the License.
 
 
-# Coyote - Advanced Pilot Training Aircraft
+## Coyote - Advanced Pilot Training Aircraft
 # Plane designed for the AIAA 2017-2018 aircraft design competition
 # by Raphael Dubois, Thibault Laurent, Bao Long Le Van, Nayan Levaux,
 #    Guillaume Noiset, Axel Piret, Arthur Scheffer, and Vincent Schmitz
diff --git a/flow/cases/n0012.py b/flow/cases/n0012.py
index 9979bdbfe44a8c65a9d5062789147c3f34faeacb..be2da32978ae60696d2ea2a8c6f6103efb309087 100644
--- a/flow/cases/n0012.py
+++ b/flow/cases/n0012.py
@@ -16,8 +16,8 @@
 # limitations under the License.
 
 
-# @brief NACA0012 airfoil configuration file for flow polar script
-# @authors Adrien Crovato
+## NACA 0012 airfoil configuration file for flow polar script
+# Adrien Crovato
 
 def getParam():
     p = {}
diff --git a/flow/cases/n0012_3.py b/flow/cases/n0012_3.py
index e5564d354aaf567f74783e6029c4e6681407d9bb..0505e4362577727eed00876c41a42dc42b008c15 100644
--- a/flow/cases/n0012_3.py
+++ b/flow/cases/n0012_3.py
@@ -16,8 +16,8 @@
 # limitations under the License.
 
 
-# @brief NACA0012 wing configuration file for flow polar script
-# @authors Adrien Crovato
+## NACA 0012 wing configuration file for flow polar script
+# Adrien Crovato
 
 def getParam():
     p = {}
@@ -29,7 +29,7 @@ def getParam():
     p['Dim'] = 3 # Problem dimension
     p['Format'] = 'vtk' # Save format (vtk or gmsh)
     p['Slice'] = [0.01, 0.25, 0.5, 0.75, 0.95] # array of (y) coordinates to perform slice along the span (empty if none)
-    p['TagId'] = 5 # tag number of physical group to be sliced
+    p['TagId'] = 4 # tag number of physical group to be sliced
     # Markers
     p['Fluid'] = 'field' # Name of physical group containing the fluid
     p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of name of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
@@ -69,4 +69,4 @@ def main():
     polar.disp()
 
 if __name__ == "__main__":
-    main()
\ No newline at end of file
+    main()
diff --git a/flow/cases/n64A410.py b/flow/cases/n64A410.py
index b020ff4394111192d7bbd43414df8b03acbe6bd7..219095dccfb8b87592e084813b88468e34aa7941 100644
--- a/flow/cases/n64A410.py
+++ b/flow/cases/n64A410.py
@@ -16,8 +16,8 @@
 # limitations under the License.
 
 
-# @brief NACA64A210 airfoil configuration file for flow trim script
-# @authors Adrien Crovato
+## NACA 64A210 airfoil configuration file for flow trim script
+# Adrien Crovato
 
 def getParam():
     p = {}
diff --git a/flow/cases/rae2822.py b/flow/cases/rae2822.py
index ec5061d157877f752de234875380ac38dd4a7759..e234637ce4f165626a6bc987b21478e1777c9f3d 100644
--- a/flow/cases/rae2822.py
+++ b/flow/cases/rae2822.py
@@ -16,8 +16,8 @@
 # limitations under the License.
 
 
-# @brief RAE2822 airfoil configuration file for flow polar script
-# @authors Adrien Crovato
+## RAE 2822 airfoil configuration file for flow polar script
+# Adrien Crovato
 
 def getParam():
     p = {}
diff --git a/flow/cases/wbht.py b/flow/cases/wbht.py
index 1e25c020977aaa8391b3b12fb55c7897df6c2136..39723b8c6895918ee1fbe42680c09f23b9c8868b 100644
--- a/flow/cases/wbht.py
+++ b/flow/cases/wbht.py
@@ -16,8 +16,8 @@
 # limitations under the License.
 
 
-# @brief Wing-body-tail configuration file for flow trim script
-# @authors Adrien Crovato
+## Wing-body-tail configuration file for flow trim script
+# Adrien Crovato
 #
 # CAUTION
 # Gmsh might not generate the wing/tail surfaces properly. Results might not be accurate 
@@ -94,4 +94,4 @@ def main():
     trim.disp()
 
 if __name__ == "__main__":
-    main()
\ No newline at end of file
+    main()
diff --git a/flow/cases/wht.py b/flow/cases/wht.py
index 19e6de7b7a09af4ddba6d86e9f2cee1316f22eef..3b11cb2600a87fcc71189446d78fda2283c231fb 100644
--- a/flow/cases/wht.py
+++ b/flow/cases/wht.py
@@ -16,8 +16,8 @@
 # limitations under the License.
 
 
-# @brief Wing-tail configuration file for flow trim script
-# @authors Adrien Crovato
+## Wing-tail configuration file for flow trim script
+# Adrien Crovato
 
 def getParam():
     p = {}
@@ -85,4 +85,4 @@ def main():
     trim.disp()
 
 if __name__ == "__main__":
-    main()
\ No newline at end of file
+    main()
diff --git a/flow/default.py b/flow/default.py
index ebc7aeecf54bf2949fbc38fa207d613cee110d8d..7b38f3c3feefe7f3689ce1f678aef1ee5e3ae8f8 100644
--- a/flow/default.py
+++ b/flow/default.py
@@ -17,15 +17,15 @@
 
 
 ## Default initialization for flow tests
-#
 # Adrien Crovato
 
 from fwk.wutils import parseargs
 import flow as flo
 import tbox
 
-## Initialize mesh and mesh writer
 def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
+    '''Initialize mesh and mesh writer
+    '''
     import tbox.gmsh as gmsh
     # load the mesh
     msh = gmsh.MeshLoader(file,__file__).execute(**pars)
@@ -41,8 +41,9 @@ def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
     gmshWriter.save(msh.name)
     return msh, gmshWriter
 
-## Initialize problem
 def problem(msh, dim, alpha, beta, minf, mcrit, sref, lref, xref, yref, zref, sur, fld = 'field', upstr = 'upstream', ff = 'farfield', dnstr = 'downstream', wk = 'wake', te = None, tp = None, vsc = False):
+    '''Initialize problem
+    '''
     pbl = flo.Problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref)
     if minf == 0:
         pbl.set(flo.Medium(msh, fld, flo.F0ElRhoL(), flo.F0ElMachL(), flo.F0ElCpL(), flo.F0PsPhiInf(dim, alpha, beta)))
@@ -74,8 +75,9 @@ def problem(msh, dim, alpha, beta, minf, mcrit, sref, lref, xref, yref, zref, su
         blww = None
     return pbl, dirichlet, wake, bnd, [blw, blww]
 
-## Initialize Picard solver
 def picard(pbl):
+    '''Initialize Picard solver
+    '''
     args = parseargs()
     solver = flo.Picard(pbl, tbox.Gmres())
     solver.nthreads = args.k
@@ -86,8 +88,9 @@ def picard(pbl):
     solver.relax = 0.7
     return solver
 
-## Initialize Newton solver
 def newton(pbl):
+    '''Initialize Newton solver
+    '''
     from tbox.solvers import LinearSolver
     args = parseargs()
     solver = flo.Newton(pbl, LinearSolver().pardiso())
@@ -101,8 +104,9 @@ def newton(pbl):
     solver.avThrsh = 1e-2
     return solver
 
-## Initialize mesh morpher
 def meshMorpher(msh, dim, mov, fxd = ['upstream', 'farfield', 'downstream'], fld = 'field', wk = 'wake', sym = 'symmetry'):
+    '''Initialize mesh morpher
+    '''
     args = parseargs()
     mshDef = tbox.MshDeform(msh, dim)
     mshDef.nthreads = args.k
@@ -115,8 +119,9 @@ def meshMorpher(msh, dim, mov, fxd = ['upstream', 'farfield', 'downstream'], fld
     mshDef.initialize()
     return mshDef
 
-## Initialize viewer
 def initViewer(problem):
+    '''Initialize viewer
+    '''
     args = parseargs()
     if not args.nogui:
         from tbox.viewer import GUI
diff --git a/flow/scripts/config.py b/flow/scripts/config.py
index 09ab1757388ef34f92750409773f2960cd713dfa..89306a1303b90adec7b79f503854da520fdbd374 100644
--- a/flow/scripts/config.py
+++ b/flow/scripts/config.py
@@ -16,8 +16,8 @@
 # limitations under the License.
 
 
-# @brief Base class for configuration
-# @authors Adrien Crovato
+## Base class for configuration
+# Adrien Crovato
 
 import flow as f
 import tbox
diff --git a/flow/scripts/polar.py b/flow/scripts/polar.py
index 8adb98b55279245068185a5c83ab545ce76def01..3d1ade85974322c6b09e31f7337bdd3c87eb3367 100644
--- a/flow/scripts/polar.py
+++ b/flow/scripts/polar.py
@@ -16,9 +16,10 @@
 # limitations under the License.
 
 
-# @brief Polar
-#  Compute lift and polar curves
-# @authors Adrien Crovato
+## Polar
+# Adrien Crovato
+#
+# Compute lift and polar curves
 
 import math
 import flow.utils as fU
@@ -66,15 +67,14 @@ class Polar(Config):
             self.tms["solver"].start()
             self.solver.run()
             self.tms["solver"].stop()
-            self.solver.save(ac, self.mshWriter)
+            self.solver.save(self.mshWriter, ac)
             print('\n')
             # extract Cp
-            if self.alphas[0] == self.alphas[-1]:
-                if self.dim == 2:
-                    Cp = fU.extract(self.bnd.groups[0].tag.elems, self.solver.Cp)
-                    tU.write(Cp, "Cp_airfoil.dat", "%1.5e", ", ", "x, y, z, Cp", "")
-                elif self.dim == 3 and self.format == 'vtk' and self.slice:
-                    fU.writeSlices(self.msh.name, self.slice, self.tag)
+            if self.dim == 2:
+                Cp = fU.extract(self.bnd.groups[0].tag.elems, self.solver.Cp)
+                tU.write(Cp, 'Cp_airfoil_'+str(ac)+'.dat', '%1.5e', ', ', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, z, Cp', '')
+            elif self.dim == 3 and self.format == 'vtk' and self.slice:
+                fU.writeSlices(self.msh.name, self.slice, self.tag, ac)
             # extract force coefficients
             self.Cls.append(self.solver.Cl)
             self.Cds.append(self.solver.Cd)
diff --git a/flow/scripts/trim.py b/flow/scripts/trim.py
index c9f1d954fcfdf989c48afcb14dec2c4762938f43..71e38240d49b888e0701b5067613d8d72d761fab 100644
--- a/flow/scripts/trim.py
+++ b/flow/scripts/trim.py
@@ -16,9 +16,10 @@
 # limitations under the License.
 
 
-# @brief Trim analysis
-#  Find the angle of attack to match a specified lift coefficient
-# @authors Adrien Crovato
+## Trim analysis
+# Adrien Crovato
+#
+# Find the angle of attack to match a specified lift coefficient
 
 import math
 import flow.utils as fU
@@ -79,7 +80,7 @@ class Trim(Config):
                 it += 1
             
         # save results
-        self.solver.save(0, self.mshWriter)
+        self.solver.save(self.mshWriter)
         print('\n')
         # extract Cp
         if self.dim == 2:
diff --git a/flow/src/wAdjoint.cpp b/flow/src/wAdjoint.cpp
index b02ef4ca864254b2f1f64178d75459955f403653..c71122416efdf10fb94d5e79e31dd7f298158d51 100644
--- a/flow/src/wAdjoint.cpp
+++ b/flow/src/wAdjoint.cpp
@@ -413,7 +413,7 @@ void Adjoint::buildWake(WakeElement *&we, Eigen::MatrixXd &Kupup, Eigen::MatrixX
  * @brief Write the results
  * @authors Adrien Crovato
  */
-void Adjoint::save(int n, MshExport *mshWriter)
+void Adjoint::save(MshExport *mshWriter, int n)
 {
     // Write files
     std::cout << "Saving files... " << std::endl;
diff --git a/flow/src/wAdjoint.h b/flow/src/wAdjoint.h
index 54cd797025cde497b7f90557348c1950d9f312bc..e0943ba6ca0864445dc23a330e2cef9081bca4df 100644
--- a/flow/src/wAdjoint.h
+++ b/flow/src/wAdjoint.h
@@ -49,7 +49,7 @@ public:
     virtual ~Adjoint() { std::cout << "~Adjoint()\n"; }
 
     void run();
-    void save(int n, tbox::MshExport *mshWriter);
+    void save(tbox::MshExport *mshWriter, int n = 0);
 
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;
diff --git a/flow/src/wSolver.cpp b/flow/src/wSolver.cpp
index 686ccad51845ce14c5f1706b4e24a47788e0a8a0..0fda55d6ea4cb30d00029be4ac97d2bb0f67db89 100644
--- a/flow/src/wSolver.cpp
+++ b/flow/src/wSolver.cpp
@@ -57,7 +57,7 @@ Solver::Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<LinearSolver> _lin
     std::cout << "**                          | __| | | /  \\ \\ \\_/ /                           **" << std::endl;
     std::cout << "**                          |_|   |_| \\__/  \\/ \\/                            **" << std::endl;
     std::cout << "*******************************************************************************" << std::endl;
-    std::cout << "** Hi! My name is Flow v1.6-2004                                             **" << std::endl;
+    std::cout << "** Hi! My name is Flow v1.7-2011                                             **" << std::endl;
     std::cout << "** Adrien Crovato & Romain Boman                                             **" << std::endl;
     std::cout << "** ULiege 2018-2020                                                          **" << std::endl;
     std::cout << "*******************************************************************************" << std::endl
@@ -126,7 +126,7 @@ bool Solver::run()
  * @brief Write the results
  * @authors Adrien Crovato
  */
-void Solver::save(int n, MshExport *mshWriter)
+void Solver::save(MshExport *mshWriter, int n)
 {
     // Write files
     std::cout << "Saving files... " << std::endl;
diff --git a/flow/src/wSolver.h b/flow/src/wSolver.h
index ba3cf3a75b02e77493898520affb125f59fc4c39..0a176700041d77797f7ebb949d98d21432bf9be6 100644
--- a/flow/src/wSolver.h
+++ b/flow/src/wSolver.h
@@ -67,7 +67,7 @@ public:
     virtual ~Solver();
 
     virtual bool run();
-    void save(int n, tbox::MshExport *mshWriter);
+    void save(tbox::MshExport *mshWriter, int n = 0);
 
 protected:
     void computeFlow();
diff --git a/flow/tests/adjoint.py b/flow/tests/adjoint.py
index 2241d156e3e40e12875f140f7d91b3a4e7e1da3f..58d0fe779feab267efd213d024b398f0668848ee 100644
--- a/flow/tests/adjoint.py
+++ b/flow/tests/adjoint.py
@@ -16,8 +16,9 @@
 # limitations under the License.
 
 
-# @brief Compute adjoint solution of lifting (linear or nonlinear) potential flow around a naca0012
-# @authors Adrien Crovato
+## Compute adjoint solution of lifting (linear or nonlinear) potential flow around a NACA 0012
+# Adrien Crovato
+#
 # Test the adjoint solver
 #
 # CAUTION
@@ -63,14 +64,14 @@ def main():
     tms['solver'].start()
     solver = floD.newton(pbl)
     solver.run()
-    solver.save(0, gmshWriter)
+    solver.save(gmshWriter)
     tms['solver'].stop()
 
     # solve adjoint problem
     tms['adjoint'].start()
     adjoint = flo.Adjoint(solver)
     adjoint.run()
-    adjoint.save(0, gmshWriter)
+    adjoint.save(gmshWriter)
     tms['adjoint'].stop()
 
     # extract Cp
diff --git a/flow/tests/bli.py b/flow/tests/bli.py
index c6e8cbc220fad581f4e3fa50f8fec8f19744f8a2..4c936f19f4e9114fa2494eaed0cdaa52f69e00cc 100644
--- a/flow/tests/bli.py
+++ b/flow/tests/bli.py
@@ -17,8 +17,8 @@
 
 
 ## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
-#
 # Amaury Bilocq
+#
 # Test the viscous-inviscid interaction scheme
 # Reference to the master's thesis: http://hdl.handle.net/2268/252195
 # Reference test cases with Naca0012 (different from master's thesis): 
@@ -132,4 +132,4 @@ def main():
     print('')
 
 if __name__ == "__main__":
-    main()
\ No newline at end of file
+    main()
diff --git a/flow/tests/cylinder.py b/flow/tests/cylinder.py
index 44ddfd1d4c98292aced25f2df6c6d09ec0798ec0..c5d09d07ad74567d6a36ff76e563943dc9a07e71 100644
--- a/flow/tests/cylinder.py
+++ b/flow/tests/cylinder.py
@@ -16,8 +16,9 @@
 # limitations under the License.
 
 
-# @brief Compute the (non)linear flow around cylinder, with various B.C.
-# @authors Adrien Crovato
+## Compute the (non)linear flow around cylinder, with various B.C.
+# Adrien Crovato
+#
 # Test solver convergence and boundary conditions
 #
 # CAUTION
@@ -67,12 +68,12 @@ def main():
     tms['picard'].start()
     solver0 = floD.picard(pbl)
     cnvrgd0 = solver0.run()
-    solver0.save(0, gmshWriter)
+    solver0.save(gmshWriter)
     tms['picard'].stop()
     tms['newton'].start()
     solver1 = floD.newton(pbl)
     cnvrgd1 = solver1.run()
-    solver1.save(0, gmshWriter)
+    solver1.save(gmshWriter)
     tms['newton'].stop()
 
     # compute mean error in Boundary, Neumann boundary condition (works only for Line2 elements)
diff --git a/flow/tests/cylinder2D5.py b/flow/tests/cylinder2D5.py
index 76ef1e207b8b9bb87b3a0c6cb557870816306365..087f5dfb1b1896ccc7adf9d2b91abcef83b07f8f 100644
--- a/flow/tests/cylinder2D5.py
+++ b/flow/tests/cylinder2D5.py
@@ -16,8 +16,9 @@
 # limitations under the License.
 
 
-# @brief Compute the (non)linear flow around extruded 2D cylinder, with various B.C.
-# @authors Adrien Crovato
+## Compute the (non)linear flow around extruded 2D cylinder, with various B.C.
+# Adrien Crovato
+#
 # Test solver convergence and boundary conditions
 #
 # CAUTION
@@ -69,7 +70,7 @@ def main():
     tms['solver'].start()
     solver = floD.newton(pbl)
     cnvrgd = solver.run()
-    solver.save(0, gmshWriter)
+    solver.save(gmshWriter)
     tms['solver'].stop()
 
     # compute mean error in Boundary, Neumann boundary condition (works only for Tri3 elements)
diff --git a/flow/tests/cylinder3.py b/flow/tests/cylinder3.py
index 38462d4e4c3ea933cce5526282bfb3ba46453d6f..e618daf76f457e3957bb1e40341bb607dbb25a94 100644
--- a/flow/tests/cylinder3.py
+++ b/flow/tests/cylinder3.py
@@ -16,8 +16,9 @@
 # limitations under the License.
 
 
-# @brief Compute the (non)linear flow around 3D cylinder, with various B.C.
-# @authors Adrien Crovato
+## Compute the (non)linear flow around 3D cylinder, with various B.C.
+# Adrien Crovato
+#
 # Test solver convergence and boundary conditions
 #
 # CAUTION
@@ -72,7 +73,7 @@ def main():
     tms['solver'].start()
     solver = floD.newton(pbl)
     cnvrgd = solver.run()
-    solver.save(0, gmshWriter)
+    solver.save(gmshWriter)
     tms['solver'].stop()
 
     # compute mean error in Boundary, Neumann boundary condition (works only for Tri3 elements)
diff --git a/flow/tests/lift.py b/flow/tests/lift.py
index ca223a02c0ad918c1c42f4006dfb12b43bc1313b..0c4e5dc179c17396eff02d35793dadf30bad4f17 100644
--- a/flow/tests/lift.py
+++ b/flow/tests/lift.py
@@ -16,8 +16,9 @@
 # limitations under the License.
 
 
-# @brief Compute lifting (linear or nonlinear) potential flow around a naca0012
-# @authors Adrien Crovato
+## Compute lifting (linear or nonlinear) potential flow around a NACA 0012
+# Adrien Crovato
+#
 # Test the nonlinear shock-capturing capability and the Kutta condition
 #
 # CAUTION
@@ -61,7 +62,7 @@ def main():
     tms['solver'].start()
     solver = floD.newton(pbl)
     solver.run()
-    solver.save(0, gmshWriter)
+    solver.save(gmshWriter)
     tms['solver'].stop()
 
     # extract Cp
diff --git a/flow/tests/lift3.py b/flow/tests/lift3.py
index bbbc9c44bd503cbe41c691a6235b29f603cdfe6c..5d754533ccbd64782851810ea3f70963f3727771 100644
--- a/flow/tests/lift3.py
+++ b/flow/tests/lift3.py
@@ -16,8 +16,9 @@
 # limitations under the License.
 
 
-# @brief Compute lifting (linear or nonlinear) potential flow around a naca0012 rectangular wing
-# @authors Adrien Crovato
+## Compute lifting (linear or nonlinear) potential flow around a NACA 0012 rectangular wing
+# Adrien Crovato
+#
 # Test the nonlinear shock-capturing capability and the Kutta condition
 #
 # CAUTION
@@ -69,7 +70,7 @@ def main():
     tms['solver'].start()
     solver = floD.newton(pbl)
     solver.run()
-    solver.save(0, gmshWriter)
+    solver.save(gmshWriter)
     tms['solver'].stop()
 
     # display timers
diff --git a/flow/tests/meshDef.py b/flow/tests/meshDef.py
index c92090da69252cc661504299c2871e0dc796c2bc..a5d37315f9f66eccb295619410601644ba3fdc40 100644
--- a/flow/tests/meshDef.py
+++ b/flow/tests/meshDef.py
@@ -16,8 +16,9 @@
 # limitations under the License.
 
 
-# @brief Compute the flow on a deformed mesh
-# @authors Adrien Crovato
+## Compute the flow around a NACA 0012 on a deformed mesh
+# Adrien Crovato
+#
 # Test the mesh deformation process
 #
 # CAUTION
@@ -99,7 +100,7 @@ def main():
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
     tms['solver1'].start()
     solver.run()
-    solver.save(0, gmshWriter)
+    solver.save(gmshWriter)
     tms['solver1'].stop()
 
     # compute reference solution
@@ -107,7 +108,7 @@ def main():
     tms['solver_ref'].start()
     solver_ref = floD.newton(pbl_ref)
     solver_ref.run()
-    solver_ref.save(0, gmshWriter_ref)
+    solver_ref.save(gmshWriter_ref)
     tms['solver_ref'].stop()
 
     # extract Cp
diff --git a/flow/tests/meshDef3.py b/flow/tests/meshDef3.py
index 1997450e37d4b7b59648fc453e94cc85aca3628b..cbfab2768d9a5339b8c4b1ade01a1558842e82a3 100644
--- a/flow/tests/meshDef3.py
+++ b/flow/tests/meshDef3.py
@@ -16,8 +16,9 @@
 # limitations under the License.
 
 
-# @brief Compute the flow on a deformed mesh
-# @authors Adrien Crovato
+## Compute the flow around a NACA 0012 wing on a deformed mesh
+# Adrien Crovato
+#
 # Test the mesh deformation process
 #
 # CAUTION
@@ -102,7 +103,7 @@ def main():
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
     tms['solver1'].start()
     solver.run()
-    solver.save(0, gmshWriter)
+    solver.save(gmshWriter)
     tms['solver1'].stop()
 
     # display results
diff --git a/flow/tests/nonlift.py b/flow/tests/nonlift.py
index d9fd64eba9c9358927a6e0a0ad18325f30508aaf..db2282732c635c4c26074d4f8af071c7d1388e71 100644
--- a/flow/tests/nonlift.py
+++ b/flow/tests/nonlift.py
@@ -16,8 +16,9 @@
 # limitations under the License.
 
 
-# @brief Compute non-lifting (linear or nonlinear) potential flow around a naca0012
-# @authors Adrien Crovato
+## Compute non-lifting (linear or nonlinear) potential flow around a NACA 0012
+# Adrien Crovato
+#
 # Test the nonlinear shock-capturing capability
 #
 # CAUTION
@@ -62,7 +63,7 @@ def main():
     tms['solver'].start()
     solver = floD.newton(pbl)
     solver.run()
-    solver.save(0, gmshWriter)
+    solver.save(gmshWriter)
     tms['solver'].stop()
 
     # extract Cp
diff --git a/flow/utils.py b/flow/utils.py
index f0e0a0996817e2bf6677ed695a6fff9d3331455b..359c7cc95739563f4e99bfb998e25ef96977c9ae 100644
--- a/flow/utils.py
+++ b/flow/utils.py
@@ -16,17 +16,14 @@
 # limitations under the License.
 
 
-"""@package docstring
-Utilities for flow
-Adrien Crovato
-"""
+## Utilities for flow
+# Adrien Crovato
 
 def computeAeroCoef(_x, _y, _Cp, _alpha):
-    """Compute 2D aerodynamic coefficients
+    '''Compute 2D aerodynamic coefficients
     The coefficients are computed from Cp averaged at nodes, which might leads to innacurate results
     For accurate results, use coefficients from solver or boundary objects
-    Adrien Crovato
-    """
+    '''
     import math
     i = 0
     Cy = 0
@@ -44,9 +41,8 @@ def computeAeroCoef(_x, _y, _Cp, _alpha):
     return (Cl, Cd, Cm)
 
 def convex_sort(vNodes, vData):
-    """Sort point cloud data forming a convex hull in counterclowise order, from TE
-    Adrien Crovato
-    """
+    '''Sort point cloud data forming a convex hull in counterclowise order, from TE
+    '''
     import numpy as np
     # store map in arrays
     data = np.zeros((vNodes.size(),4))
@@ -76,9 +72,8 @@ def convex_sort(vNodes, vData):
     return data
 
 def extract(vElems, vData):
-    """Extract data at first node of elements
-    Adrien Crovato
-    """
+    '''Extract data at first node of elements
+    '''
     import numpy as np
     data = np.zeros((vElems.size()+1,4))
     i = 0
@@ -90,10 +85,9 @@ def extract(vElems, vData):
         i += 1
     return data
 
-def writeSlices(mshName, ys, wId):
-    """Write slice data for each ys coordinate along the span (only works with VTK)
-    Adrien Crovato
-    """
+def writeSlices(mshName, ys, wId, n = 0):
+    '''Write slice data for each ys coordinate along the span (only works with VTK)
+    '''
     import numpy as np
     try:
         import vtk
@@ -103,7 +97,10 @@ def writeSlices(mshName, ys, wId):
     import tboxVtk.reader as vtkR
     import tboxVtk.cutter as vtkC
     reader = vtkR.Reader()
-    reader.open(mshName)
+    if n == 0:
+        reader.open(mshName)
+    else:
+        reader.open(mshName+'_'+str(n))
     cutter = vtkC.Cutter(reader.reader.GetOutputPort())
     for i in range(0, len(ys)):
         pts, elems, vals = cutter.extract(cutter.cut(wId, [0., ys[i], 0.], [0., 1., 0.]), 2, ['Cp'])
@@ -113,4 +110,7 @@ def writeSlices(mshName, ys, wId):
         data = np.hstack((data, vals['Cp']))
         tU.sort(elems, data)
         data = np.vstack((data, data[0,:]))
-        tU.write(data, 'slice_'+str(i)+'.dat', "%1.5e", ", ", "x, y, z, x/c, Cp", "")
+        if n == 0:
+            tU.write(data, 'slice_'+str(i)+'.dat', '%1.5e', ', ', 'x, y, z, x/c, Cp', '')
+        else:
+            tU.write(data, 'slice_'+str(n)+'_'+str(i)+'.dat', '%1.5e', ', ', 'x, y, z, x/c, Cp', '')
diff --git a/flow/validation/agard.py b/flow/validation/agard.py
index 2a8d488adc70ac98e19f85655108d4c54732a32f..6fd48c4617ef8687075f0384ec5c456b71526e9c 100644
--- a/flow/validation/agard.py
+++ b/flow/validation/agard.py
@@ -16,8 +16,8 @@
 # limitations under the License.
 
 
-# @brief Compute flow around the Agard445 wing at 1 degrees AOA and Mach 0.80
-# @authors Adrien Crovato
+## Compute flow around the Agard 445 wing at 1 degrees AOA and Mach 0.80
+# Adrien Crovato
 
 import numpy as np
 import flow.utils as floU
@@ -78,7 +78,7 @@ def main():
     tms['solver'].start()
     solver = floD.newton(pbl)
     solver.run()
-    solver.save(0, writer)
+    solver.save(writer)
     tms['solver'].stop()
 
     # post process
diff --git a/flow/validation/onera.py b/flow/validation/onera.py
index e53482cbeaccce4bec392628600dadcf45c3db64..afedc0941657e463423dd5b09371a69465b27617 100644
--- a/flow/validation/onera.py
+++ b/flow/validation/onera.py
@@ -16,8 +16,8 @@
 # limitations under the License.
 
 
-# @brief Compute flow around the Onera M6 wing at 3 degrees AOA and Mach 0.84
-# @authors Adrien Crovato
+## Compute flow around the Onera M6 wing at 3 degrees AOA and Mach 0.84
+# Adrien Crovato
 
 import numpy as np
 import flow.utils as floU
@@ -76,7 +76,7 @@ def main():
     tms['solver'].start()
     solver = floD.newton(pbl)
     solver.run()
-    solver.save(0, writer)
+    solver.save(writer)
     tms['solver'].stop()
 
     # post process
diff --git a/flow/viscous/airfoil.py b/flow/viscous/airfoil.py
index 15fc43c7374a45c43548d3dc7478e28cb9c6e96c..e76a60931a4157457273261f3ce2f2ead4a534ef 100644
--- a/flow/viscous/airfoil.py
+++ b/flow/viscous/airfoil.py
@@ -17,7 +17,6 @@
 
 
 ## Airfoil around which the boundary layer is computed
-#
 # Amaury Bilocq
 
 from flow.viscous.boundary import Boundary
@@ -129,4 +128,4 @@ class Airfoil(Boundary):
         uData = data[0:np.argmax(data[:,0])+1]
         lData = data[np.argmax(data[:,0])+1:None]
         lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
-        return connectListNodes, connectListElems,[uData, lData]
\ No newline at end of file
+        return connectListNodes, connectListElems,[uData, lData]
diff --git a/flow/viscous/boundary.py b/flow/viscous/boundary.py
index 3a8e862bb038441ce3b6d8398dc22a13f4247a59..773cf7f9d17fef94a120d9f7ea463e93363d09d1 100644
--- a/flow/viscous/boundary.py
+++ b/flow/viscous/boundary.py
@@ -17,7 +17,6 @@
 
 
 ## Base class representing a physical boundary
-#
 # Amaury Bilocq
 
 import numpy as np
diff --git a/flow/viscous/coupler.py b/flow/viscous/coupler.py
index 4a672c20c63881acd8b26dfafe98ca1c1348e3e4..605bc4beb072b076ee1c688d69401c9cbb465240 100644
--- a/flow/viscous/coupler.py
+++ b/flow/viscous/coupler.py
@@ -17,7 +17,6 @@
 
 
 ## Viscous-inviscid coupler (quasi-simultaneous coupling)
-#
 # Amaury Bilocq
 
 import numpy as np
@@ -80,6 +79,6 @@ class Coupler:
             it += 1
             self.vsolver.it += 1     
         # save results
-        self.isolver.save(0, self.writer)
+        self.isolver.save(self.writer)
         self.vsolver.writeFile()
-        print('\n')
\ No newline at end of file
+        print('\n')
diff --git a/flow/viscous/newton.py b/flow/viscous/newton.py
index bed39ec64ca0a43df345cce476da8c4c5e8da6c2..6c068686bc56fbd992b028a4f1aec8e6b13a9b57 100644
--- a/flow/viscous/newton.py
+++ b/flow/viscous/newton.py
@@ -17,7 +17,6 @@
 
 
 ## Newton raphson method coupled with linear solver
-#
 # Amaury Bilocq
 
 import numpy as np
diff --git a/flow/viscous/solver.py b/flow/viscous/solver.py
index 54f97a45fd379b73e28d80e2a4112c54e86e8457..da3a9914b9b6d96dc1f663c35f95e753eb575ca8 100644
--- a/flow/viscous/solver.py
+++ b/flow/viscous/solver.py
@@ -17,7 +17,6 @@
 
 
 ## Integral boundary layer equations viscous solver
-#
 # Amaury Bilocq
 #
 # TODO list
@@ -475,4 +474,4 @@ class Solver:
         # Sort the following results in reference frame
         group.deltaStar = group.deltaStar[group.connectListNodes.argsort()] 
         group.xx = group.xx[group.connectListNodes.argsort()]
-        group.u = blwVel[group.connectListElems.argsort()]
\ No newline at end of file
+        group.u = blwVel[group.connectListElems.argsort()]
diff --git a/flow/viscous/wake.py b/flow/viscous/wake.py
index 0f3cd5e5c3cef24ebca096ab6d7f8c439890245f..c925c2eacdfa8303dc45172c07e428d67987b06d 100644
--- a/flow/viscous/wake.py
+++ b/flow/viscous/wake.py
@@ -17,7 +17,6 @@
 
 
 ## Wake behind airfoil (around which the boundary layer is computed)
-#
 # Amaury Bilocq
 
 from flow.viscous.boundary import Boundary
@@ -76,4 +75,4 @@ class Wake(Boundary):
         data[:,9] = data[connectListNodes,9]
         # Separated upper and lower part 
         data = np.delete(data,0,1)
-        return connectListNodes, connectListElems, data
\ No newline at end of file
+        return connectListNodes, connectListElems, data