diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/find1dcells.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/find1dcells.f90
index 8b7d9301df..4cbcad8b03 100644
--- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/find1dcells.f90
+++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/find1dcells.f90
@@ -49,7 +49,8 @@ module m_find1dcells
 !>    it is assumed that kc has been allocated
 !>    it is assumed that findcells has already been called (for 2d cells)
    subroutine find1dcells()
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
 
       implicit none
 
diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_pol_to_cellmask.F90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_pol_to_cellmask.F90
new file mode 100644
index 0000000000..2891fab09e
--- /dev/null
+++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_pol_to_cellmask.F90
@@ -0,0 +1,159 @@
+ !----- AGPL --------------------------------------------------------------------
+!
+!  Copyright (C)  Stichting Deltares, 2017-2026.
+!
+!  This file is part of Delft3D (D-Flow Flexible Mesh component).
+!
+!  Delft3D is free software: you can redistribute it and/or modify
+!  it under the terms of the GNU Affero General Public License as
+!  published by the Free Software Foundation version 3.
+!
+!  Delft3D  is distributed in the hope that it will be useful,
+!  but WITHOUT ANY WARRANTY; without even the implied warranty of
+!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+!  GNU Affero General Public License for more details.
+!
+!  You should have received a copy of the GNU Affero General Public License
+!  along with Delft3D.  If not, see <http://www.gnu.org/licenses/>.
+!
+!  contact: delft3d.support@deltares.nl
+!  Stichting Deltares
+!  P.O. Box 177
+!  2600 MH Delft, The Netherlands
+!
+!  All indications and logos of, and references to, "Delft3D",
+!  "D-Flow Flexible Mesh" and "Deltares" are registered trademarks of Stichting
+!  Deltares, and remain the property of Stichting Deltares. All rights reserved.
+!
+!-------------------------------------------------------------------------------
+
+!> Wrapper around cellmask_from_polygon_set that uses OpenMP to parallelize the loop over all points if not in MPI mode
+module m_pol_to_cellmask
+   use precision, only: dp
+   use m_cellmask_from_polygon_set, only: cellmask_from_polygon_set_init, cellmask_from_polygon_set_cleanup, cellmask_from_polygon_set
+
+   implicit none
+
+   private
+
+   public :: pol_to_cellmask, init_cell_geom_as_polylines, find_cells_crossed_by_polyline
+
+
+contains
+
+   subroutine pol_to_cellmask()
+      use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl
+      use m_alloc, only: realloc
+
+      integer :: k
+
+      if (NPL == 0) then
+         return
+      end if
+
+      call realloc(cellmask, nump1d2d, keepexisting=.false., fill=0)
+
+      call cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl)
+
+      !> Dynamic scheduling in case of unequal work, chunksize guided
+      !$OMP PARALLEL DO SCHEDULE(GUIDED)
+      do k = 1, nump
+         cellmask(k) = cellmask_from_polygon_set(xzw(k), yzw(k))
+      end do
+      !$OMP END PARALLEL DO
+
+      call cellmask_from_polygon_set_cleanup()
+
+   end subroutine pol_to_cellmask
+
+!> Initialize xpl, ypl, zpl arrays with all netcell geometries (called once)
+   subroutine init_cell_geom_as_polylines()
+      use network_data
+      use m_alloc
+      USE m_cellmask_from_polygon_set, ONLY: cellmask_from_polygon_set_init, cellmask_initialized
+      use m_missing, only: dmiss
+
+      integer :: k, n, k1, total_points, ipoint
+
+      if (cellmask_initialized) then !> reuse cellmask cache boolean
+         return
+      end if
+
+      call savepol()
+
+      ! calculate total points needed: sum(netcell(k)%n + 1) for all cells
+      ! +1 for dmiss separator after each polygon
+      total_points = 0
+      do k = 1, nump
+         total_points = total_points + netcell(k)%n + 1 ! +1 for dmiss
+      end do
+
+      ! allocate or reallocate xpl, ypl, zpl
+      call realloc(xpl, total_points, keepexisting=.false.)
+      call realloc(ypl, total_points, keepexisting=.false.)
+      call realloc(zpl, total_points, keepexisting=.false.)
+
+      ! fill arrays with netcell geometry
+      ipoint = 0
+      do k = 1, nump
+         do n = 1, netcell(k)%n
+            ipoint = ipoint + 1
+            k1 = netcell(k)%nod(n)
+            xpl(ipoint) = xk(k1)
+            ypl(ipoint) = yk(k1)
+            zpl(ipoint) = real(k, dp) ! store cell index as z-value
+         end do
+
+         ! add separator
+         ipoint = ipoint + 1
+         xpl(ipoint) = dmiss
+         ypl(ipoint) = dmiss
+         zpl(ipoint) = dmiss
+      end do
+
+      npl = ipoint
+
+      ! initialize the cellmask module with these polygons
+      ! this builds bounding boxes and polygon indices
+      call cellmask_from_polygon_set_init(npl, xpl, ypl, zpl)
+
+   end subroutine init_cell_geom_as_polylines
+
+!> Find all cells crossed by polyline using brute force on cached geometry. The routine is inclusive of edge cases (touching edges or vertices).
+   subroutine find_cells_crossed_by_polyline(xpoly, ypoly, crossed_cells, error)
+      use m_alloc, only: realloc
+      use network_data, only: cellmask, nump
+      use m_missing, only: dmiss
+      use m_cellmask_from_polygon_set, only: cellmask_initialized, find_cells_for_segment
+
+      real(kind=dp), intent(in) :: xpoly(:) !< Polyline x-coordinates
+      real(kind=dp), intent(in) :: ypoly(:) !< Polyline y-coordinates
+      integer, allocatable, intent(out) :: crossed_cells(:) !> Indices of crossed cells in network_data::netcells
+      character, dimension(:), allocatable, intent(out) :: error !> Error message, empty if no error, to be handled at call site
+
+      integer :: npoly, i
+
+      error = ''
+
+      npoly = size(xpoly)
+      if (any(xpoly == dmiss) .or. any(ypoly == dmiss) .or. npoly < 2) then
+         error = 'Invalid polyline input'
+         return
+      end if
+
+      if (.not. cellmask_initialized) then
+         call init_cell_geom_as_polylines()
+      end if
+
+      call realloc(cellmask, nump, keepexisting=.false., fill=0)
+
+      ! Process each segment and put the result in cellmask
+      do i = 1, npoly - 1
+         call find_cells_for_segment(xpoly(i), ypoly(i), xpoly(i + 1), ypoly(i + 1), cellmask)
+      end do
+
+      crossed_cells = pack([(i, i=1, nump)], mask=(cellmask == 1))
+
+   end subroutine find_cells_crossed_by_polyline
+
+end module m_pol_to_cellmask
diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_gui/m_samples_to_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_samples_to_cellmask.f90
similarity index 94%
rename from src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_gui/m_samples_to_cellmask.f90
rename to src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_samples_to_cellmask.f90
index cce9ce5c07..47ded9a5e0 100644
--- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_gui/m_samples_to_cellmask.f90
+++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_samples_to_cellmask.f90
@@ -39,7 +39,8 @@ subroutine samples_to_cellmask()
 
       use network_data, only: cellmask, nump1d2d
       use m_samples, only: ns, xs, ys
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
       use m_alloc, only: realloc
 
       integer :: i, k
diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90
deleted file mode 100644
index 53aa7697dc..0000000000
--- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90
+++ /dev/null
@@ -1,68 +0,0 @@
- !----- AGPL --------------------------------------------------------------------
-!
-!  Copyright (C)  Stichting Deltares, 2017-2026.
-!
-!  This file is part of Delft3D (D-Flow Flexible Mesh component).
-!
-!  Delft3D is free software: you can redistribute it and/or modify
-!  it under the terms of the GNU Affero General Public License as
-!  published by the Free Software Foundation version 3.
-!
-!  Delft3D  is distributed in the hope that it will be useful,
-!  but WITHOUT ANY WARRANTY; without even the implied warranty of
-!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-!  GNU Affero General Public License for more details.
-!
-!  You should have received a copy of the GNU Affero General Public License
-!  along with Delft3D.  If not, see <http://www.gnu.org/licenses/>.
-!
-!  contact: delft3d.support@deltares.nl
-!  Stichting Deltares
-!  P.O. Box 177
-!  2600 MH Delft, The Netherlands
-!
-!  All indications and logos of, and references to, "Delft3D",
-!  "D-Flow Flexible Mesh" and "Deltares" are registered trademarks of Stichting
-!  Deltares, and remain the property of Stichting Deltares. All rights reserved.
-!
-!-------------------------------------------------------------------------------
-
-!> Wrapper around cellmask_from_polygon_set that uses OpenMP to parallelize the loop over all points if not in MPI mode
-module m_pol_to_cellmask
-   use precision, only: dp
-   use m_cellmask_from_polygon_set, only: cellmask_from_polygon_set_init, cellmask_from_polygon_set_cleanup, cellmask_from_polygon_set
-
-   implicit none
-
-   private
-
-   public :: pol_to_cellmask
-
-contains
-
-   subroutine pol_to_cellmask()
-      use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl
-      use m_alloc, only: realloc
-
-      integer :: k
-
-      if (NPL == 0) then
-         return
-      end if
-
-      call realloc(cellmask, nump1d2d, keepexisting=.false., fill=0)
-
-      call cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl)
-
-      !> Dynamic scheduling in case of unequal work, chunksize guided
-      !$OMP PARALLEL DO SCHEDULE(GUIDED)
-      do k = 1, nump
-         cellmask(k) = cellmask_from_polygon_set(xzw(k), yzw(k))
-      end do
-      !$OMP END PARALLEL DO
-
-      call cellmask_from_polygon_set_cleanup()
-
-   end subroutine pol_to_cellmask
-
-end module m_pol_to_cellmask
diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/timespace/fm_external_forcings.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/timespace/fm_external_forcings.f90
index a5209cb350..0ab5d6b7a4 100644
--- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/timespace/fm_external_forcings.f90
+++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/timespace/fm_external_forcings.f90
@@ -2509,7 +2509,8 @@ subroutine finalize()
       use m_filez, only: doclose
       use m_physcoef, only: constant_dicoww, dicoww
       use m_array_or_scalar, only: realloc
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
 
       integer :: j, k, ierr, l, n, itp, kk, k1, k2, kb, kt, nstor, i, ja
       integer :: imba, needextramba, needextrambar
diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90
index 7db987e4a2..5f248c4276 100644
--- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90
+++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90
@@ -3,7 +3,8 @@ module test_pol_to_cellmask
    use precision, only: dp
    use m_missing, only: dmiss
    use network_data, only: cellmask, npl, nump, xzw, yzw, xpl, ypl, zpl
-   use m_cellmask_from_polygon_set, only: cellmask_from_polygon_set_init, cellmask_from_polygon_set, cellmask_from_polygon_set_cleanup
+   use m_cellmask_from_polygon_set
+   use m_pol_to_cellmask
    use geometry_module, only: pinpok_legacy, pinpok_raycast
 
    implicit none(external)
@@ -406,7 +407,8 @@ end subroutine test_pinpok_raycast_complex
    subroutine test_incells_basic_functionality() bind(C)
       ! Test basic incells functionality: point inside/outside netcells
       use gridoperations, only: incells
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
       use network_data, only: netcell, nump, xk, yk
       use m_alloc, only: realloc
 
@@ -471,7 +473,8 @@ end subroutine test_incells_basic_functionality
    subroutine test_incells_complex_geometry() bind(C)
       ! Test incells with more complex netcell geometries (triangles, pentagons, hexagons)
       use gridoperations, only: incells
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
       use network_data, only: netcell, nump, xk, yk
 
       integer :: kin_old, kin_new
@@ -524,7 +527,8 @@ end subroutine test_incells_complex_geometry
    subroutine test_incells_large_grid() bind(C)
       ! Test incells with a larger grid (performance sanity check)
       use gridoperations, only: incells
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
       use network_data, only: netcell, nump
 
       integer :: kin_old, kin_new, i, mismatches
@@ -566,7 +570,8 @@ end subroutine test_incells_large_grid
    subroutine test_incells_cache_consistency() bind(C)
       ! Test that cache initialization produces consistent results
       use gridoperations, only: incells
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
       use network_data, only: nump
 
       integer :: kin1, kin2, kin_old
@@ -614,7 +619,8 @@ end subroutine test_incells_cache_consistency
    subroutine test_incells_edge_cases() bind(C)
       ! Test edge cases: empty grid, single cell, point at infinity
       use gridoperations, only: incells
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
       use network_data, only: netcell, nump
 
       integer :: kin_old, kin_new
