diff --git a/terrainlib_lua/rivermapper.lua b/terrainlib_lua/rivermapper.lua
index 38ebb7b..92ab11f 100644
--- a/terrainlib_lua/rivermapper.lua
+++ b/terrainlib_lua/rivermapper.lua
@@ -93,27 +93,6 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 	local links = {}
 	local basin_links
 
-	-- Function to analyse a link between two nodes
-	local function add_link(i1, i2, b1, isY)
-		-- i1, i2: coordinates of two nodes
-		-- b1: basin that contains i1
-		-- isY: whether the link is in Y direction
-		-- Note that basin number #0 represents the outside of the map; or if the coordinate is inside the map, means that the basin number is uninitialized.
-		local b2 = i2 == 0 and 0 or basin_id[i2]
-		local elev = i2 == 0 and dem[i1] or mmax(dem[i1], dem[i2]) -- Elevation of the highest of the two sides of the link (or only i1 if b2 is map outside)
-		local l2 = basin_links[b2]
-		if not l2 then
-			l2 = {b2, b1, elev=elev, i=mmax(i1,i2), is_y=isY}
-			basin_links[b2] = l2
-		elseif l2.elev > elev then -- If this link is lower than the lowest registered link between these two basins, register it as the new lowest pass
-			l2.elev = elev
-			l2.i = mmax(i1,i2)
-			l2.is_y = isY
-			l2[1] = b2
-			l2[2] = b1
-		end
-	end
-
 	for i=1, X*Y do
 		basin_id[i] = 0
 	end
@@ -139,10 +118,32 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 			-- If no river is coming from the East, we might be at the limit of two basins, thus we need to test adjacency.
 			elseif i%X > 0 then
 				if basin_id[i+1] ~= ib and basin_id[i+1] ~= 0 then
-					add_link(i, i+1, ib, false)
+					local b2 = basin_id[i+1]
+					local elev = mmax(dem[i], dem[i+1]) -- Elevation of the highest of the two sides of the link (or only i1 if b2 is map outside)
+					local l2 = basin_links[b2]
+					if not l2 then
+						l2 = {b2, ib, elev=elev, i=i+1, is_y=false}
+						basin_links[b2] = l2
+					elseif l2.elev > elev then -- If this link is lower than the lowest registered link between these two basins, register it as the new lowest pass
+						l2.elev = elev
+						l2.i = i+1
+						l2.is_y = false
+						l2[1] = b2
+						l2[2] = ib
+					end
 				end
 			else -- If the eastern neighbour is outside the map
-				add_link(i, 0, ib, false)
+				local l2 = basin_links[0]
+				if not l2 then
+					l2 = {0, ib, elev=dem[i], i=i, is_y=false}
+					basin_links[0] = l2
+				elseif l2.elev > dem[i] then
+					l2.elev = dem[i]
+					l2.i = i
+					l2.is_y = false
+					l2[1] = 0
+					l2[2] = ib
+				end
 			end
 
 			if d >= 4 then -- River coming from the South
@@ -151,10 +152,32 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 				singular[cur] = i+X
 			elseif i <= X*(Y-1) then
 				if basin_id[i+X] ~= ib and basin_id[i+X] ~= 0 then
-					add_link(i, i+X, ib, true)
+					local b2 = basin_id[i+X]
+					local elev = mmax(dem[i], dem[i+X])
+					local l2 = basin_links[b2]
+					if not l2 then
+						l2 = {b2, ib, elev=elev, i=i+X, is_y=true}
+						basin_links[b2] = l2
+					elseif l2.elev > elev then
+						l2.elev = elev
+						l2.i = i+X
+						l2.is_y = true
+						l2[1] = b2
+						l2[2] = ib
+					end
 				end
 			else
-				add_link(i, 0, ib, true)
+				local l2 = basin_links[0]
+				if not l2 then
+					l2 = {0, ib, elev=dem[i], i=i, is_y=true}
+					basin_links[0] = l2
+				elseif l2.elev > dem[i] then
+					l2.elev = dem[i]
+					l2.i = i
+					l2.is_y = true
+					l2[1] = 0
+					l2[2] = ib
+				end
 			end
 
 			if d >= 2 then -- River coming from the West
@@ -163,10 +186,32 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 				singular[cur] = i-1
 			elseif i%X ~= 1 then
 				if basin_id[i-1] ~= ib and basin_id[i-1] ~= 0 then
-					add_link(i, i-1, ib, false)
+					local b2 = basin_id[i-1]
+					local elev = mmax(dem[i], dem[i-1])
+					local l2 = basin_links[b2]
+					if not l2 then
+						l2 = {b2, ib, elev=elev, i=i, is_y=false}
+						basin_links[b2] = l2
+					elseif l2.elev > elev then
+						l2.elev = elev
+						l2.i = i
+						l2.is_y = false
+						l2[1] = b2
+						l2[2] = ib
+					end
 				end
 			else
-				add_link(i, 0, ib, false)
+				local l2 = basin_links[0]
+				if not l2 then
+					l2 = {0, ib, elev=dem[i], i=i, is_y=false}
+					basin_links[0] = l2
+				elseif l2.elev > dem[i] then
+					l2.elev = dem[i]
+					l2.i = i
+					l2.is_y = false
+					l2[1] = 0
+					l2[2] = ib
+				end
 			end
 
 			if d >= 1 then -- River coming from the North
@@ -174,10 +219,32 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 				singular[cur] = i-X
 			elseif i > X then
 				if basin_id[i-X] ~= ib and basin_id[i-X] ~= 0 then
-					add_link(i, i-X, ib, true)
+					local b2 = basin_id[i-X]
+					local elev = mmax(dem[i], dem[i-X])
+					local l2 = basin_links[b2]
+					if not l2 then
+						l2 = {b2, ib, elev=elev, i=i, is_y=true}
+						basin_links[b2] = l2
+					elseif l2.elev > elev then
+						l2.elev = elev
+						l2.i = i
+						l2.is_y = true
+						l2[1] = b2
+						l2[2] = ib
+					end
 				end
 			else
-				add_link(i, 0, ib, true)
+				local l2 = basin_links[0]
+				if not l2 then
+					l2 = {0, ib, elev=dem[i], i=i, is_y=true}
+					basin_links[0] = l2
+				elseif l2.elev > dem[i] then
+					l2.elev = dem[i]
+					l2.i = i
+					l2.is_y = true
+					l2[1] = 0
+					l2[2] = ib
+				end
 			end
 	end
 	dirs2 = nil