diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridOrthogonalization.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridOrthogonalization.hpp
index 274166081..97406b620 100644
--- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridOrthogonalization.hpp
+++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridOrthogonalization.hpp
@@ -82,11 +82,17 @@ namespace meshkernel
         /// @return A matrix whose coefficients are true if the node is an invalid boundary node, false otherwise
         [[nodiscard]] lin_alg::Matrix<bool> ComputeInvalidVerticalBoundaryNodes() const;
 
-        /// @brief Compute and assign grid points for grid line section
-        void ComputePointsForGridLine(const UInt m,
-                                      const UInt n,
-                                      const UInt startN,
-                                      const int nextVertical);
+        /// @brief Compute and assign grid points for grid line section in n-direction
+        void ComputePointsForGridLineN(const UInt m,
+                                       const UInt n,
+                                       const UInt startN,
+                                       const int nextVertical);
+
+        /// @brief Compute and assign grid points for grid line section in m-direction
+        void ComputePointsForGridLineM(const UInt m,
+                                       const UInt n,
+                                       const UInt startM,
+                                       const int nextHorizontal);
 
         OrthogonalizationParameters m_orthogonalizationParameters; ///< The orthogonalization parameters
 
diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridOrthogonalization.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridOrthogonalization.cpp
index 800f8a6aa..9d83b9894 100644
--- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridOrthogonalization.cpp
+++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridOrthogonalization.cpp
@@ -97,10 +97,10 @@ meshkernel::UndoActionPtr CurvilinearGridOrthogonalization::Compute()
     return undoAction;
 }
 
-void CurvilinearGridOrthogonalization::ComputePointsForGridLine(const UInt m,
-                                                                const UInt n,
-                                                                const UInt startN,
-                                                                const int nextVertical)
+void CurvilinearGridOrthogonalization::ComputePointsForGridLineN(const UInt m,
+                                                                 const UInt n,
+                                                                 const UInt startN,
+                                                                 const int nextVertical)
 {
 
     for (auto nn = startN + 1; nn < n; ++nn)
@@ -150,6 +150,60 @@ void CurvilinearGridOrthogonalization::ComputePointsForGridLine(const UInt m,
     }
 }
 
+void CurvilinearGridOrthogonalization::ComputePointsForGridLineM(const UInt m,
+                                                                 const UInt n,
+                                                                 const UInt startM,
+                                                                 const int nextHorizontal)
+{
+
+    for (auto mm = startM + 1; mm < m; ++mm)
+    {
+
+        if (mm < m_lowerLeft.m_m || mm > m_upperRight.m_m || n < m_lowerLeft.m_n || n > m_upperRight.m_n)
+        {
+            continue;
+        }
+        if (m_grid.GetNodeType(n, m) == NodeType::Invalid)
+        {
+            continue;
+        }
+        const auto bottomNode = m_grid.GetNode(n, mm - 1);
+        const auto horizontalNode = m_grid.GetNode(n + nextHorizontal, mm);
+        const auto upperNode = m_grid.GetNode(n, mm + 1);
+
+        Point boundaryNode;
+        if (nextHorizontal == 1)
+        {
+            const auto qb = 1.0 / m_orthoEqTerms.atp(n, mm - 1);
+            const auto qc = 1.0 / m_orthoEqTerms.atp(n, mm);
+            const auto qbc = m_orthoEqTerms.atp(n, mm - 1) + m_orthoEqTerms.atp(n, mm);
+            const auto rn = qb + qc + qbc;
+            boundaryNode.x = (bottomNode.x * qb + horizontalNode.x * qbc + upperNode.x * qc + bottomNode.y - upperNode.y) / rn;
+            boundaryNode.y = (bottomNode.y * qb + horizontalNode.y * qbc + upperNode.y * qc + upperNode.x - bottomNode.x) / rn;
+        }
+
+        if (nextHorizontal == -1)
+        {
+            const auto qb = 1.0 / m_orthoEqTerms.atp(n - 1, mm - 1);
+            const auto qc = 1.0 / m_orthoEqTerms.atp(n - 1, mm);
+            const auto qbc = m_orthoEqTerms.atp(n - 1, mm - 1) + m_orthoEqTerms.atp(n - 1, mm);
+            const auto rn = qb + qc + qbc;
+            boundaryNode.x = (bottomNode.x * qb + horizontalNode.x * qbc + upperNode.x * qc + upperNode.y - bottomNode.y) / rn;
+            boundaryNode.y = (bottomNode.y * qb + horizontalNode.y * qbc + upperNode.y * qc + bottomNode.x - upperNode.x) / rn;
+        }
+
+        // Vertical spline index
+        const auto splineIndex = m_grid.NumM() + n;
+
+        const auto node_value = m_splines.ComputeClosestPointOnSplineSegment(splineIndex,
+                                                                             startM,
+                                                                             m,
+                                                                             boundaryNode);
+
+        m_grid.GetNode(n, mm) = node_value;
+    }
+}
+
 void CurvilinearGridOrthogonalization::ProjectHorizontalBoundaryGridNodes()
 {
     // m grid lines (horizontal)
@@ -182,7 +236,7 @@ void CurvilinearGridOrthogonalization::ProjectHorizontalBoundaryGridNodes()
                 (nodeType == NodeType::BottomRight || nodeType == NodeType::UpperRight) &&
                 nextVertical != 0)
             {
-                ComputePointsForGridLine(m, n, startN, nextVertical);
+                ComputePointsForGridLineN(m, n, startN, nextVertical);
             }
         }
     }
