unst 8559: discharge boundaries

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_data/fm_external_forcings_data.f90

@@ -87,6 +87,7 @@ module fm_external_forcings_data
                                                         !! 1,*=boundary typ, see type indicator kbndz(4,*) below
    integer, allocatable :: itpeu(:) !< hulp (numl) edge oriented,
                                                         !! 1,*=boundary typ, see type indicator kbndz(4,*) below
+   integer, allocatable :: itpeu_ib(:) !< hulp (numl) for storing boundary segment number
    integer, allocatable :: itpenz(:) !< temp (numl) edge oriented,
                                                         !! 1,*=boundary number (nopenbndsect), see type indicator kbndz(5,*) below
    integer, allocatable :: itpenu(:) !< temp (numl) edge oriented,

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/compute/transport.f90

@@ -153,10 +153,10 @@ subroutine transport()
       if (jased > 0 .and. jased < 4) then
 
          if (nbndz > 0) then
-            call setequilibriumsedimentbnds(nbndz, 6, kbndz, kbanz, 0)
+            call setequilibriumsedimentbnds(nbndz, 8, kbndz, kbanz, 0)
          end if
          if (nbndu > 0) then
-            call setequilibriumsedimentbnds(nbndu, 6, kbndu, kbanu, 1)
+            call setequilibriumsedimentbnds(nbndu, 8, kbndu, kbanu, 1)
          end if
 
          do k = 1, nbndsd ! set prescribed sediment boundary conditions

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