@@ -819,7 +825,6 @@ end subroutine test_samples_to_cellmask_edge_cases
 !$f90tw TESTCODE(TEST, test_pol_to_cellmask, test_find_cells_crossed_by_polyline_simple, test_find_cells_crossed_by_polyline_simple,
    subroutine test_find_cells_crossed_by_polyline_simple() bind(C)
       ! Test find_cells_crossed_by_polyline: polyline crosses some cells, misses others
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, find_cells_crossed_by_polyline, cleanup_cell_geom_polylines
       use network_data, only: nump
       use m_alloc, only: realloc
 
@@ -888,7 +893,6 @@ end subroutine test_find_cells_crossed_by_polyline_simple
 !$f90tw TESTCODE(TEST, test_pol_to_cellmask, test_find_cells_crossed_by_polyline_edge_cases, test_find_cells_crossed_by_polyline_edge_cases,
    subroutine test_find_cells_crossed_by_polyline_edge_cases() bind(C)
       ! Test find_cells_crossed_by_polyline: polyline crosses some cells, misses others
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, find_cells_crossed_by_polyline, cleanup_cell_geom_polylines
       use network_data, only: nump
       use m_alloc, only: realloc
 
@@ -952,6 +956,50 @@ end subroutine test_find_cells_crossed_by_polyline_edge_cases
 ! Helper subroutines for setting up test geometries
 ! ============================================================================
 
+!$f90tw TESTCODE(TEST, test_pol_to_cellmask, test_dbpinpol_basic_functionality, test_dbpinpol_basic_functionality,
+   subroutine test_dbpinpol_basic_functionality() bind(C)
+      ! Test dbpinpol with known inputs/outputs (no legacy comparison)
+      use geometry_module, only: dbpinpol
+      use m_polygon, only: xpl, ypl, zpl, npl
+      use m_alloc, only: realloc
+
+      real(kind=dp) :: xp, yp
+      integer :: in_result
+
+      ! Setup simple rectangular polygon
+      npl = 5
+      call realloc(xpl, npl, keepexisting=.false.)
+      call realloc(ypl, npl, keepexisting=.false.)
+      call realloc(zpl, npl, keepexisting=.false.)
+
+      xpl = [0.0_dp, 10.0_dp, 10.0_dp, 0.0_dp, 0.0_dp]
+      ypl = [0.0_dp, 0.0_dp, 10.0_dp, 10.0_dp, 0.0_dp]
+      zpl = [1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp]
+
+      ! Initialize AND test first point
+      in_result = -1
+      xp = 0.0_dp
+      yp = 0.0_dp
+      call dbpinpol(xp, yp, in_result, dmiss, 1, npl, xpl, ypl, zpl)
+      call f90_expect_eq(in_result, 1, "Point (0,0) should be inside polygon")
+
+      ! Test second point (reusing cached data)
+      xp = 5.0_dp
+      yp = 5.0_dp
+      call dbpinpol(xp, yp, in_result, dmiss, 1, npl, xpl, ypl, zpl)
+      call f90_expect_eq(in_result, 1, "Point (5,5) should be inside polygon")
+
+      ! Test outside point
+      xp = 15.0_dp
+      yp = 15.0_dp
+      call dbpinpol(xp, yp, in_result, dmiss, 1, npl, xpl, ypl, zpl)
+      call f90_expect_eq(in_result, 0, "Point (15,15) should be outside polygon")
+
+      deallocate (xpl, ypl, zpl)
+
+   end subroutine
+!$f90tw)
+
    ! Helper subroutine: setup 5 cells in a row (0-10, 10-20, 20-30, 30-40, 40-50)
    subroutine setup_row_netcells()
       use network_data, only: netcell, nump, xk, yk, numk, nump1d2d
diff --git a/src/utils_lgpl/deltares_common/packages/deltares_common/src/geometry_module.f90 b/src/utils_lgpl/deltares_common/packages/deltares_common/src/geometry_module.f90
index c6171c35f2..63bea731ea 100644
--- a/src/utils_lgpl/deltares_common/packages/deltares_common/src/geometry_module.f90
+++ b/src/utils_lgpl/deltares_common/packages/deltares_common/src/geometry_module.f90
@@ -61,7 +61,6 @@ module geometry_module
    public :: cart3Dtospher
    public :: dbpinpol
    public :: cross
-   public :: dbpinpol_optinside_perpol
    public :: get_startend
    public :: pinpok3D
    public :: cross3D
@@ -101,6 +100,17 @@ module geometry_module
       module procedure clockwise_hp
    end interface clockwise
 
+   ! Interface for dbpinpol (implementation in submodule)
+   interface
+      module subroutine dbpinpol(xp, yp, in, dmiss, JINS, NPL, xpl, ypl, zpl)
+         real(kind=dp), intent(in) :: xp, yp
+         integer, intent(inout) :: in
+         real(kind=dp), intent(in) :: dmiss
+         integer, intent(in) :: JINS, NPL
+         real(kind=dp), optional, intent(in) :: xpl(:), ypl(:), zpl(:)
+      end subroutine dbpinpol
+   end interface
+
 contains
 
    !> projects a point to a polyline and finds the closest link.
@@ -275,76 +285,69 @@ pure subroutine pinpok(xl, yl, n, x, y, inside, jins_dummy, dmiss_dummy)
 
 !> optimized ray-casting point-in-polygon test.
 !! pure function that works with array slices or full arrays.
-   pure function pinpok_raycast(xl, yl, x, y, n) result(is_inside)
+pure function pinpok_raycast(xl, yl, x, y, n) result(is_inside)
 
-      real(kind=dp), intent(in) :: xl, yl !< point coordinates to test
-      integer, intent(in) :: n !< number of polygon points
-      real(kind=dp), dimension(n), intent(in) :: x, y !< polygon coordinates (at least n elements)
-      logical :: is_inside !< result: true if inside (respecting jins mode)
+    real(kind=dp), intent(in) :: xl, yl !< point coordinates to test
+    integer, intent(in) :: n !< number of polygon points
+    real(kind=dp), dimension(n), intent(in) :: x, y !< polygon coordinates (at least n elements)
+    logical :: is_inside !< result: true if inside (respecting jins mode)
 
-      ! locals
-      integer :: i, j, crossings
-      real(kind=dp) :: x_intersect
+    ! locals
+    integer :: i, j, crossings
+    real(kind=dp) :: x_intersect
+    real(kind=dp), parameter :: edge_tolerance = 1.0e-10_dp
 
-      is_inside = .false.
+    is_inside = .false.
 
-      ! degenerate polygon check
-      if (n <= 2) then
-         is_inside = .true.
-         if (jins == 0) then
+    ! degenerate polygon check
+    if (n <= 2) then
+        is_inside = .true.
+        if (jins == 0) then
             is_inside = .not. is_inside
-         end if
-         return
-      end if
+        end if
+        return
+    end if
 
-      ! ray-casting algorithm: count crossings of horizontal ray from point to +infinity
-      crossings = 0
+    ! ray-casting algorithm: count crossings of horizontal ray from point to +infinity
+    crossings = 0
       j = n ! In order to check the entire polygon, start the first link which lies between the last point (n) and the first point
 
-      do i = 1, n
-         ! check for missing value (polygon separator)
-         if (x(i) == dmiss) then
+    do i = 1, n
+        ! check for missing value (polygon separator)
+        if (x(i) == dmiss) then
             exit
-         end if
-
-         ! check if point is exactly on a vertex
-         if (xl == x(j) .and. yl == y(j)) then
-            is_inside = .true.
-            if (jins == 0) then
-               is_inside = .not. is_inside
-            end if
-            return
-         end if
+        end if
 
-         ! check if ray crosses this edge
-         ! edge crosses horizontal line through test point if one endpoint is above and one below
-         if ((y(j) > yl) .neqv. (y(i) > yl)) then
+        ! check if ray crosses this edge
+        ! edge crosses horizontal line through test point if one endpoint is above and one below
+        if ((y(j) > yl) .neqv. (y(i) > yl)) then
             ! compute x-coordinate of edge-ray intersection
             x_intersect = x(j) + (yl - y(j)) * (x(i) - x(j)) / (y(i) - y(j))
-            if (xl < x_intersect) then
-               ! ray crosses edge to the right of point
-               crossings = crossings + 1
-            else if (xl == x_intersect) then
-               ! point is exactly on the edge
-               is_inside = .true.
-               if (jins == 0) then
-                  is_inside = .not. is_inside
-               end if
-               return
+            
+            if (xl < x_intersect - edge_tolerance) then
+                ! ray crosses edge to the right of point
+                crossings = crossings + 1
+            else if (abs(xl - x_intersect) <= edge_tolerance) then
+                ! point is on or very near the edge (including vertices)
+                is_inside = .true.
+                if (jins == 0) then
+                    is_inside = .not. is_inside
+                end if
+                return
             end if
-         end if
+        end if
          j = i ! current point becomes previous for next iteration
-      end do
+    end do
 
-      ! odd number of crossings = inside, even = outside
-      is_inside = (mod(crossings, 2) == 1)
+    ! odd number of crossings = inside, even = outside
+    is_inside = (mod(crossings, 2) == 1)
 
-      ! respect jins mode
-      if (jins == 0) then
-         is_inside = .not. is_inside
-      end if
+    ! respect jins mode
+    if (jins == 0) then
+        is_inside = .not. is_inside
+    end if
 
-   end function pinpok_raycast
+end function pinpok_raycast
    !
    ! PINPOK
    !
@@ -865,182 +868,6 @@ function matprod(A, b)
       matprod = (/(A(i, 1) * b(1) + A(i, 2) * b(2) + A(i, 3) * b(3), i=1, 3)/)
    end function matprod
 
-   subroutine dbpinpol(xp, yp, in, dmiss, JINS, NPL, xpl, ypl, zpl) ! ALS JE VOOR VEEL PUNTEN MOET NAGAAN OF ZE IN POLYGON ZITTEN
-      implicit none
-      real(kind=dp), intent(in) :: xp, yp !< point coordinates
-      integer, intent(inout) :: in !< in(-1): initialization, out(0): outside polygon, out(1): inside polygon
-      integer :: num
-      real(kind=dp), intent(in) :: dmiss
-      integer, intent(in) :: JINS, NPL
-      real(kind=dp), optional, intent(in) :: xpl(:), ypl(:), zpl(:)
-
-      if (NPL > 0 .and. present(xpl)) then
-         call dbpinpol_optinside_perpol(xp, yp, 0, 0, in, num, dmiss, JINS, NPL, xpl, ypl, zpl)
-      else
-         call dbpinpol_optinside_perpol(xp, yp, 0, 0, in, num, dmiss, JINS, NPL)
-      end if
-   end subroutine dbpinpol
-
-   !> The original dbpinpol routine, extended with an optional per-polygon-specified inside-mode.
-      !! Used this for checking for many points whether they are inside the global polygons.
-      !! Optionally, the global jins=1/other:inside/outside-mode can be replaced by an
-      !! inside/outside mode per polygon: that should then be stored as a +1/-1 in the first
-      !! zpl(istart) point of each polygon.
-   subroutine dbpinpol_optinside_perpol(xp, yp, inside_perpol, iselect, in, numselect, dmiss, JINS, NPL, xpl, ypl, zpl) ! ALS JE VOOR VEEL PUNTEN MOET NAGAAN OF ZE IN POLYGON ZITTEN
-
-      use m_alloc
-
-      implicit none
-
-      real(kind=dp), intent(in) :: xp, yp !< point coordinates
-      integer, intent(in) :: inside_perpol !< Specify whether or not (1/0) to use each polygon's first point zpl-value as the jins(ide)-option (only 0 or 1 allowed), or use the global JINS variable.
-      integer, intent(in) :: iselect !< use all polygons (0), only first-zpl<0 polygons (-1), or all but first-zpl<0 polygons (1)
-      integer, intent(inout) :: in !< in(-1): initialization, out(0): outside polygon, out(1): inside polygon
-      integer, intent(inout) :: numselect !< number of polygons of "iselect" type considered
-
-      integer :: MAXPOLY = 1000 ! will grow if needed
-
-      real(kind=dp), allocatable, save :: xpmin(:), ypmin(:), xpmax(:), ypmax(:)
-      integer, save :: Npoly
-      integer, allocatable, save :: iistart(:), iiend(:)
-
-      integer :: ipoint ! points to first part of a polygon-subsection in polygon array
-      integer :: istart, iend ! point to start and and node of a polygon in polygon array respectively
-      integer :: ipoly ! polygon number
-
-      logical :: Linit ! initialization of polygon bounds, and start and end nodes respectively
-
-      integer :: jins_opt !< The actual used jins-mode (either global, or per poly)
-      real(kind=dp), intent(in) :: dmiss
-      integer, intent(in) :: JINS, NPL
-      real(kind=dp), optional, intent(in) :: xpl(NPL), ypl(NPL), zpl(NPL)
-      integer :: count
-      numselect = 0
-
-      if (NPL == 0) then
-         in = 1
-         return
-      end if
-
-      Linit = (in < 0)
-
-      in = 0
-
-      !     initialization
-      if (Linit) then
-         !         write(6,"('dbpinpol: init... ', $)")
-         ipoint = 1
-         ipoly = 0
-         call realloc(xpmin, maxpoly, keepExisting=.false.)
-         call realloc(xpmax, maxpoly, keepExisting=.false.)
-         call realloc(ypmin, maxpoly, keepExisting=.false.)
-         call realloc(ypmax, maxpoly, keepExisting=.false.)
-         call realloc(iistart, maxpoly, keepExisting=.false.)
-         call realloc(iiend, maxpoly, keepExisting=.false.)
-
-         do while (ipoint < NPL)
-            ipoly = ipoly + 1
-            if (ipoly > maxpoly) then
-               maxpoly = ceiling(maxpoly * 1.1)
-               call realloc(xpmin, maxpoly, keepExisting=.true.)
-               call realloc(xpmax, maxpoly, keepExisting=.true.)
-               call realloc(ypmin, maxpoly, keepExisting=.true.)
-               call realloc(ypmax, maxpoly, keepExisting=.true.)
-               call realloc(iistart, maxpoly, keepExisting=.true.)
-               call realloc(iiend, maxpoly, keepExisting=.true.)
-            end if
-
-            !           get polygon start and end pointer respectively
-            call get_startend(NPL - ipoint + 1, xpl(ipoint:NPL), ypl(ipoint:NPL), istart, iend, dmiss)
-            istart = istart + ipoint - 1
-            iend = iend + ipoint - 1
-
-            if (istart >= iend .or. iend > NPL) exit ! done
-
-            xpmin(ipoly) = minval(xpl(istart:iend))
-            xpmax(ipoly) = maxval(xpl(istart:iend))
-            ypmin(ipoly) = minval(ypl(istart:iend))
-            ypmax(ipoly) = maxval(ypl(istart:iend))
-
-            iistart(ipoly) = istart
-            iiend(ipoly) = iend
-
-            !           advance pointer
-            ipoint = iend + 2
-         end do ! do while ( ipoint < NPL .and. ipoly < MAXPOLY )
-         Npoly = ipoly
-
-         !         write(6,"('done, Npoly=', I4)") Npoly
-      end if
-      count = 0
-      do ipoly = 1, Npoly
-         istart = iistart(ipoly)
-         iend = iiend(ipoly)
-
-         !         write(6,"('dbpinpol: ipoly=', I4, ', istart=', I16, ', iend=', I16)") ipoly, istart, iend
-
-         if (istart >= iend .or. iend > NPL) exit ! done
-
-         if (iselect == -1 .and. (zpl(istart) == DMISS .or. zpl(istart) >= 0)) cycle
-         if (iselect == 1 .and. (zpl(istart) /= DMISS .and. zpl(istart) < 0)) cycle
-
-         numselect = numselect + 1
-
-         if (inside_perpol == 1 .and. zpl(istart) /= dmiss) then ! only if third column was actually supplied
-            jins_opt = int(zpl(istart)) ! Use inside-option per each polygon.
-         else
-            jins_opt = JINS ! Use global inside-option.
-         end if
-
-         if (jins_opt == 1) then ! inside polygon
-            if (xp >= xpmin(ipoly) .and. xp <= xpmax(ipoly) .and. &
-                yp >= ypmin(ipoly) .and. yp <= ypmax(ipoly)) then
-               call PINPOK(Xp, Yp, iend - istart + 1, xpl(istart), ypl(istart), IN, jins, dmiss)
-               if (jins_opt > 0 .neqv. JINS > 0) then ! PINPOK has used global jins, but polygon asked the exact opposite, so negate the result here.
-                  IN = 1 - in ! IN-1
-               end if
-
-               if (in == 1) then
-                  count = count + 1
-               end if
-            end if
-         else ! outside polygon
-            if (xp >= xpmin(ipoly) .and. xp <= xpmax(ipoly) .and. &
-                yp >= ypmin(ipoly) .and. yp <= ypmax(ipoly)) then
-               call PINPOK(Xp, Yp, iend - istart + 1, xpl(istart), ypl(istart), IN, jins, dmiss)
-               if (jins_opt > 0 .neqv. JINS > 0) then ! PINPOK has used global jins, but polygon asked the exact opposite, so negate the result here.
-                  IN = 1 - in ! IN-1
-               end if
-
-               if (in == 1) then
-                  continue
-               else
-                  count = count + 1
-               end if
-               !else
-               !   in = 1 ! outside check succeeded (completely outside of polygon's bounding box), return 'true'.
-               !   exit
-            end if
-         end if
-      end do ! do ipoly=1,Npoly
-
-      if (jins == 1) then !< NB: here, it is checked based on jins, while in the loop based on jins_opt, need to be fixed
-         if (mod(count, 2) == 1) then
-            in = 1
-         else
-            in = 0
-         end if
-      else
-         if (mod(count, 2) == 1) then
-            in = 0
-         else
-            in = 1
-         end if
-      end if
-
-      return
-   end subroutine dbpinpol_optinside_perpol
-
    !>  get the start and end index of the first enclosed non-DMISS subarray
    subroutine get_startend(num, x, y, jstart, jend, dmiss)
 
diff --git a/src/utils_lgpl/deltares_common/packages/deltares_common/src/geometry_module_dbpinpol.f90 b/src/utils_lgpl/deltares_common/packages/deltares_common/src/geometry_module_dbpinpol.f90
new file mode 100644
index 0000000000..a1abc2854d
--- /dev/null
+++ b/src/utils_lgpl/deltares_common/packages/deltares_common/src/geometry_module_dbpinpol.f90
@@ -0,0 +1,83 @@
+!----- LGPL --------------------------------------------------------------------
+!
+!  Copyright (C)  Stichting Deltares, 2011-2026.
+!
+!  This library is free software; you can redistribute it and/or
+!  modify it under the terms of the GNU Lesser General Public
+!  License as published by the Free Software Foundation version 2.1.
+!
+!  This library is distributed in the hope that it will be useful,
+!  but WITHOUT ANY WARRANTY; without even the implied warranty of
+!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+!  Lesser General Public License for more details.
+!
+!  You should have received a copy of the GNU Lesser General Public
+!  License along with this library; if not, see <http://www.gnu.org/licenses/>.
+!
+!  contact: delft3d.support@deltares.nl
+!  Stichting Deltares
+!  P.O. Box 177
+!  2600 MH Delft, The Netherlands
+!
+!  All indications and logos of, and references to, "Delft3D" and "Deltares"
+!  are registered trademarks of Stichting Deltares, and remain the property of
+!  Stichting Deltares. All rights reserved.
+!
+!-------------------------------------------------------------------------------
+!
+submodule (geometry_module) geometry_module_dbpinpol
+   
+   implicit none
+   
+contains
+   
+   module subroutine dbpinpol(xp, yp, in, dmiss, JINS, NPL, xpl, ypl, zpl)
+      use m_cellmask_from_polygon_set, only: cellmask_from_polygon_set_init, &
+                                             cellmask_from_polygon_set, &
+                                             cellmask_from_polygon_set_cleanup
+      implicit none
+
+      real(kind=dp), intent(in) :: xp, yp
+      integer, intent(inout) :: in
+      real(kind=dp), intent(in) :: dmiss
+      integer, intent(in) :: JINS, NPL
+      real(kind=dp), optional, intent(in) :: xpl(:), ypl(:), zpl(:)
+
+      integer :: num
+      logical, save :: initialized = .false.
+
+      ! Special case: NPL == 0 means "no polygons to check, everything is considered inside"
+      if (NPL == 0) then
+         in = 1
+         return
+      end if
+
+      ! Initialization phase (when in < 0)
+      if (in < 0) then
+         ! Clean up any previous initialization
+         if (initialized) then
+            call cellmask_from_polygon_set_cleanup()
+         end if
+
+         ! Build optimized spatial index
+         if (present(xpl)) then
+            call cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl)
+            initialized = .true.
+         end if
+
+         in = 0 ! Reset for subsequent queries
+      end if
+
+      ! Query phase (when in >= 0)
+      if (.not. initialized) then
+         ! Safety: if someone forgot to initialize
+         in = 0
+         return
+      end if
+
+      ! Use your optimized point-in-polygon with bounding boxes
+      in = cellmask_from_polygon_set(xp, yp)
+
+   end subroutine dbpinpol
+
+end submodule geometry_module_dbpinpol
\ No newline at end of file
diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 b/src/utils_lgpl/deltares_common/packages/deltares_common/src/m_cellmask_from_polygon_set.f90
similarity index 82%
rename from src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90
rename to src/utils_lgpl/deltares_common/packages/deltares_common/src/m_cellmask_from_polygon_set.f90
index 35c5bce6ff..c13ddf055d 100644
--- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90
+++ b/src/utils_lgpl/deltares_common/packages/deltares_common/src/m_cellmask_from_polygon_set.f90
@@ -32,13 +32,14 @@ module m_cellmask_from_polygon_set
    use precision, only: dp
    use m_polygon, only: xpl, ypl, zpl, npl, maxpol, restorepol, savepol
 
-   implicit none
+   implicit none(external)
 
    private
 
    public :: cellmask_from_polygon_set_init, cellmask_from_polygon_set_cleanup, cellmask_from_polygon_set, pinpok_elemental
-   public :: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
-   public :: find_cells_crossed_by_polyline
+   public :: point_find_netcell, cleanup_cell_geom_polylines
+   public :: cellmask_initialized
+   public :: find_cells_for_segment, line_segments_intersect, point_to_line_distance
 
    integer :: polygons = 0 !< Number of polygons stored in module arrays xpl, ypl, zpl
    real(kind=dp), allocatable :: x_poly_min(:), y_poly_min(:) !< Polygon bounding box min coordinates, (dim = polygons)
@@ -242,57 +243,6 @@ elemental function pinpok_elemental(x, y, i_poly) result(is_inside)
 
    end function pinpok_elemental
 
-   !> Initialize xpl, ypl, zpl arrays with all netcell geometries (called once)
-   subroutine init_cell_geom_as_polylines()
-      use network_data
-      use m_alloc
-
-      integer :: k, n, k1, total_points, ipoint
-
-      if (cellmask_initialized) then !> reuse cellmask cache boolean
-         return
-      end if
-
-      call savepol()
-
-      ! calculate total points needed: sum(netcell(k)%n + 1) for all cells
-      ! +1 for dmiss separator after each polygon
-      total_points = 0
-      do k = 1, nump
-         total_points = total_points + netcell(k)%n + 1 ! +1 for dmiss
-      end do
-
-      ! allocate or reallocate xpl, ypl, zpl
-      call realloc(xpl, total_points, keepexisting=.false.)
-      call realloc(ypl, total_points, keepexisting=.false.)
-      call realloc(zpl, total_points, keepexisting=.false.)
-
-      ! fill arrays with netcell geometry
-      ipoint = 0
-      do k = 1, nump
-         do n = 1, netcell(k)%n
-            ipoint = ipoint + 1
-            k1 = netcell(k)%nod(n)
-            xpl(ipoint) = xk(k1)
-            ypl(ipoint) = yk(k1)
-            zpl(ipoint) = real(k, dp) ! store cell index as z-value
-         end do
-
-         ! add separator
-         ipoint = ipoint + 1
-         xpl(ipoint) = dmiss
-         ypl(ipoint) = dmiss
-         zpl(ipoint) = dmiss
-      end do
-
-      npl = ipoint
-
-      ! initialize the cellmask module with these polygons
-      ! this builds bounding boxes and polygon indices
-      call cellmask_from_polygon_set_init(npl, xpl, ypl, zpl)
-
-   end subroutine init_cell_geom_as_polylines
-
    !> call general polygon cleanup and restore previous polygon data
    subroutine cleanup_cell_geom_polylines()
       call cellmask_from_polygon_set_cleanup()
@@ -333,50 +283,10 @@ elemental function point_find_netcell(x, y) result(k)
 
    end function point_find_netcell
 
-!> Find all cells crossed by polyline using brute force on cached geometry. The routine is inclusive of edge cases (touching edges or vertices).
-   subroutine find_cells_crossed_by_polyline(xpoly, ypoly, crossed_cells, error)
-      use m_alloc, only: realloc
-      use network_data, only: cellmask, nump
-      use m_missing, only: dmiss
-
-      implicit none
-
-      real(kind=dp), intent(in) :: xpoly(:) !< Polyline x-coordinates
-      real(kind=dp), intent(in) :: ypoly(:) !< Polyline y-coordinates
-      integer, allocatable, intent(out) :: crossed_cells(:) !> Indices of crossed cells in network_data::netcells
-      character, dimension(:), allocatable, intent(out) :: error !> Error message, empty if no error, to be handled at call site
-
-      integer :: npoly, i
-
-      error = ''
-
-      npoly = size(xpoly)
-      if (any(xpoly == dmiss) .or. any(ypoly == dmiss) .or. npoly < 2) then
-         error = 'Invalid polyline input'
-         return
-      end if
-
-      if (.not. cellmask_initialized) then
-         call init_cell_geom_as_polylines()
-      end if
-
-      call realloc(cellmask, nump, keepexisting=.false., fill=0)
-
-      ! Process each segment and put the result in cellmask
-      do i = 1, npoly - 1
-         call find_cells_for_segment(xpoly(i), ypoly(i), xpoly(i + 1), ypoly(i + 1), cellmask)
-      end do
-
-      crossed_cells = pack([(i, i=1, nump)], mask=(cellmask == 1))
-
-   end subroutine find_cells_crossed_by_polyline
-
 !> Find all cells that a segment crosses and mark them in cellmask
    subroutine find_cells_for_segment(xa, ya, xb, yb, cellmask)
       use m_polygon, only: xpl, ypl
 
-      implicit none
-
       real(kind=dp), intent(in) :: xa, ya, xb, yb !< Segment endpoints
       integer, intent(inout) :: cellmask(:) !< Cell mask array: 1=crossed, 0=not crossed
 
diff --git a/src/utils_lgpl/deltares_common/packages/deltares_common/src/m_polygon.f90 b/src/utils_lgpl/deltares_common/packages/deltares_common/src/m_polygon.f90
new file mode 100644
index 0000000000..bbb7533a7c
--- /dev/null
+++ b/src/utils_lgpl/deltares_common/packages/deltares_common/src/m_polygon.f90
@@ -0,0 +1,155 @@
+module m_polygon
+
+   implicit none
+
+   double precision, allocatable :: XPL(:), YPL(:), ZPL(:), DZL(:), DZR(:), DCREST(:), DTL(:), DTR(:), DVEG(:)
+   double precision, allocatable, private :: XPH(:), YPH(:), ZPH(:)
+   integer, allocatable :: IWEIRT(:)
+   integer :: NPL, NPH, MAXPOL, MP, MPS, jakol45 = 0
+   character(len=64), allocatable :: nampli(:) ! Names of polylines, set in reapol,
+   ! not shifted/updated during editpol.
+   double precision :: dxuni = 40d0 ! uniform spacing
+   integer :: MAXPOLY = 1000 ! will grow if needed
+   double precision, allocatable :: xpmin(:), ypmin(:), xpmax(:), ypmax(:), zpmin(:), zpmax(:)
+   integer :: Npoly
+   integer, allocatable :: iistart(:), iiend(:)
+   integer, allocatable :: ipsection(:)
+
+contains
+   !> Increase size of global polyline array.
+      !! Specify new size and whether existing points need to be maintained.
+   subroutine increasepol(N, jaKeepExisting)
+      use m_missing
+      use m_alloc
+      implicit none
+      integer :: n !< Desired new minimum size
+      integer :: jaKeepExisting !< Whether or not (1/0) to keep existing points.
+      logical :: jakeep
+      integer :: maxpolcur
+      integer :: ierr
+
+      maxpolcur = size(xpl)
+      if (N <= maxpolcur) then
+         return
+      end if
+      MAXPOL = max(100000, int(5d0 * N))
+
+      jakeep = jaKeepExisting == 1
+
+      call realloc(xpl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
+      call realloc(ypl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
+      call realloc(zpl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
+
+      if (jakol45 == 1) then
+         call realloc(dzl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
+         call realloc(dzr, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
+      else if (jakol45 == 2) then
+         call realloc(dcrest, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
+         call realloc(dzl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
+         call realloc(dzr, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
+         call realloc(dtl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
+         call realloc(dtr, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
+         call realloc(dveg, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
+         call realloc(iweirt, maxpol, keepExisting=jakeep, stat=ierr)
+      end if
+
+      !     make sure nampli is allocated
+      if (.not. allocated(nampli)) then
+         allocate (nampli(0))
+      end if
+
+   end subroutine increasepol
+
+   !> Copies the global polygon into the backup polygon arrays.
+   subroutine SAVEPOL()
+
+      use m_alloc
+      use m_missing
+      implicit none
+
+      if (NPL > 0) then
+         call realloc(xph, maxpol, keepExisting=.false.)
+         call realloc(yph, maxpol, keepExisting=.false.)
+         call realloc(zph, maxpol, keepExisting=.false.)
+         XPH(1:NPL) = XPL(1:NPL)
+         YPH(1:NPL) = YPL(1:NPL)
+         ZPH(1:NPL) = ZPL(1:NPL)
+      end if
+
+      MPS = MP
+      NPH = NPL
+
+      return
+   end subroutine savepol
+
+   !> Puts back a previously saved backup polygon into the global polygon arrays.
+   subroutine RESTOREPOL()
+      use m_alloc
+      use m_missing
+      implicit none
+
+      maxpol = max(maxpol, nph)
+      call realloc(xpl, maxpol, keepExisting=.false.)
+      call realloc(ypl, maxpol, keepExisting=.false.)
+      call realloc(zpl, maxpol, keepExisting=.false.)
+
+      if (NPH > 0) then
+         XPL(1:NPH) = XPH(1:NPH)
+         YPL(1:NPH) = YPH(1:NPH)
+         ZPL(1:NPH) = ZPH(1:NPH)
+         deallocate (XPH, YPH, ZPH)
+      end if
+      MP = MPS
+      NPL = NPH
+
+      return
+   end subroutine restorepol
+
+   function m_polygon_destructor() result(ierr)
+
+      implicit none
+
+      integer :: ierr
+
+      ierr = 0
+
+      if (allocated(XPL) .and. ierr == 0) deallocate (XPL, stat=ierr)
+      if (allocated(YPL) .and. ierr == 0) deallocate (YPL, stat=ierr)
+      if (allocated(ZPL) .and. ierr == 0) deallocate (ZPL, stat=ierr)
+      if (allocated(XPH) .and. ierr == 0) deallocate (XPH, stat=ierr)
+
+      if (allocated(YPH) .and. ierr == 0) deallocate (YPH, stat=ierr)
+      if (allocated(ZPH) .and. ierr == 0) deallocate (ZPH, stat=ierr)
+      if (allocated(ZPH) .and. ierr == 0) deallocate (ZPH, stat=ierr)
+
+      if (allocated(DZL) .and. ierr == 0) deallocate (DZL, stat=ierr)
+      if (allocated(DZR) .and. ierr == 0) deallocate (DZR, stat=ierr)
+      if (allocated(DCREST) .and. ierr == 0) deallocate (DCREST, stat=ierr)
+      if (allocated(DTL) .and. ierr == 0) deallocate (DTL, stat=ierr)
+      if (allocated(DTR) .and. ierr == 0) deallocate (DTR, stat=ierr)
+      if (allocated(DVEG) .and. ierr == 0) deallocate (DVEG, stat=ierr)
+      if (allocated(IWEIRT) .and. ierr == 0) deallocate (IWEIRT, stat=ierr)
+
+      if (allocated(xpmin) .and. ierr == 0) deallocate (xpmin, stat=ierr)
+      if (allocated(ypmin) .and. ierr == 0) deallocate (ypmin, stat=ierr)
+      if (allocated(xpmax) .and. ierr == 0) deallocate (xpmax, stat=ierr)
+      if (allocated(ypmax) .and. ierr == 0) deallocate (ypmax, stat=ierr)
+      if (allocated(zpmin) .and. ierr == 0) deallocate (zpmin, stat=ierr)
+      if (allocated(zpmax) .and. ierr == 0) deallocate (zpmax, stat=ierr)
+      if (allocated(iistart) .and. ierr == 0) deallocate (iistart, stat=ierr)
+      if (allocated(iiend) .and. ierr == 0) deallocate (iiend, stat=ierr)
+      if (allocated(ipsection) .and. ierr == 0) deallocate (ipsection, stat=ierr)
+
+      jakol45 = 0
+      dxuni = 40d0
+      MAXPOLY = 1000
+      NPL = 0
+      NPH = 0
+      MAXPOL = 0
+      MP = 0
+      MPS = 0
+      Npoly = 0
+
+   end function m_polygon_destructor
+
+end module m_polygon
diff --git a/src/utils_lgpl/gridgeom/packages/gridgeom/src/generalmodules.f90 b/src/utils_lgpl/gridgeom/packages/gridgeom/src/generalmodules.f90
index 0b8eaea174..69d6da4e31 100644
--- a/src/utils_lgpl/gridgeom/packages/gridgeom/src/generalmodules.f90
+++ b/src/utils_lgpl/gridgeom/packages/gridgeom/src/generalmodules.f90
@@ -99,185 +99,6 @@ subroutine default_sferic()
    end subroutine default_sferic
 end module m_sferic
 
-module m_polygon
-   use precision, only: dp
-
-   implicit none
-
-   real(kind=dp), allocatable, dimension(:) :: xpl
-   real(kind=dp), allocatable, dimension(:) :: ypl
-   real(kind=dp), allocatable, dimension(:) :: zpl
-   real(kind=dp), allocatable, dimension(:) :: dzl
-   real(kind=dp), allocatable, dimension(:) :: dzr
-   real(kind=dp), allocatable, dimension(:) :: dcrest
-   real(kind=dp), allocatable, dimension(:) :: dtl
-   real(kind=dp), allocatable, dimension(:) :: dtr
-   real(kind=dp), allocatable, dimension(:) :: dveg
-   real(kind=dp), allocatable, dimension(:), private :: xph
-   real(kind=dp), allocatable, dimension(:), private :: yph
-   real(kind=dp), allocatable, dimension(:), private :: zph
-   integer, allocatable, dimension(:) :: iweirt
-   integer :: npl
-   integer :: nph
-   integer :: maxpol
-   integer :: mp
-   integer :: mps
-   integer :: jakol45 = 0
-   character(len=64), allocatable, dimension(:) :: nampli ! Names of polylines, set in reapol,
-
-   ! not shifted/updated during editpol.
-   real(kind=dp) :: dxuni = 40d0 ! uniform spacing
-   integer :: maxpoly = 1000 ! will grow if needed
-   real(kind=dp), allocatable, dimension(:) :: xpmin
-   real(kind=dp), allocatable, dimension(:) :: ypmin
-   real(kind=dp), allocatable, dimension(:) :: xpmax
-   real(kind=dp), allocatable, dimension(:) :: ypmax
-   real(kind=dp), allocatable, dimension(:) :: zpmin
-   real(kind=dp), allocatable, dimension(:) :: zpmax
-   integer :: Npoly
-   integer, allocatable, dimension(:) :: iistart
-   integer, allocatable, dimension(:) :: iiend
-   integer, allocatable, dimension(:) :: ipsection
-
-contains
-   !> Increase size of global polyline array.
-      !! Specify new size and whether existing points need to be maintained.
-   subroutine increasepol(N, jaKeepExisting)
-      use m_missing
-      use m_alloc
-      implicit none
-      integer :: n !< Desired new minimum size
-      integer :: jaKeepExisting !< Whether or not (1/0) to keep existing points.
-      logical :: jakeep
-      integer :: maxpolcur
-      integer :: ierr
-
-      maxpolcur = size(xpl)
-      if (N <= maxpolcur) then
-         return
-      end if
-      MAXPOL = max(100000, int(5d0 * N))
-
-      jakeep = jaKeepExisting == 1
-
-      call realloc(xpl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
-      call realloc(ypl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
-      call realloc(zpl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
-
-      if (jakol45 == 1) then
-         call realloc(dzl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
-         call realloc(dzr, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
-      else if (jakol45 == 2) then
-         call realloc(dcrest, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
-         call realloc(dzl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
-         call realloc(dzr, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
-         call realloc(dtl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
-         call realloc(dtr, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
-         call realloc(dveg, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
-         call realloc(iweirt, maxpol, keepExisting=jakeep, stat=ierr)
-      end if
-
-      !     make sure nampli is allocated
-      if (.not. allocated(nampli)) then
-         allocate (nampli(0))
-      end if
-
-   end subroutine increasepol
-
-   !> Copies the global polygon into the backup polygon arrays.
-   subroutine SAVEPOL()
-
-      use m_alloc
-      use m_missing
-      implicit none
-
-      if (NPL > 0) then
-         call realloc(xph, maxpol, keepExisting=.false.)
-         call realloc(yph, maxpol, keepExisting=.false.)
-         call realloc(zph, maxpol, keepExisting=.false.)
-         XPH(1:NPL) = XPL(1:NPL)
-         YPH(1:NPL) = YPL(1:NPL)
-         ZPH(1:NPL) = ZPL(1:NPL)
-      end if
-
-      MPS = MP
-      NPH = NPL
-
-      return
-   end subroutine savepol
-
-   !> Puts back a previously saved backup polygon into the global polygon arrays.
-   subroutine RESTOREPOL()
-      use m_alloc
-      use m_missing
-      implicit none
-
-      maxpol = max(maxpol, nph)
-      call realloc(xpl, maxpol, keepExisting=.false.)
-      call realloc(ypl, maxpol, keepExisting=.false.)
-      call realloc(zpl, maxpol, keepExisting=.false.)
-
-      if (NPH > 0) then
-         XPL(1:NPH) = XPH(1:NPH)
-         YPL(1:NPH) = YPH(1:NPH)
-         ZPL(1:NPH) = ZPH(1:NPH)
-         deallocate(XPH, YPH, ZPH)
-      end if
-      MP = MPS
-      NPL = NPH
-
-      return
-   end subroutine restorepol
-
-   function m_polygon_destructor() result(ierr)
-
-      implicit none
-
-      integer :: ierr
-
-      ierr = 0
-
-      if (allocated(XPL) .and. ierr == 0) deallocate (XPL, stat=ierr)
-      if (allocated(YPL) .and. ierr == 0) deallocate (YPL, stat=ierr)
-      if (allocated(ZPL) .and. ierr == 0) deallocate (ZPL, stat=ierr)
-      if (allocated(XPH) .and. ierr == 0) deallocate (XPH, stat=ierr)
-
-      if (allocated(YPH) .and. ierr == 0) deallocate (YPH, stat=ierr)
-      if (allocated(ZPH) .and. ierr == 0) deallocate (ZPH, stat=ierr)
-      if (allocated(ZPH) .and. ierr == 0) deallocate (ZPH, stat=ierr)
-
-      if (allocated(DZL) .and. ierr == 0) deallocate (DZL, stat=ierr)
-      if (allocated(DZR) .and. ierr == 0) deallocate (DZR, stat=ierr)
-      if (allocated(DCREST) .and. ierr == 0) deallocate (DCREST, stat=ierr)
-      if (allocated(DTL) .and. ierr == 0) deallocate (DTL, stat=ierr)
-      if (allocated(DTR) .and. ierr == 0) deallocate (DTR, stat=ierr)
-      if (allocated(DVEG) .and. ierr == 0) deallocate (DVEG, stat=ierr)
-      if (allocated(IWEIRT) .and. ierr == 0) deallocate (IWEIRT, stat=ierr)
-
-      if (allocated(xpmin) .and. ierr == 0) deallocate (xpmin, stat=ierr)
-      if (allocated(ypmin) .and. ierr == 0) deallocate (ypmin, stat=ierr)
-      if (allocated(xpmax) .and. ierr == 0) deallocate (xpmax, stat=ierr)
-      if (allocated(ypmax) .and. ierr == 0) deallocate (ypmax, stat=ierr)
-      if (allocated(zpmin) .and. ierr == 0) deallocate (zpmin, stat=ierr)
-      if (allocated(zpmax) .and. ierr == 0) deallocate (zpmax, stat=ierr)
-      if (allocated(iistart) .and. ierr == 0) deallocate (iistart, stat=ierr)
-      if (allocated(iiend) .and. ierr == 0) deallocate (iiend, stat=ierr)
-      if (allocated(ipsection) .and. ierr == 0) deallocate (ipsection, stat=ierr)
-
-      jakol45 = 0
-      dxuni = 40d0
-      MAXPOLY = 1000
-      NPL = 0
-      NPH = 0
-      MAXPOL = 0
-      MP = 0
-      MPS = 0
-      Npoly = 0
-
-   end function m_polygon_destructor
-
-end module m_polygon
-
 !
 ! Stores the coordinates of the cells
 !