diff --git a/dart/api/internal/polar.py b/dart/api/internal/polar.py
index 198b57ea6c265e7036defadc64fa372350c2c7ab..d1f49acbeeb978d40ca85aa7f03b70c0f095dafa 100644
--- a/dart/api/internal/polar.py
+++ b/dart/api/internal/polar.py
@@ -98,7 +98,7 @@ class Polar:
             # extract Cp
             if self.dim == 2:
                 Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
-                tU.write(Cp, f'Cp_{self.msh.name}_airfoil{acs}.dat', '%1.5e', ', ', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, z, Cp', '')
+                tU.write(Cp, f'Cp_{self.msh.name}_airfoil{acs}.dat', '%1.4e', ',', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, cp', '')
             elif self.dim == 3 and self.format == 'vtk' and self.slice:
                 dU.write_slices(self.msh.name, self.slice, self.tag, acs)
             # extract force coefficients
diff --git a/dart/api/internal/trim.py b/dart/api/internal/trim.py
index 31c1028047cad01c86212c0451384db2ad901ea6..9e0c263a6d0ac2845ebdc3a223659e4f6e0563a7 100644
--- a/dart/api/internal/trim.py
+++ b/dart/api/internal/trim.py
@@ -110,7 +110,7 @@ class Trim:
         # extract Cp
         if self.dim == 2:
             Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
-            tU.write(Cp, f'Cp_{self.msh.name}_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+            tU.write(Cp, f'Cp_{self.msh.name}_airfoil.dat', '%1.4e', ',', 'alpha = '+str(self.alpha*180/math.pi)+' deg\nx, y, cp', '')
         elif self.dim == 3 and self.format == 'vtk' and self.slice:
             dU.write_slices(self.msh.name, self.slice, self.tag)
 
diff --git a/dart/tests/adjoint.py b/dart/tests/adjoint.py
index 0834f241870e465993754a7d35a3e07e5efe0ac7..c2795d617e40402443b1bc53b3d6d7c9d54ae442 100644
--- a/dart/tests/adjoint.py
+++ b/dart/tests/adjoint.py
@@ -94,9 +94,9 @@ def main():
 
     # extract Cp and drag sensitivities
     Cp = floU.extract(bnd.groups[0].tag.elems, solver.Cp)
-    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.4e', ',', 'x, y, cp', '')
     dCd = floU.extract(bnd.groups[0].tag.elems, adjoint.dCdMsh)
-    tboxU.write(dCd, 'dCd_airfoil.dat', '%1.5e', ', ', 'x, y, z, dCd_dx, dCd_dy, dCd_dz', '')
+    tboxU.write(dCd, 'dCd_airfoil.dat', '%1.4e', ',', 'x, y, dCd_dx, dCd_dy, dCd_dz', '')
     # display results
     print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
     print('       M    alpha       Cl       Cd   dCl_da   dCd_da')
diff --git a/dart/tests/lift.py b/dart/tests/lift.py
index a541a8368e1d8d04883bd9fa0cfa5c86661bc27b..981eb9f95cf16f2efdae8c2c04fcd296495a08d2 100644
--- a/dart/tests/lift.py
+++ b/dart/tests/lift.py
@@ -66,7 +66,7 @@ def main():
 
     # extract Cp
     Cp = floU.extract(bnd.groups[0].tag.elems, solver.Cp)
-    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.4e', ',', 'x, y, cp', '')
     # display results
     print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
     print('       M    alpha       Cl       Cd       Cm')
@@ -80,22 +80,22 @@ def main():
 
     # visualize solution and plot results
     floD.initViewer(pbl)
-    floD.plot(Cp[:,0], Cp[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
+    floD.plot(Cp[:,0], Cp[:,-1], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
 
     # check results
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
     if M_inf == 0 and alpha == 3*math.pi/180:
-        tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.1, 1.5e-1))
+        tests.add(CTest('min(Cp)', min(Cp[:,-1]), -1.1, 1.5e-1))
         tests.add(CTest('Cl', solver.Cl, 0.36, 5e-2))
         tests.add(CTest('Cm', solver.Cm, 0.0, 1e-2))
     elif M_inf == 0.6 and alpha == 2*math.pi/180:
-        tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.03, 1e-1))
+        tests.add(CTest('min(Cp)', min(Cp[:,-1]), -1.03, 1e-1))
         tests.add(CTest('Cl', solver.Cl, 0.315, 5e-2))
         tests.add(CTest('Cm', solver.Cm, 0.0, 1e-2))
     elif M_inf == 0.7 and alpha == 2*math.pi/180:
         tests.add(CTest('iteration count', solver.nIt, 12, 3, forceabs=True))
-        tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.28, 5e-2))
+        tests.add(CTest('min(Cp)', min(Cp[:,-1]), -1.28, 5e-2))
         tests.add(CTest('Cl', solver.Cl, 0.391, 5e-2))
         tests.add(CTest('Cd', solver.Cd, 0.0013, 5e-4, forceabs=True))
         tests.add(CTest('Cm', solver.Cm, 0.0, 1e-2))