@@ -43,7 +43,7 @@ subroutine find_nearest_flownodes(n, xx, yy, names, node_nrs_nearest, jakdtree,
       use messagehandling, only: msgbuf, msg_flush, idlen
       use m_partitioninfo, only: jampi, reduce_int_max, reduce_kobs
       use m_flowgeom, only: nd, xz, yz
-      use m_GlobalParameters, only: INDTP_1D, INDTP_ALL
+      use m_GlobalParameters, only: INDTP_1D, INDTP_ALL, INDTP_2D
       use kdtree2Factory, only: treeglob
       use m_inflowcell, only: inflowcell
       use geometry_module, only: dbdistance
@@ -115,6 +115,10 @@ subroutine find_nearest_flownodes(n, xx, yy, names, node_nrs_nearest, jakdtree,
                   node_nrs_nearest(i) = k
                end if
             end if
+            if (jaoutside == 1 .and. iLocTp == INDTP_2D) then
+               call find_nearest_2D_internal_bruteforce(xx(i), yy(i), node_nrs_nearest(i))
+            end if
+
          end do
       end if
 
@@ -325,5 +329,37 @@ subroutine find_nearest_1D_or_boundary_flownode_bruteforce(x, y, node_nr_nearest
       end if
 
    end subroutine find_nearest_1D_or_boundary_flownode_bruteforce
+
+   !> Find the 2-D internal flownode with the shortest distance to the point [x, y]
+   !! Brute-force approach: simply check all flownodes in the entire grid
+   subroutine find_nearest_2D_internal_bruteforce(x, y, node_nr_nearest)
+      use stdlib_kinds, only: dp
+      use m_flowgeom, only: xz, yz,ndxi
+      use geometry_module, only: dbdistance
+      use m_sferic, only: jsferic, jasfer3D
+      use m_missing, only: dmiss
+
+      implicit none
+
+      real(dp), intent(in) :: x !< x-Coordinate of input point.
+      real(dp), intent(in) :: y !< y-Coordinate of input point.
+      integer, intent(out) :: node_nr_nearest !< Node number of nearest 1D or boundary flow node.
+
+      real(dp)             :: dist, distmin
+      integer              :: i_node
+
+      node_nr_nearest = 0
+
+      distmin = huge(1.0_dp)
+
+      do i_node = 1, ndxi
+          dist = dbdistance(x, y, xz(i_node), yz(i_node), jsferic, jasfer3D, dmiss)
+          if (dist < distmin) then
+              distmin = dist
+              node_nr_nearest = i_node
+          end if
+      end do
+
+   end subroutine find_nearest_2D_internal_bruteforce
 
 end module m_find_flownode

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/timespace/fm_external_forcings.f90

@@ -579,7 +579,7 @@ subroutine findexternalboundarypoints() ! find external boundary points
       if (allocated(ftpet)) then
          deallocate (ftpet)
       end if
-      allocate (kce(nx), ke(nx), kez(nx), keu(nx), kes(nx), ketm(nx), kesd(nx), keuxy(nx), ket(nx), ken(nx), ke1d2d(nx), keg(nx), ked(nx), kep(nx), kedb(nx), keklep(nx), kevalv(nx), kegs(nx), kegen(nx), itpez(nx), itpenz(nx), itpeu(nx), itpenu(nx), kew(nx), ftpet(nx), stat=ierr)
+      allocate (kce(nx), ke(nx), kez(nx), keu(nx), kes(nx), ketm(nx), kesd(nx), keuxy(nx), ket(nx), ken(nx), ke1d2d(nx), keg(nx), ked(nx), kep(nx), kedb(nx), keklep(nx), kevalv(nx), kegs(nx), kegen(nx), itpez(nx), itpenz(nx), itpeu(nx), itpeu_ib(nx), itpenu(nx), kew(nx), ftpet(nx), stat=ierr)
       call aerr('kce(nx), ke(nx), kez(nx), keu(nx), kes(nx), ketm(nx), kesd(nx), keuxy(nx), ket(nx), ken(nx), ke1d2d(nx), keg(nx), ked(nx), kep(nx), kedb(nx), keklep(nx), kevalv(nx), kegs(nx), kegen(nx), itpez(nx), itpenz(nx), itpeu(nx) , itpenu(nx), kew(nx), ftpet(nx)', ierr, 17 * nx)
       kce = 0
       ke = 0
@@ -602,6 +602,7 @@ subroutine findexternalboundarypoints() ! find external boundary points
       itpez = 0
       itpenz = 0
       itpeu = 0
+      itpeu_ib = 0
       itpenu = 0
       kew = 0
       ftpet = 1e6_dp
@@ -953,12 +954,13 @@ subroutine read_location_files_from_boundary_blocks(filename, nx, kce, num_bc_in
             call prop_get(node_ptr, '', 'bndBlDepth', blDepth)
 
             if (group_ok) then
+               num_bc_ini_blocks = num_bc_ini_blocks + 1
                if (rrtolb > 0.0_dp) then
-                  call processexternalboundarypoints(quantity, location_file, filetype, return_time, nx, kce, numz, numu, nums, numtm, numsd, numt, numuxy, numn, num1d2d, numqh, numw, numtr, numsf, rrtolrel=(1 + 2 * rrtolb) / (1 + 2 * rrtol), tfc=transformcoef, width1D=width1D, blDepth=blDepth)
+                  call processexternalboundarypoints(quantity, location_file, filetype, return_time, nx, kce, numz, numu, nums, numtm, numsd, numt, numuxy, numn, num1d2d, numqh, numw, numtr, numsf, rrtolrel=(1 + 2 * rrtolb) / (1 + 2 * rrtol), tfc=transformcoef, width1D=width1D, blDepth=blDepth,num_bc_ini_blocks = num_bc_ini_blocks)
                else
-                  call processexternalboundarypoints(quantity, location_file, filetype, return_time, nx, kce, numz, numu, nums, numtm, numsd, numt, numuxy, numn, num1d2d, numqh, numw, numtr, numsf, rrtolrel=1.0_dp, tfc=transformcoef, width1D=width1D, blDepth=blDepth)
+                  call processexternalboundarypoints(quantity, location_file, filetype, return_time, nx, kce, numz, numu, nums, numtm, numsd, numt, numuxy, numn, num1d2d, numqh, numw, numtr, numsf, rrtolrel=1.0_dp, tfc=transformcoef, width1D=width1D, blDepth=blDepth,num_bc_ini_blocks = num_bc_ini_blocks)
                end if
-               num_bc_ini_blocks = num_bc_ini_blocks + 1
+
             end if
 
             file_ok = file_ok .and. group_ok
@@ -1008,7 +1010,7 @@ end subroutine appendrettime
    subroutine processexternalboundarypoints(qid, filename, filetype, return_time, nx, kce, &
                                             numz, numu, nums, numtm, numsd, numt, numuxy, numn, num1d2d, &
                                             numqh, numw, numtr, numsf, rrtolrel, tfc, &
-                                            width1D, blDepth) ! helper for finding external boundary points
+                                            width1D, blDepth, num_bc_ini_blocks) ! helper for finding external boundary points
       use m_netw
       use m_flow, qid_flow => qid, filetype_flow => filetype
       use m_flowgeom
@@ -1040,13 +1042,14 @@ subroutine processexternalboundarypoints(qid, filename, filetype, return_time, n
       real(kind=dp), intent(in) :: return_time
       integer, intent(in) :: numz, numu, nums, numtm, numsd, & !
                              numt, numuxy, numn, num1d2d, numw, numtr, numsf !
+
       integer, intent(inout) :: numqh
       real(kind=dp), intent(in) :: rrtolrel !< To enable a more strict rrtolerance value than the global rrtol. Measured w.r.t. global rrtol.
 
       real(kind=dp), dimension(NUMGENERALKEYWRD), optional, intent(in) :: tfc
       real(kind=dp), optional, intent(in) :: width1D !< Optional custom width for boundary flow link.
       real(kind=dp), optional, intent(in) :: blDepth !< Optional custom bed level depths below water level boundaries's initial value for boundary points.
-
+      integer      , optional, intent(in) :: num_bc_ini_blocks 
       character(len=256) :: qidfm !
       integer :: itpbn
       character(len=NAMTRACLEN) :: tracnam, sfnam, qidnam
@@ -1169,7 +1172,8 @@ subroutine processexternalboundarypoints(qid, filename, filetype, return_time, n
             itpbn = 9
          end if
 
-         itpeu(nbndu + 1:nbndu + numu) = itpbn
+         itpeu   (nbndu + 1:nbndu + numu) = itpbn
+         if (present(num_bc_ini_blocks) ) itpeu_ib(nbndu + 1:nbndu + numu) = num_bc_ini_blocks
 
          call addopenbndsection(numu, keu(nbndu + 1:nbndu + numu), filename, IBNDTP_U)
 
@@ -1890,7 +1894,7 @@ subroutine setup(iresult)
 
       integer :: ierr
       integer :: k, L, LF, KB, KBI, N, K2, iad, numnos, isf, mx, itrac
-      integer, parameter :: N4 = 6
+      integer, parameter :: N4 = 8
       character(len=256) :: rec
       integer :: tmp_nbndu, tmp_nbndt, tmp_nbndn
 
@@ -2101,8 +2105,9 @@ subroutine setup(iresult)
             kbndu(3, k) = Lf
             kbndu(4, k) = itpeu(k)
             kbndu(5, k) = itpenu(k)
-            kbndu(6, k) = ftpet(k) ! riemann relaxation time
-
+            kbndu(6, k) = ftpet(k)        ! riemann relaxation time
+            kbndu(7, k) = itpeu_ib(k)     ! boundary segment number
+            kbndu(8, k) = 1               ! factor for multiplification of discharge boundaries 
             lnxbnd(Lf - lnxi) = itpenu(k)
 
             do n = 1, nd(kbi)%lnx
@@ -3147,4 +3152,69 @@ subroutine allocatewindarrays()
 
    end subroutine allocatewindarrays
 
+   subroutine det_net2pli_position(pliFile,position)
+     ! Fuction: determines position of network relative to a boundary pli
+     !          position =  1: network on the right/on top of a pli
+     !          position = -1: network on the left /below  of a pli
+     use m_filez,            only: oldfil, doclose
+     use m_polygon
+     use m_reapol,           only: reapol
+     use m_missing
+     use m_GlobalParameters, only: INDTP_2D
+     use m_find_flownode,    only: find_nearest_flownodes
+     use m_flowgeom,         only: xz, yz
+     implicit none
+
+      character(len=*), intent(in) :: pliFile
+      integer         , intent(out):: position
+
+      integer                        :: mpli, i_first, i_last, jakdtree ! i_bndpnt
+      real(kind=dp)                  :: xpli_centre, ypli_centre, xbnd_centre, ybnd_centre, vx,vy, wx, wy,cross ! dist, distmin,    
+
+      character(len=5), dimension(1) :: tmpname
+      integer         , dimension(1) :: kbnd_centre
+      real(kind=dp)   , dimension(1) :: x_tmp, y_tmp   
+
+      ! Read the boundary pli
+      call oldfil(mpli, pliFile)
+      call reapol(mpli, 0)
+      call doclose(mpli)
+
+      ! determine centre point of boundary pli 
+      i_first = floor(npl/2.0_dp)
+      i_last  = i_first
+      if (mod(npl,2) == 0)  i_last = i_first + 1
+      xpli_centre = (xpl(i_first) + xpl(i_last))/2.0_dp
+      ypli_centre = (ypl(i_first) + ypl(i_last))/2.0_dp
+
+      ! Determine closest boundary point (probably more ellegant to use find_nearest_flownode, for now keep it simple)
+      tmpname(1)      = 'dummy'
+      x_tmp(1)        = xpli_centre
+      y_tmp(1)        = ypli_centre
+      kbnd_centre(1)  = 0
+      jakdtree        = 0
+      call find_nearest_flownodes(1, x_tmp, y_tmp, tmpname(1),kbnd_centre(1), jakdtree, 1, INDTP_2D)
+      ! Only look for internal 2D points, Boundary pli might be within network
+      xbnd_centre = xz(kbnd_centre(1))
+      ybnd_centre = yz(kbnd_centre(1))
+
+      ! Determine cross product of pli and centre point bnd
+      ! orientation vector pli
+      vx = xpl(i_first + 1) - xpl(i_first)
+      vy = ypl(i_first + 1) - ypl(i_first)
+      ! orientation vector from centre bnd point to firts pli point
+      wx = xbnd_centre - xpl(i_first)
+      wy = ybnd_centre - ypl(i_first)
+      ! cross product
+      cross = (vx*wy) - (vy*wx)
+
+      ! cross < 0 ==> Model right of pli or obove pli
+      if (cross < 0.0_dp) then
+         position = 1
+      else if (cross > 0.0_dp) then 
+         position = -1
+      end if
+
+   end subroutine det_net2pli_position
+
 end module fm_external_forcings

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/timespace/fm_external_forcings_init.f90

@@ -365,6 +365,7 @@ function init_boundary_forcings(block_ptr, base_dir, file_name, group_name, itpe
       use string_module, only: strcmpi
       use properties, only: prop_get
       use unstruc_files, only: resolvePath
+      use fm_external_forcings_data, only: kbndu      
 
       type(tree_data), pointer, intent(in) :: block_ptr !< Pointer to boundary block in extforce file; child node of the extforce file tree
       character(len=*), intent(in) :: base_dir !< Base directory of the ext file
@@ -380,7 +381,7 @@ function init_boundary_forcings(block_ptr, base_dir, file_name, group_name, itpe
       character(len=INI_VALUE_LEN) :: location_file, quantity, forcing_file, property_name, property_value
       type(tree_data), pointer :: key_value_ptr
       character(len=300) :: error_message
-      integer :: operand
+      integer :: operand, ibndu, facdis
       logical :: is_successful
       integer :: method, num_items_in_block, j
 
@@ -502,6 +503,19 @@ function init_boundary_forcings(block_ptr, base_dir, file_name, group_name, itpe
                else
                   is_successful = addtimespacerelation_boundaries(quantity, location_file, filetype=filetype, method=method, &
                                                                   operand=operand, forcing_file=forcing_file)
+                  ! TK_temp: Determine sign needed for discharge boundaries from history file
+                  if (quantity == 'dischargebnd') then
+                     facdis = 1
+                     if (index(trim(forcing_file)//'|', '_his.nc|') > 0 ) then
+                        call det_net2pli_position(location_file,facdis)
+                     end if
+                     ! fill kbndu(8,:) with multiplification discharge boundaries
+                     do ibndu = 1, size(kbndu,2)
+                        if (kbndu(7,ibndu) == ib) then
+                            kbndu(8,ibndu) = facdis
+                        end if
+                     end do
+                 end if
                end if
                res = res .and. is_successful ! Remember any previous errors.
                operand = OPERAND_UNKNOWN

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/timespace/fm_external_forcings_update_boundaries.f90

@@ -58,13 +58,15 @@ module subroutine set_external_forcings_boundaries(time, iresult)
       use m_oned_functions
       use m_obs_on_flowgeom, only: obs_on_flowgeom
       use unstruc_messages, only: callback_msg
+!      use m_filez
 
       implicit none
 
       real(kind=dp), intent(in) :: time !< Current simulation time (s)
       integer, intent(out) :: iresult !< Integer error status
 
       integer :: i, n, k2, kb, L, itrac, isf
+!      integer, save :: minp                         
       real(kind=dp) :: dQ
 
       iresult = DFM_EXTFORCERROR
@@ -165,6 +167,7 @@ module subroutine set_external_forcings_boundaries(time, iresult)
 
       if (item_dischargebnd /= ec_undef_int) then
          success = ec_gettimespacevalue(ecInstancePtr, item_dischargebnd, irefdate, tzone, tunit, time)
+         zbndq = zbndq*kbndu(8,:)
          if (.not. success) then
             goto 888
          end if
@@ -175,6 +178,15 @@ module subroutine set_external_forcings_boundaries(time, iresult)
          if (.not. success) then
             goto 888
          end if
+
+!         if (time0 == 0.0) then
+!            call newfil(minp,' compare_bigSmall_bndForcing/fromEcModule.tim')
+!         end if
+
+!         write(minp,'(7f12.3)') time/60.0,zbnds(4),zbnds(74),zbnds(5),zbnds(75),zbnds(1),zbnds(71)
+!         if (time == tstop_user) then
+!            call doclose(minp)
+!         end if
       end if
 
       if (nbndTM > 0) then