@@ -220,52 +274,7 @@ void CurvilinearGridOrthogonalization::ProjectVerticalBoundariesGridNodes()
                 nextHorizontal != 0 &&
                 startM != constants::missing::uintValue)
             {
-                for (auto mm = startM + 1; mm < m; ++mm)
-                {
-
-                    if (mm < m_lowerLeft.m_m || mm > m_upperRight.m_m || n < m_lowerLeft.m_n || n > m_upperRight.m_n)
-                    {
-                        continue;
-                    }
-                    if (m_grid.GetNodeType(n, m) == NodeType::Invalid)
-                    {
-                        continue;
-                    }
-                    const auto bottomNode = m_grid.GetNode(n, mm - 1);
-                    const auto horizontalNode = m_grid.GetNode(n + nextHorizontal, mm);
-                    const auto upperNode = m_grid.GetNode(n, mm + 1);
-
-                    Point boundaryNode;
-                    if (nextHorizontal == 1)
-                    {
-                        const auto qb = 1.0 / m_orthoEqTerms.atp(n, mm - 1);
-                        const auto qc = 1.0 / m_orthoEqTerms.atp(n, mm);
-                        const auto qbc = m_orthoEqTerms.atp(n, mm - 1) + m_orthoEqTerms.atp(n, mm);
-                        const auto rn = qb + qc + qbc;
-                        boundaryNode.x = (bottomNode.x * qb + horizontalNode.x * qbc + upperNode.x * qc + bottomNode.y - upperNode.y) / rn;
-                        boundaryNode.y = (bottomNode.y * qb + horizontalNode.y * qbc + upperNode.y * qc + upperNode.x - bottomNode.x) / rn;
-                    }
-
-                    if (nextHorizontal == -1)
-                    {
-                        const auto qb = 1.0 / m_orthoEqTerms.atp(n - 1, mm - 1);
-                        const auto qc = 1.0 / m_orthoEqTerms.atp(n - 1, mm);
-                        const auto qbc = m_orthoEqTerms.atp(n - 1, mm - 1) + m_orthoEqTerms.atp(n - 1, mm);
-                        const auto rn = qb + qc + qbc;
-                        boundaryNode.x = (bottomNode.x * qb + horizontalNode.x * qbc + upperNode.x * qc + upperNode.y - bottomNode.y) / rn;
-                        boundaryNode.y = (bottomNode.y * qb + horizontalNode.y * qbc + upperNode.y * qc + bottomNode.x - upperNode.x) / rn;
-                    }
-
-                    // Vertical spline index
-                    const auto splineIndex = m_grid.NumM() + n;
-
-                    const auto node_value = m_splines.ComputeClosestPointOnSplineSegment(splineIndex,
-                                                                                         startM,
-                                                                                         m,
-                                                                                         boundaryNode);
-
-                    m_grid.GetNode(n, mm) = node_value;
-                }
+                ComputePointsForGridLineM(m, n, startM, nextHorizontal);
             }
         }
     }