Fm/task/unst 9257 snap longculvert 2

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

@@ -37,11 +37,14 @@ module m_longculverts_data
 
    private
 
+   public :: is_2D2D_longculvertlink
+
    !> Type definition for longculvert data.
    type, public :: t_longculvert
       character(len=IdLen) :: id
       character(len=IdLen) :: branchid !< if newculverts, corresponding network branch
       character(len=IdLen) :: csdefid !< if newculverts, corresponding network crossdef
+      character(len=IdLen) :: contactId !< if newculverts, corresponding 2D-2D contact, cannot exist toghether with branchid
       integer :: numlinks !< Number of links of the long culvert
       integer, dimension(:), allocatable :: netlinks !< Net link numbers of the long culvert
       integer, dimension(:), allocatable :: flowlinks !< Flow link numbers of the long culvert
@@ -67,4 +70,24 @@ module m_longculverts_data
    integer, public :: nlongculverts !< Number of longculverts
    logical, public :: newculverts
 
+contains
+
+   !> simple routine which checks if a given flowlink L is part of a 2D-2D longculvert. Lives here to avoid cyclic dependency.
+   elemental subroutine is_2D2D_longculvertlink(L, res, i)
+      integer, intent(in) :: L !< Flowlink number
+      logical, intent(out) :: res !< Result: true if L is part of a 2D-2D longculvert
+      integer, intent(out) :: i !< Index of the longculvert in longculverts derived type array
+      res = .false.
+      do i = 1, nlongculverts
+         if (allocated(longculverts(i)%netlinks)) then
+            if (size(longculverts(i)%netlinks) == 1) then
+               if ((longculverts(i)%netlinks(1) == L)) then
+                  res = .true.
+                  exit
+               end if
+            end if
+         end if
+      end do
+   end subroutine is_2D2D_longculvertlink
+
 end module m_longculverts_data

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

@@ -41,6 +41,7 @@ module m_save_ugrid_state
    character(len=ug_idsLongNamesLen), allocatable :: nbranchlongnames(:), nnodelongnames(:), nodelongnames(:)
    character(len=255) :: network1dname, mesh2dname, mesh1dname, contactname_1D2D, contactname_2D2D !MAXSTRLEN = 255
    character(len=ug_idsLen), allocatable :: mesh1dNodeIds(:)
+   character(len=ug_idsLen), dimension(:), allocatable :: contactids_2D2D ! Long names of 2D-2D contacts
    integer, allocatable, dimension(:) :: mesh1dUnmergedToMerged(:)
    !integer, allocatable, dimension(:)                 :: mesh1dMergedToUnMerged(:)
    integer :: numMesh1dBeforeMerging