diff --git a/dart/tests/morpher.py b/dart/tests/morpher.py
index ff109a79c72bc2e30467e1aa3495e416320fa059..62400804b6599d1b6de077c23c340f4f0e4cee49 100644
--- a/dart/tests/morpher.py
+++ b/dart/tests/morpher.py
@@ -112,9 +112,9 @@ def main():
 
     # extract Cp
     Cp = floU.extract(bnd.groups[0].tag.elems, solver.Cp)
-    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.4e', ',', 'x, y, cp', '')
     Cp_ref = floU.extract(bnd_ref.groups[0].tag.elems, solver_ref.Cp)
-    tboxU.write(Cp, 'Cp_airfoil_ref.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    tboxU.write(Cp, 'Cp_airfoil_ref.dat', '%1.4e', ',', 'x, y, cp', '')
     # display results
     print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
     print('  case        M    alpha       Cl       Cd       Cm')
@@ -128,15 +128,15 @@ def main():
     print(tms)
 
     # visualize solution and plot results
-    floD.plot(Cp[:,0], Cp[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
-    floD.plot(Cp_ref[:,0], Cp_ref[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver_ref.Cl, solver_ref.Cd, solver_ref.Cm, 4), 'invert': True})
+    floD.plot(Cp[:,0], Cp[:,-1], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
+    floD.plot(Cp_ref[:,0], Cp_ref[:,-1], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver_ref.Cl, solver_ref.Cd, solver_ref.Cm, 4), 'invert': True})
 
     # check results
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
     tests.add(CTest('solver: iteration count', solver.nIt, 5, 3, forceabs=True)) # todo: tolerance should be reset to 1 as soon as gmsh4 is used
     tests.add(CTest('solver_ref: iteration count', solver_ref.nIt, 5, 1, forceabs=True))
-    tests.add(CTest('min(Cp)', min(Cp[:,3]), min(Cp_ref[:,3]), 5e-2))
+    tests.add(CTest('min(Cp)', min(Cp[:,-1]), min(Cp_ref[:,-1]), 5e-2))
     tests.add(CTest('Cl', solver.Cl, solver_ref.Cl, 5e-2))
     tests.add(CTest('Cm', solver.Cm, solver_ref.Cm, 5e-2))
     tests.run()
diff --git a/dart/tests/nonlift.py b/dart/tests/nonlift.py
index 17effbc1737c4a60b0d9977f257ef9e94ac31db8..cbde07046bda95a03dacf3469b403386e8ded256 100644
--- a/dart/tests/nonlift.py
+++ b/dart/tests/nonlift.py
@@ -66,7 +66,7 @@ def main():
 
     # extract Cp
     Cp = floU.extract(bnd.groups[0].tag.elems, solver.Cp)
-    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.4e', ',', 'x, y, cp', '')
     # display results
     print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
     print('       M    alpha       Cl       Cd       Cm')
@@ -80,18 +80,18 @@ def main():
 
     # visualize solution and plot results
     floD.initViewer(pbl)
-    floD.plot(Cp[:,0], Cp[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
+    floD.plot(Cp[:,0], Cp[:,-1], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
 
     # check results
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
     if M_inf == 0:
-        tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.405, 5e-2))
+        tests.add(CTest('min(Cp)', min(Cp[:,-1]), -0.405, 5e-2))
     elif M_inf == 0.7:
-        tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.63, 5e-2))
+        tests.add(CTest('min(Cp)', min(Cp[:,-1]), -0.63, 5e-2))
     elif M_inf == 0.8:
         tests.add(CTest('iteration count', solver.nIt, 13, 3, forceabs=True))
