diff --git a/src/lib/grid.cpp b/src/lib/grid.cpp
index cbd7541..bbd82c6 100644
--- a/src/lib/grid.cpp
+++ b/src/lib/grid.cpp
@@ -186,7 +186,7 @@ void wrap_gwdt_computecosts(
 void wrap_flow_routing_d8_carve(
         py::array_t<ptrdiff_t> node, py::array_t<uint8_t> direction,
         py::array_t<float> dem, py::array_t<float> dist,
-        py::array_t<int32_t> flats, std::tuple<ptrdiff_t, ptrdiff_t> dims){
+        py::array_t<int32_t> flats, std::tuple<ptrdiff_t, ptrdiff_t> dims, int order){
 
     ptrdiff_t *node_ptr = node.mutable_data();
     uint8_t *direction_ptr = direction.mutable_data();
@@ -197,8 +197,6 @@ void wrap_flow_routing_d8_carve(
     std::array<ptrdiff_t, 2> dims_array = {std::get<0>(dims), std::get<1>(dims)};
     ptrdiff_t *dims_ptr = dims_array.data();
 
-    unsigned int order = (direction.flags() & py::array::c_style) ? 1 : 0;
-
     flow_routing_d8_carve(node_ptr, direction_ptr, dem_ptr, dist_ptr, flats_ptr, dims_ptr, order);
 }
 
@@ -215,7 +213,7 @@ ptrdiff_t wrap_flow_routing_d8_edgelist(
         py::array_t<ptrdiff_t> target,
         py::array_t<ptrdiff_t> node,
         py::array_t<uint8_t> direction,
-        std::tuple<ptrdiff_t,ptrdiff_t> dims) {
+        std::tuple<ptrdiff_t,ptrdiff_t> dims, int order) {
 
     ptrdiff_t *source_ptr = source.mutable_data();
     ptrdiff_t *target_ptr = target.mutable_data();
@@ -226,8 +224,6 @@ ptrdiff_t wrap_flow_routing_d8_edgelist(
     std::array<ptrdiff_t, 2> dims_array = {std::get<0>(dims), std::get<1>(dims)};
     ptrdiff_t *dims_ptr = dims_array.data();
 
-    unsigned int order = (direction.flags() & py::array::c_style) ? 1 : 0;
-
     return flow_routing_d8_edgelist(source_ptr, target_ptr, node_ptr, direction_ptr, dims_ptr, order);
 }
 
diff --git a/src/topotoolbox/flow_object.py b/src/topotoolbox/flow_object.py
index db93b86..a5bc5dd 100644
--- a/src/topotoolbox/flow_object.py
+++ b/src/topotoolbox/flow_object.py
@@ -134,18 +134,20 @@ def _d8_carve(grid: GridObject,
     process, which could lead to issues when using very large DEMs.
     """
     dims = grid.dims
+    order = (0 if grid.z.flags.f_contiguous else 1)
+
 
     (aux, filled_dem, flats) = grid.auxiliary_topography(bc, hybrid)
 
     node = np.zeros_like(aux, dtype=np.int64)  # node: dtype=np.int64
     direction = np.zeros_like(aux, dtype=np.uint8)
-    _grid.flow_routing_d8_carve(node, direction, filled_dem, aux, flats, dims)
+    _grid.flow_routing_d8_carve(node, direction, filled_dem, aux, flats, dims, order)
 
     # ravel is used here to flatten the arrays. The memory order should not matter
     # because we only need a block of contiguous memory interpreted as a 1D array.
     source = np.zeros(aux.size, dtype=np.int64)  # source: dtype=int64
     target = np.zeros(aux.size, dtype=np.int64)       # target: dtype=int64
-    edge_count = _grid.flow_routing_d8_edgelist(source, target, node, direction, dims)
+    edge_count = _grid.flow_routing_d8_edgelist(source, target, node, direction, dims, order)
 
     return (direction,
             node,