@@ -49,7 +50,7 @@ module m_save_ugrid_state
    integer :: contactnlinks !< Total number of links in all mesh contacts (typically we'll have one mesh contact with many netlinks part of it)
    integer, allocatable :: contact_cell_idx(:, :) !< cell numbers this contact connects (2D for 2D-2D links, 1D for 1D-2D links). Dimension(2, contactnlinks). Second row is always 2D
    type(t_hashlist) :: hashlist_contactids !< Hash list for quick search for contact ids.
-
+   
 !> Sets ALL (scalar) variables in this module to their default values.
 !! For a reinit prior to flow computation, only call reset_save_ugrid_state() instead.
 contains

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_data/partition.F90

@@ -605,10 +605,11 @@ subroutine partition_make_domain(idmn, numlay_cell, numlay_node, jacells, ierror
       use messagehandling, only: mess, level_warn
       use dfm_error, only: dfm_genericerror, dfm_noerr
       use m_polygon, only: npl
-      use network_data, only: lc, numl, kn, link_2d, lnn, lne, nump1d2d, netstat, netstat_ok, numk, cellmask, lperm, netcell, numl1d
+      use network_data, only: lc, numl, kn, link_2d, lnn, lne, nump1d2d, nump, netstat, netstat_ok, numk, cellmask, lperm, netcell, numl1d
       use m_alloc, only: realloc
       use gridoperations, only: findcells
       use m_remove_masked_netcells, only: remove_masked_netcells
+      use m_save_ugrid_state, only: contactnetlinks, contactids_2D2D, hashlist_contactids
       implicit none
 
       integer, intent(in) :: idmn !< domain number
@@ -619,7 +620,7 @@ subroutine partition_make_domain(idmn, numlay_cell, numlay_node, jacells, ierror
 
       integer :: ic1, ic2, L
       logical :: domain_needs_cell_1, domain_needs_cell_2
-      integer :: i
+      integer :: i, icontact, i_valid_contact
       integer, dimension(:, :), allocatable :: lne_org
       integer :: i_old
       character(len=128) :: message
@@ -672,6 +673,23 @@ subroutine partition_make_domain(idmn, numlay_cell, numlay_node, jacells, ierror
          end do
       end if
 
+      if (allocated(contactnetlinks) .and. allocated(contactids_2D2D)) then
+         i_valid_contact = 0
+         call realloc(contactids_2D2D, size(contactnetlinks), keepExisting=.false.)
+         do icontact = 1, size(contactnetlinks)
+            L = contactnetlinks(icontact)
+            if (L > 0 .and. L <= numL) then
+               if (Lc(L) == 1) then !> valid link
+                  if (lne_org(1, L) < nump .and. lne_org(2, L) < nump) then !> both cells are real 2D cells
+                     i_valid_contact = i_valid_contact + 1
+                     contactids_2D2D(i_valid_contact) = hashlist_contactids%id_list(icontact)
+                  end if
+               end if
+            end if
+         end do
+         call realloc(contactids_2D2D, i_valid_contact, keepExisting=.true.)
+      end if
+
 !     physically remove nodes and links from network
       NPL = 0 ! number of polygon points, we set no polygon
       call findcells(100000) ! output link permutation array "Lperm" (only used if jacells.eq.1)

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

@@ -44,7 +44,7 @@ module unstruc_model
    use properties, only: prop_get, prop_file, tree_create, tree_destroy, max_prop_length
    use m_waveconst
 
-   implicit none
+   implicit none(external)
 
    !> The version number of the MDU File format: d.dd, [config_major].[config_minor], e.g., 1.03
     !!
@@ -428,6 +428,7 @@ subroutine loadModel(filename)
       use m_realan, only: realan
       use m_filez, only: oldfil
       use unstruc_messages, only: threshold_abort
+      use m_readCrossSections, only: readCrossSectionDefinitions
 
       character(*), intent(inout) :: filename !< Name of file to be read (in current directory or with full path).
 
@@ -479,11 +480,19 @@ subroutine loadModel(filename)
       network%sferic = jsferic == 1
 
       threshold_abort = LEVEL_FATAL
-      if (istat == 0 .and. jadoorladen == 0 .and. network%numk > 0 .and. network%numl > 0) then
-         timerHandle = 0
-         call timstrt('Read 1d attributes', timerHandle)
-         call read_1d_attributes(md_1dfiles, network)
-         call timstop(timerHandle)
+      if (istat == 0 .and. jadoorladen == 0) then
+         if (network%numk > 0 .and. network%numl > 0) then
+            timerHandle = 0
+            call timstrt('Read 1d attributes', timerHandle)
+            call read_1d_attributes(md_1dfiles, network)
+            call timstop(timerHandle)
+            ! Always read cross section definitions if specified (needed for long culverts)
+         else if (len_trim(md_1dfiles%cross_section_definitions) > 0) then
+            call readCrossSectionDefinitions(network, md_1dfiles%cross_section_definitions)
+            if (network%CSDefinitions%Count < 1) then
+               call SetMessage(LEVEL_WARN, 'No Cross_Section Definitions Found in file '//trim(md_1dfiles%cross_section_definitions))
+            end if
+         end if
       end if
 
       timerHandle = 0

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_io/unstruc_netcdf.f90

@@ -11865,6 +11865,8 @@ subroutine unc_write_net_ugrid2(ncid, id_tsp, janetcell, jaidomain, jaiglobal_s)
       use m_set_nod_adm
       use m_inquire_link_type, only: is_valid_2d2d_netlink, is_valid_1d2d_netlink, is_valid_1D_netlink, count_1D_edges, count_1D_nodes
       use m_cell_geometry, only: blcell
+      use m_longculverts_data, only: longculverts, is_2D2D_longculvertlink
+
       implicit none
 
       integer, intent(in) :: ncid !< NetCDF file id
@@ -11884,7 +11886,8 @@ subroutine unc_write_net_ugrid2(ncid, id_tsp, janetcell, jaidomain, jaiglobal_s)
       integer :: ierr
       integer :: i, k, k1, k2, numl2d, numk2d, L, Lnew, nv, n1, n2, n
       integer :: num_1d_nodes, node_index
-      logical :: jaInDefine
+      logical :: jaInDefine, is_longculvert_contact
+      integer :: longculvertindex
       integer :: id_zf
       real(kind=hp), allocatable :: xn(:), yn(:), zn(:), xe(:), ye(:), zf(:)
       integer :: n1dedges, n1d2dcontacts, n2d2dcontacts, start_index
@@ -11947,12 +11950,16 @@ subroutine unc_write_net_ugrid2(ncid, id_tsp, janetcell, jaidomain, jaiglobal_s)
       if (n2d2dcontacts > 0) then
          allocate (contacts_2D2D(2, n2d2dcontacts))
          call realloc(contacttype_2D2D, n2d2dcontacts, keepExisting=.false., fill=5)
-
+         call realloc(contactids_2D2D, n2d2dcontacts, keepExisting=.true., fill='')
          do i = 1, n2d2dcontacts
             L = temp_indices(i)
             n1 = abs(lne(1, L))
             n2 = abs(lne(2, L))
             contacts_2D2D(1:2, i) = [n1, n2]
+            call is_2D2D_longculvertlink(L, is_longculvert_contact, longculvertindex)
+            if (is_longculvert_contact) then
+               contactids_2D2D(i) = longculverts(longculvertindex)%contactID !< reuse branchid for long culvert contacts
+            end if
          end do
       end if
 
@@ -12283,8 +12290,9 @@ subroutine unc_write_net_ugrid2(ncid, id_tsp, janetcell, jaidomain, jaiglobal_s)
          if (n2d2dcontacts > 0) then
             ierr = ug_def_mesh_contact(ncid, id_tsp%meshcontact_2D2D, trim(contactname_2D2D), n2d2dcontacts, id_tsp%meshids2d, id_tsp%meshids2d, UG_LOC_FACE, UG_LOC_FACE, start_index)
             ierr = nf90_enddef(ncid)
+
             ! Put the contacts
-            ierr = ug_put_mesh_contact(ncid, id_tsp%meshcontact_2D2D, contacts_2D2D(1, :), contacts_2D2D(2, :), contacttype_2D2D)
+            ierr = ug_put_mesh_contact(ncid, id_tsp%meshcontact_2D2D, contacts_2D2D(1, :), contacts_2D2D(2, :), contacttype_2D2D, contactsids=contactids_2D2D)
             ierr = nf90_redef(ncid) ! TODO: AvD: I know that all this redef is slow. Split definition and writing soon.
          end if