-        tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.89, 5e-2))
+        tests.add(CTest('min(Cp)', min(Cp[:,-1]), -0.89, 5e-2))
         tests.add(CTest('Cd', solver.Cd, 0.0058, 5e-4, forceabs=True))
     else:
         raise Exception('Test not defined for this flow')
diff --git a/dart/utils.py b/dart/utils.py
index 6d846bfa5b8aa773a1f3e1555602a0abec7c37a5..3bf610cb16660d860d3a9b24f6ac470e3d8bc928 100644
--- a/dart/utils.py
+++ b/dart/utils.py
@@ -97,17 +97,16 @@ def extract(vElems, vData):
     else:
         raise RuntimeError('dart.utils.extract: unrecognized format for vData!')
     # store data (not efficient, but OK since only meant for small 2D cases)
-    data = np.zeros((vElems.size()+1,3+size))
+    data = np.zeros((vElems.size()+1,2+size))
     i = 0
     while i < vElems.size()+1:
         data[i,0] = vElems[i%vElems.size()].nodes[0].pos[0]
         data[i,1] = vElems[i%vElems.size()].nodes[0].pos[1]
-        data[i,2] = vElems[i%vElems.size()].nodes[0].pos[2]
         if size == 1:
-            data[i,3] = vData[vElems[i%vElems.size()].nodes[0].row]
+            data[i,2] = vData[vElems[i%vElems.size()].nodes[0].row]
         else:
             for j in range(size):
-                data[i,3+j] = vData[vElems[i%vElems.size()].nodes[0].row][j]
+                data[i,2+j] = vData[vElems[i%vElems.size()].nodes[0].row][j]
         i += 1
     return data
 
@@ -138,10 +137,14 @@ def write_slices(mshName, ys, wId, sfx=''):
     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'])
-        x_c = np.zeros((pts.shape[0],1))
-        x_c[:,0] = (pts[:,0] - min(pts[:,0])) / (max(pts[:,0]) - min(pts[:,0]))
-        data = np.hstack((pts, x_c))
-        data = np.hstack((data, vals['Cp']))
+        x_c = np.zeros((pts.shape[0], 2))
+        ile = np.argmin(pts[:,0]) # LE index
+        crd = max(pts[:,0]) - min(pts[:,0]) # chord
+        x_c[:,0] = (pts[:,0] - pts[ile,0]) / crd
+        x_c[:,1] = (pts[:,2] - pts[ile,2]) / crd
+        data = np.hstack((x_c, vals['Cp']))
         tU.sort(elems, data)
         data = np.vstack((data, data[0,:]))
-        tU.write(data, mshName+sfx+'_slice_'+str(i)+'.dat', '%1.5e', ', ', 'x, y, z, x/c, Cp', '')
+        hdr = f'y = {ys[i]}, c = {crd}\n'
+        hdr += '{:>9s}, {:>10s}, {:>10s}'.format('x/c', 'z/c', 'cp')
+        np.savetxt(mshName+sfx+'_slice_'+str(i)+'.dat', data, fmt='%+1.4e', delimiter=',', header=hdr)
diff --git a/dart/validation/agard.py b/dart/validation/agard.py
index 9d48d8e11e936401be2060413e40e6ad014cca8f..fd5452b80a5937b375c7a67909a8bd0b72c22527 100644
--- a/dart/validation/agard.py
+++ b/dart/validation/agard.py
@@ -83,9 +83,9 @@ def main():
     # post process
     tms['post'].start()
     if canPost:
