diff --git a/hysop/backend/host/python/operator/flowrate_correction.py b/hysop/backend/host/python/operator/flowrate_correction.py
index 2d6c3de245027fea453bdc35419105e9ccb7f9fb..9ed4274584801523e25207a038513728489a39ee 100644
--- a/hysop/backend/host/python/operator/flowrate_correction.py
+++ b/hysop/backend/host/python/operator/flowrate_correction.py
@@ -44,7 +44,7 @@ class PythonFlowRateCorrection(HostOperator):
         self.absorption_start = absorption_start
         if self.absorption_start is None:
             self.absorption_start = velocity.domain.end[-1]
-        ## TODO: Correction is taking an absorption start (see absorption operator) only in X dir.
+        # TODO: Correction is taking an absorption start (see absorption operator) only in X dir.
         super(PythonFlowRateCorrection, self).__init__(
             input_fields=input_fields, output_fields=output_fields,
             input_params=input_params, **kwds)
@@ -63,35 +63,38 @@ class PythonFlowRateCorrection(HostOperator):
         super(PythonFlowRateCorrection, self).discretize()
         self.dvelocity = self.input_discrete_tensor_fields[self.velocity]
         self.dvorticity = self.input_discrete_tensor_fields[self.vorticity]
-        topo = self.dvelocity[0].topology
-        mesh = topo.mesh
+        vtopo = self.dvelocity[0].topology
+        wtopo = self.dvorticity[0].topology
+        vmesh = vtopo.mesh
+        wmesh = wtopo.mesh
 
         # Compute volume and surface integration coefficients
-        spaceStep = mesh.space_step
-        lengths = topo.domain.length
+        spaceStep = vmesh.space_step
+        lengths = vtopo.domain.length
         self._inv_ds = 1. / np.prod(lengths[:-1])
-        self._inv_dvol = 1. / (lengths[0]*lengths[1]*(self.absorption_start-topo.domain.origin[-1]))
+        self._inv_dvol = 1. / (lengths[0]*lengths[1] *
+                               (self.absorption_start-vtopo.domain.origin[-1]))
         self.coeff_mean = np.prod(spaceStep) / np.prod(lengths)
 
         # Compute space coordinates from domain origin
-        x0 = mesh.global_origin[-1]
-        self.x_coord = mesh.local_mesh_coords[0] - x0
+        x0 = vmesh.global_origin[-1]
+        self.x_coord = vmesh.local_mesh_coords[0] - x0
 
         # Compute slices for volume and surface integration
-        gstart = mesh.local_origin.copy()
-        gstart[-1] = mesh.global_origin[-1]
-        if mesh.is_inside_local_domain(gstart):
-            ind4integ_v = [_ for _ in mesh.local_compute_slices]
-            ind4integ_v[-1] = mesh.point_local_indices(gstart)[-1]
+        gstart = vmesh.local_origin.copy()
+        gstart[-1] = vmesh.global_origin[-1]
+        if vmesh.is_inside_local_domain(gstart):
+            ind4integ_v = [_ for _ in vmesh.local_compute_slices]
+            ind4integ_v[-1] = vmesh.point_local_indices(gstart)[-1]
             self._ind4integ_v = tuple(ind4integ_v)
         else:
             self._ind4integ_v = None
         self._ds_in = np.prod(spaceStep[0:2])
-        ind4integ_w = [_ for _ in mesh.local_compute_slices]
-        box_end = mesh.local_origin.copy()
+        ind4integ_w = [_ for _ in wmesh.local_compute_slices]
+        box_end = wmesh.local_origin.copy()
         box_end[-1] = self.absorption_start
         ind4integ_w[-1] = slice(ind4integ_w[-1].start,
-                                mesh.point_local_indices(box_end)[-1],
+                                wmesh.point_local_indices(box_end)[-1],
                                 ind4integ_w[-1].step)
         self._ind4integ_w = tuple(ind4integ_w)
         self._ds_box = np.prod(spaceStep)
@@ -134,13 +137,13 @@ class PythonFlowRateCorrection(HostOperator):
 
         flowrate_req = self.flowrate() * self._inv_ds
         velo_shift = np.zeros((dv.dim, ))
-        velo_shift = flowrate_req - self.rates[0:dv.dim]
+        velo_shift = flowrate_req - self.rates[dv.dim-1::-1]
 
         # Shift velocity Vx
         dv.data[0][...] += velo_shift[2]
 
         # Update Vy and Vz with mean vorticity
-        dv.data[1][...] += velo_shift[1] + self.rates[dv.dim]*self.x_coord
+        dv.data[1][...] += velo_shift[1] + self.rates[dv.dim+2]*self.x_coord
         dv.data[2][...] += velo_shift[0] + self.rates[dv.dim+1]*self.x_coord
         if (self.ghost_exchanger is not None):
             self.ghost_exchanger()