diff --git a/README.md b/README.md
index 16c14f0..8166707 100644
--- a/README.md
+++ b/README.md
@@ -78,6 +78,7 @@ Other parameters can be specified by `--parameter value`. Syntax `--parameter=va
 | `niter`       | Number of iterations. Each iteration represents a time `time/niter`. | `--niter 10` |
 | `sea_level_variations` | Amplitude of sea level variations throughout the simulation (if any). | `--sea_level_variations 10` |
 | `sea_level_variations_time` | Characteristic time of variation for sea level, in the same units than `time`. Increasing it will result in slower variations between iterations. | `--sea_level_variations_time 1` |
+| `flow_method` | Algorithm used for local flow calculation. Possible values are `steepest` (every node flows toward the steepest neighbour when possible), and `semirandom` (default, flow direction is determined randomly between lower neighbours, with lowest ones having greater probability). | `--flow_method semirandom` |
 | | **Alternatives** |
 | `config`      | Another way to specify configuration file | `--config terrain_higher.conf` |
 | `output`      | Another way to specify output dir | `--output ~/.minetest/worlds/my_world/river_data` |
diff --git a/generate.py b/generate.py
index 9d28f95..783d5a5 100755
--- a/generate.py
+++ b/generate.py
@@ -87,6 +87,7 @@ sea_level = float(get_setting('sea_level', 0.0))
 sea_level_variations = float(get_setting('sea_level_variations', 0.0))
 sea_level_variations_time = float(get_setting('sea_level_variations_time', 1.0))
 flex_radius = float(get_setting('flex_radius', 20.0))
+flow_method = get_setting('flow_method', 'semirandom')
 
 time = float(get_setting('time', 10.0))
 niter = int(get_setting('niter', 10))
@@ -120,7 +121,7 @@ n = noisemap(mapsize+1, mapsize+1, **params)
 ### COMPUTE LANDSCAPE EVOLUTION
 # Initialize landscape evolution model
 print('Initializing model')
-model = terrainlib.EvolutionModel(n, K=K, m=m, d=d, sea_level=sea_level, flex_radius=flex_radius)
+model = terrainlib.EvolutionModel(n, K=K, m=m, d=d, sea_level=sea_level, flex_radius=flex_radius, flow_method=flow_method)
 terrainlib.update(model.dem, model.lakes, t=5, sea_level=model.sea_level, title='Initializing...')
 
 dt = time/niter
diff --git a/terrainlib/erosion.py b/terrainlib/erosion.py
index 6547640..9d2ed89 100644
--- a/terrainlib/erosion.py
+++ b/terrainlib/erosion.py
@@ -54,7 +54,7 @@ def diffusion(dem, time, d=1):
     return im.gaussian_filter(dem, radius, mode='reflect') # Diffusive erosion is a simple Gaussian blur
 
 class EvolutionModel:
-    def __init__(self, dem, K=1, m=0.5, d=1, sea_level=0, flow=False, flex_radius=100):
+    def __init__(self, dem, K=1, m=0.5, d=1, sea_level=0, flow=False, flex_radius=100, flow_method='semirandom'):
         self.dem = dem
         #self.bedrock = dem
         self.K = K
@@ -63,6 +63,8 @@ class EvolutionModel:
         self.sea_level = sea_level
         self.flex_radius = flex_radius
         self.define_isostasy()
+        self.flow_method = flow_method
+        #set_flow_method(flow_method)
         if flow:
             self.calculate_flow()
         else:
@@ -72,7 +74,7 @@ class EvolutionModel:
             self.flow_uptodate = False
 
     def calculate_flow(self):
-        self.dirs, self.lakes, self.rivers = flow(self.dem)
+        self.dirs, self.lakes, self.rivers = flow(self.dem, method=self.flow_method)
         self.flow_uptodate = True
 
     def advection(self, time):
diff --git a/terrainlib/rivermapper.py b/terrainlib/rivermapper.py
index 2de45b5..3c93bd5 100644
--- a/terrainlib/rivermapper.py
+++ b/terrainlib/rivermapper.py
@@ -12,7 +12,19 @@ from collections import defaultdict
 # The algorithm here makes use of most of the paper's concepts, including the Planar Boruvka algorithm.
 # Only flow_local and accumulate_flow are custom algorithms.
 
-def flow_local(plist):
+# Define two different method for local flow routing
+def flow_local_steepest(plist):
+    vmax = 0.0
+    imax = 0.0
+    for i, p in enumerate(plist):
+        if p > vmax:
+            vmax = p
+            imax = i
+    if vmax > 0.0:
+        return imax+1
+    return 0
+
+def flow_local_semirandom(plist):
     """
     Determines a flow direction based on denivellation for every neighbouring node.
     Denivellation must be positive for downward and zero for flat or upward:
@@ -27,7 +39,16 @@ def flow_local(plist):
             return i+1
         r -= p
 
-def flow(dem):
+flow_local_methods = {
+    'steepest' : flow_local_steepest,
+    'semirandom' : flow_local_semirandom,
+}
+
+def flow(dem, method='semirandom'):
+    if method in flow_local_methods:
+        flow_local = flow_local_methods[method]
+    else:
+        raise KeyError('Flow method \'{}\' does not exist'.format(method))
 
     # Flow locally
     dirs1 = np.zeros(dem.shape, dtype=int)