Convex hull

In the previously posted flownet generation program, there is a seemingly insurmountable problem. After equipotential points for a given head have been extracted from between the solved finite difference grid points (here, using linear interpolation) comes the hard part: properly joining them. I have outlined how I approached the problem in the first post about the program. There, I did point out that in equipotential points for a given head, discovered to the extreme left and right of the grid, a jumpy line springs into existence because of the limited space the grid actually represents which obviously clips the isobars. However, at the time, I hadn't realized that there was an even more sinister problem with the program and the 'connect the dots' logic - one which doesn't, at this point, seem solvable because of the weird ways in which these equipotential lines naturally seem to be able to bend. Take for example the following output with
contourInterval=5
numDys=13
numDxs=16
DxsTillDam=2
DamWidthInDxs=4
reservoirHead=100
downstreamBoundaryHead=50

The problem is visible for the 15m isobar. Starting from underneath the dam, the point near the 16.29m grid node should have been the one to be joined first rather than the point between the 11.06m and 19.26m grid nodes. It is easy enough for us to say via a simple visual inspection. But the 'order in descending Y coordinate and then by shortest consecutive distance' logic falls flat in this case since clearly, the program is doing exactly that. The point actually being joined is clearly closer to the equipotential point on dam than the 'right' point which happens to be a little too left.
I thought of trying all orderings of the points to minimize total length but quickly concluded that it would be computationally unviable - the number of equipotential points for a head regularly crossed 20 when I was testing the program - that's 20!, or about 10^18, number of combinations to try for drawing each isobar. What's more, I wasn't sure it would even produce the desired result. So, I scrapped the idea.

Then, I researched ways to do this online. Something called 'Convex hull' came up. Basically, it is a way of drawing a polygon around a bunch of points while passing through the right subset of the said points so that no point is left outside the polygon. The trick is to test each possible line segment (defined by two unique points) for convexity i.e. whether all other points lie towards the same side (of the line segment); if they do, the line segment being tested is part of the polygon - the 'convex hull'. The idea behind this test - whether a point lies toward one side of a line or the other - comes from checking the point's y-coordinate against that the line's, at the point's x-coordinate. If all the points lie above (the idea of a point being 'above' a line can be extended from the obvious case where the line is horizontal, right until it isn't completely vertical) the line, their y-coordinates are going to be greater than the y-coordinates of the line segment at their(the points') x-coordinates. Anyway, this thread's first answer is what I implemented - but not into the Fortran program itself. I was skeptical you see. I couldn't quite picture how this thing worked or what it would look like in action. So, I first made a small JavaScript program using p5.js so that I could have something concrete to look at and potentially tweak.

Here is the code for it.
Then, it became clear that no amount of tweaking the convex hull algorithm could manage to connect the equipotential points 'properly'. There would be too many exceptions - non-convex points that would need to be connected, but would not be. Leaving such points out would produce smooth and nice equipotential lines but that would also be misleading and far from accurate.
So, there's where I basically left the flownet program, without progress.
But at least I got to learn something new in the process.