-        floU.write_slices(msh.name,[0.01, 0.37, 0.75],5)
-        data = tbxU.read('agard445_slice_1.dat')
-        floD.plot(data[:,3], data[:,4], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': True})
+        floU.write_slices(msh.name,[0.01, 0.37, 0.75], 5)
+        data = np.loadtxt(f'{msh.name}_slice_1.dat', delimiter=',', skiprows=2)
+        floD.plot(data[:,0], data[:,-1], {'xlabel': 'x/c', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': True})
     tms['post'].stop()
 
     # display timers
diff --git a/dart/validation/crm.py b/dart/validation/crm.py
index 27ebf2bd417ff7e7aca0582a29d6f6498c1b17cb..12b0d0a2d05aabc823137886f5194175f503c74b 100644
--- a/dart/validation/crm.py
+++ b/dart/validation/crm.py
@@ -106,8 +106,8 @@ def main():
     if hasVTK:
         utils.write_slices(msh.name, [3.81973, 5.8765, 8.2271, 11.753, 14.6913, 17.6295, 19.0986, 21.4492, 24.9751, 27.91338], 12)
         for i in range(10):
-            data = np.loadtxt(f'{msh.name}_slice_{i}.dat', delimiter=',', skiprows=1) # read
-            plt.plot(data[:,3], data[:,4], lw=2) # plot results
+            data = np.loadtxt(f'{msh.name}_slice_{i}.dat', delimiter=',', skiprows=2) # read
+            plt.plot(data[:,0], data[:,-1], lw=2) # plot results
         font = {'family': 'serif', 'color': 'darkred', 'weight': 'normal', 'size': 16}
         plt.xlabel('$x/c$', fontdict=font)
         plt.ylabel('$C_p$', fontdict=font)
diff --git a/dart/validation/onera.py b/dart/validation/onera.py
index 8eb19ebc22912e9d473f18778d724f08630cb5f9..cc1c3d33b8d2e2bfbd50b96eb440addd9cb99253 100644
--- a/dart/validation/onera.py
+++ b/dart/validation/onera.py
@@ -81,9 +81,9 @@ def main():
     # post process
     tms['post'].start()
     if canPost:
-        floU.write_slices(msh.name,[0.01, 0.239, 0.526, 0.777, 0.957, 1.076, 1.136, 1.184],5)
-        data = tbxU.read('oneraM6_slice_2.dat')
-        floD.plot(data[:,3], data[:,4], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': True})
+        floU.write_slices(msh.name,[0.01, 0.239, 0.526, 0.777, 0.957, 1.076, 1.136, 1.184], 5)
+        data = np.loadtxt(f'{msh.name}_slice_2.dat', delimiter=',', skiprows=2)
+        floD.plot(data[:,0], data[:,-1], {'xlabel': 'x/c', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': True})
     tms['post'].stop()
 
     # display timers
diff --git a/vii/interfaces/dart/DartInterface.py b/vii/interfaces/dart/DartInterface.py
index 70b2514a78ae89a04a65a984161406420ff90ddd..c347fdb1160b2dd908f8635da87c25c23bae5540 100644
--- a/vii/interfaces/dart/DartInterface.py
+++ b/vii/interfaces/dart/DartInterface.py
@@ -54,7 +54,7 @@ class DartInterface:
         """ Write Cp distribution around the airfoil on file.
         """
         Cp = floU.extract(self.group[0].boundary.tag.elems, self.solver.Cp)
-        tboxU.write(Cp, 'Cp'+sfx+'.dat', '%1.5e',', ', 'x, y, z, Cp', '')
+        tboxU.write(Cp, 'Cp'+sfx+'.dat', '%1.4e', ',', 'x, y, cp', '')
     
     def save(self, writer):
         self.solver.save(writer)
diff --git a/vii/tests/bli2D.py b/vii/tests/bli2D.py
index e4456e222a8674c6d14904c26465161de6ca2f32..6d0534528a823080aca89438f6af7fa196434c65 100644
--- a/vii/tests/bli2D.py
+++ b/vii/tests/bli2D.py
@@ -84,7 +84,7 @@ def main():
     # Save results.
     iSolverAPI.save(gmshWritter)
     Cp = invUtils.extract(bnd.groups[0].tag.elems, iSolverAPI.solver.Cp)
-    tboxU.write(Cp, 'Cp.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    tboxU.write(Cp, 'Cp.dat', '%1.4e', ',', 'x, y, cp', '')
 
     print('')
     # Display results.
@@ -105,7 +105,7 @@ def main():
     print(Coupler.tms)
 
     # Plot results.
-    invDefault.plot(Cp[:,0], Cp[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(iSolverAPI.getCl(), vSolver.Cdt, iSolverAPI.getCm(), 4), 'invert': True})
+    invDefault.plot(Cp[:,0], Cp[:,-1], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(iSolverAPI.getCl(), vSolver.Cdt, iSolverAPI.getCm(), 4), 'invert': True})
     invDefault.plot(vSolution['x/c'], vSolution['Cf'], {'xlabel': '$x/c$', 'ylabel': '$c_f$', 'title': 'Cdt = {0:.{3}f}, Cdf = {1:.{3}f}, Cdp = {2:.{3}f}'.format(vSolver.Cdt, vSolver.Cdf, vSolver.Cdp, 4), 'xlim': [0, 1]})
 
     # Check results.