Index: DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs
===================================================================
diff -u -r5268 -r5276
--- DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs	(.../PlLinesToWaternetConverter.cs)	(revision 5268)
+++ DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs	(.../PlLinesToWaternetConverter.cs)	(revision 5276)
@@ -29,6 +29,9 @@
 using Deltares.DamEngine.Data.Geometry;
 using Deltares.DamEngine.Data.Geotechnics;
 using Deltares.DamEngine.Data.Standard;
+using Deltares.MacroStability.CSharpWrapper.Input;
+using CharacteristicPointType = Deltares.DamEngine.Data.Geotechnics.CharacteristicPointType;
+using SoilProfile = Deltares.DamEngine.Data.Geotechnics.SoilProfile;
 
 namespace Deltares.DamEngine.Calculators.KernelWrappers.MacroStabilityCommon;
 
@@ -557,43 +560,36 @@
             }
 
             List<Point2D> intersectionPoints = outerLoop.IntersectionXzPointsWithGeometryString(plLineInGeometryPoint);
-            // Only the layers that are below and do not intersect the phreatic line are kept.
-            if ((intersectionPoints.Count == 0) && (outerLoop[0].Z < plLineInGeometryPoint.GetZatX(outerLoop[0].X)))
+            // Only the layers that intersect the phreatic line are kept.
+            if (intersectionPoints.Count > 0)
             {
                 var geometrySurface = new GeometrySurface(outerLoop);
-                GeometryPointString top = geometrySurface.DetermineTopGeometrySurface();
-                top.SyncPoints();
-                top.SortPointsByXAscending();
-                curves.Add(top);
-            }
-
-            // For lowest layers that intersect the phreatic line the bottom of these layers is kept
-            if (intersectionPoints.Count > 0 && outerLoop.Points.Any(p => Math.Abs(p.Z - soilProfile2D.GetSoilProfile1D(p.X).BottomLevel) < toleranceAlmostEqual))
-            {
-                var geometrySurface = new GeometrySurface(outerLoop);
                 GeometryPointString bottom = geometrySurface.DetermineBottomGeometrySurface();
                 bottom.SyncPoints();
                 bottom.SortPointsByXAscending();
                 curves.Add(bottom);
             }
+            
+            // Add also the surface line in case of free water
+            curves.Add(soilProfile2D.Geometry.SurfaceLine);
         }
 
-        // The waternet line is the polyline formed with the highest lines connected to each other.
+        // The waternet line is the polyline formed with the lowest lines connected to each other.
         List<Line> lines = DivideCurvesIntoLines(curves);
         IEnumerable<Point2D> allPoints = curves.SelectMany(curve => curve.CalcPoints);
         double[] xCoordinates = allPoints.Select(point => point.X).OrderBy(x => x).Distinct().ToArray();
         List<Line> splitLines = SplitLinesAtXCoordinates(xCoordinates, lines);
         for (var i = 0; i < xCoordinates.Length - 1; i++)
         {
-            List<Line> relevantLines = splitLines.FindAll(line => Math.Abs(line.BeginPoint.X - xCoordinates[i]) < toleranceAlmostEqual && Math.Abs(line.BeginPoint.X - line.EndPoint.X) > 0);
-            double maxZBegin = relevantLines.Select(line => line.BeginPoint.Z).Max();
-            List<Line> relevantLines2 = relevantLines.FindAll(line => Math.Abs(line.BeginPoint.Z - maxZBegin) < toleranceAlmostEqual);
-            double maxZEnd = relevantLines2.Select(line => line.EndPoint.Z).Max();
-            Line relevantLine = relevantLines2.Find(line => Math.Abs(line.EndPoint.Z - maxZEnd) < toleranceAlmostEqual);
+            List<Line> relevantLinesAtX = splitLines.FindAll(line => Math.Abs(line.BeginPoint.X - xCoordinates[i]) < toleranceAlmostEqual && Math.Abs(line.BeginPoint.X - line.EndPoint.X) > 0);
+            double minZBegin = relevantLinesAtX.Select(line => line.BeginPoint.Z).Min();
+            List<Line> relevantLinesAtXWithZBeginMax = relevantLinesAtX.FindAll(line => Math.Abs(line.BeginPoint.Z - minZBegin) < toleranceAlmostEqual);
+            double minZEnd = relevantLinesAtXWithZBeginMax.Select(line => line.EndPoint.Z).Min();
+            Line relevantLine = relevantLinesAtXWithZBeginMax.Find(line => Math.Abs(line.EndPoint.Z - minZEnd) < toleranceAlmostEqual);
             waternetLine.Points.Add(new GeometryPoint(relevantLine.BeginPoint.X, relevantLine.BeginPoint.Z));
             waternetLine.Points.Add(new GeometryPoint(relevantLine.EndPoint.X, relevantLine.EndPoint.Z));
         }
-        
+
         // Make sure the waternet line for PL 1 is not in the bottom aquifer.
         foreach (GeometryPoint waternetLinePoint in waternetLine.Points)
         {