SOFTSUP-39123: Incorrect waterlevels

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/setbobs.f90

@@ -30,12 +30,12 @@
 !
 !
 module m_set_bobs
-   use m_get2dnormal, only: get2dnormal
-   use m_get1ddir, only: get1ddir
    use m_duikerstoprofs
 
-   implicit none
+   implicit none(type, external)
+
 contains
+
    subroutine setbobs() ! and set blu, weigthed depth at u point
       use precision, only: dp
       use m_netw, only: netcell, zk, zkuni, numl, numl1d, lne, kn, nmk
@@ -53,8 +53,7 @@ subroutine setbobs() ! and set blu, weigthed depth at u point
       use unstruc_model, only: md_convertlongculverts
 
       integer :: L, k1, k2, n1, n2, n, k, k3, LL, kk, Ls, Lf, mis, i, numcoords, ibotL
-      real(kind=dp) :: bl1, bl2, bedlevel_at_link, bln, zn1, zn2, zn3, wn, alf, skewn, xt, yt, xn, yn
-      ! real(kind=dp), external :: skewav
+      real(kind=dp) :: bl1, bl2, bedlevel_at_link, bln, zn1, zn2, zn3, wn, alf
 
       ! First, prepare bed levels at pressure points:
 
@@ -240,23 +239,10 @@ subroutine setbobs() ! and set blu, weigthed depth at u point
          end if
 
          if (kcu(L) == 3) then ! 1D2D internal link, bobs at minimum
-            if (kcs(n1) == 21) then
-               bedlevel_at_link = bl(n1)
-               call get2Dnormal(n1, xn, yn) ! xn, yn = 2D land normal vector pointing upward, both zero = flat
-               call get1Ddir(n2, xt, yt) ! xt, yt = 1D river tangential normal vector
-            end if
-            if (kcs(n2) == 21) then
-               bedlevel_at_link = bl(n2)
-               call get2Dnormal(n2, xn, yn)
-               call get1Ddir(n1, xt, yt)
-            end if
-            skewn = abs(xn * xt + yn * yt)
-            bob(1, L) = bedlevel_at_link
-            bob(2, L) = bedlevel_at_link ! revisit later+ wu(L)*skewn ! TODO: HK: why wu here? Why not dx(L) or something similar?
-            bob0(1, L) = bedlevel_at_link
-            bob0(2, L) = bedlevel_at_link ! revisit later+ wu(L)*skewn ! TODO: HK: why wu here? Why not dx(L) or something similar?
-            bl(n1) = min(bl(n1), bedlevel_at_link)
-            bl(n2) = min(bl(n2), bedlevel_at_link)
+            bob(1, L) = max(bl(n1), bl(n2))
+            bob(2, L) = max(bl(n1), bl(n2)) ! revisit later+ wu(L)*skewn ! TODO: HK: why wu here? Why not dx(L) or something similar?
+            bob0(1, L) = max(bl(n1), bl(n2))
+            bob0(2, L) = max(bl(n1), bl(n2)) ! revisit later+ wu(L)*skewn ! TODO: HK: why wu here? Why not dx(L) or something similar?
          else if (kcu(L) == BEDLEV_TYPE_MIN) then ! left right
             bedlevel_at_link = min(zn1, zn2)
             bob(1, L) = zn1
@@ -338,9 +324,8 @@ subroutine setbobs() ! and set blu, weigthed depth at u point
                bl(n1) = bl1
                bob(1, L) = max(bl(n1), bl(n2))
                bob(2, L) = bob(1, L)
-               !elseif (bl(n1) == 1d30 .or. bl(n2) == 30) then
             else if (.not. network%loaded) then
-!          SPvdP: previous expression is problematic when zk(k2) and/or zk(k3) have missing values
+!               SPvdP: previous expression is problematic when zk(k2) and/or zk(k3) have missing values
                zn2 = zk(k2)
                if (zn2 == dmiss) then
                   zn2 = zkuni
@@ -411,6 +396,7 @@ subroutine setbobs() ! and set blu, weigthed depth at u point
       jaupdbndbl = 0 ! after first run of setbobs set to 0 = no update
 
    end subroutine setbobs
+
    !> calculate bed level at a link based on the type of bed level definition and the input parameters.
    function get_bedlevel_at_link(n1, n2, k1, k2, blu, ibot) result(bedlevel_at_link)