diff --git a/src/tools_gpl/waq_tools/waqmerge/waqmerge.f90 b/src/tools_gpl/waq_tools/waqmerge/waqmerge.f90 index 87903a9442..4512ef2657 100644 --- a/src/tools_gpl/waq_tools/waqmerge/waqmerge.f90 +++ b/src/tools_gpl/waq_tools/waqmerge/waqmerge.f90 @@ -35,6 +35,7 @@ program waqmerge use waqmerge_version_module, only: getfullversionstring_waqmerge use m_date_time_utils_external, only: fill_in_date_time use m_file_path_utils, only: extract_file_extension + use m_merge_domains implicit none diff --git a/src/utils_lgpl/io_hyd/packages/io_hyd/src/merge_domains.f90 b/src/utils_lgpl/io_hyd/packages/io_hyd/src/merge_domains.f90 index 2987102686..88b7d93cb2 100644 --- a/src/utils_lgpl/io_hyd/packages/io_hyd/src/merge_domains.f90 +++ b/src/utils_lgpl/io_hyd/packages/io_hyd/src/merge_domains.f90 @@ -26,6 +26,11 @@ !------------------------------------------------------------------------------- ! ! + module m_merge_domains + implicit none + + contains + subroutine merge_domains(hyd, domain_hyd_coll) ! function : merge the domains, make pointers to final domain @@ -40,15 +45,15 @@ subroutine merge_domains(hyd, domain_hyd_coll) ! declaration of the arguments type(t_hydrodynamics) :: hyd ! description of the hydrodynamics - type(t_hydrodynamics_collection) :: domain_hyd_coll ! description of all domain hydrodynamics + type(t_hydrodynamics_collection) :: domain_hyd_coll ! description of all domain hydrodynamics ! local declarations integer :: n_domain ! number of domains integer :: i_domain ! index in collection integer :: idmn ! flow like domain index (0:n_domain-1) - type(t_hydrodynamics), pointer :: d_hyd ! description of one domain hydrodynamics - type(t_hydrodynamics), pointer :: l_hyd ! description of a linked domain hydrodynamics + type(t_hydrodynamics), pointer :: d_hyd ! description of one domain hydrodynamics + type(t_hydrodynamics), pointer :: l_hyd ! description of a linked domain hydrodynamics integer :: nosegl ! total number of segments per layer integer :: num_boundary_conditions ! total number of boundaries integer :: nobndl ! total number of boundaries per layer @@ -68,9 +73,9 @@ subroutine merge_domains(hyd, domain_hyd_coll) integer :: isect ! index of section integer :: no_bnd ! number of boundaries in section integer :: i_bnd ! index of boundary - type(t_openbnd_section), pointer :: openbndsect ! single section - type(t_open_boundary_line),pointer :: openbndlin ! single open boundary lin - type(t_openbnd_section) :: new_sect ! single section new + type(t_openbnd_section), pointer :: openbndsect ! single section + type(t_open_boundary_line),pointer :: openbndlin ! single open boundary lin + type(t_openbnd_section) :: new_sect ! single section new logical :: bnd_active ! if a boundary is active integer :: iret ! return value integer :: ik ! node counter @@ -113,7 +118,7 @@ subroutine merge_domains(hyd, domain_hyd_coll) hyd%crs = d_hyd%crs hyd%conv_type = d_hyd%conv_type hyd%conv_version = d_hyd%conv_version - + ! copy waqgeom layer information from first domain hyd%waqgeom%num_layers = d_hyd%waqgeom%num_layers hyd%waqgeom%numtopsig = d_hyd%waqgeom%numtopsig @@ -453,6 +458,9 @@ subroutine merge_domains(hyd, domain_hyd_coll) end do end do + ! determine the stable numbering for the boundary cells + call sort_and_number_boundaries( domain_hyd_coll%hyd_pnts(1:domain_hyd_coll%current_size), num_boundary_conditions ) + ! exchanges ! gather the exchanges per layer, and in such a way that the internal exchanges come first, and overlap with the edge face table num_exchanges_u_dir = 0 @@ -489,6 +497,7 @@ subroutine merge_domains(hyd, domain_hyd_coll) do iq = 1, d_hyd%num_exchanges_u_dir ip1 = d_hyd%ipoint(1,iq) ip2 = d_hyd%ipoint(2,iq) + if ( ip1 .lt. 0 .and. ip2 .ge. min_seg .and. ip2 .le. max_seg) then ! ip1 is a boundary and ip2 is in this layer if(d_hyd%idomain(ip2) .eq. idmn .and. (.not.d_hyd%ispoint_bnd(abs(ip1)))) then @@ -497,8 +506,9 @@ subroutine merge_domains(hyd, domain_hyd_coll) d_hyd%iglobal_link(iq) = num_exchanges_u_dir if (abs(ip1) .le. d_hyd%nobndl) then num_boundary_conditions = num_boundary_conditions + 1 - d_hyd%iglobal_bnd(-ip1) = -num_boundary_conditions - call renum_bnd(d_hyd%openbndsect_coll,ip1,-num_boundary_conditions) + !d_hyd%iglobal_bnd(-ip1) = -num_boundary_conditions + !call renum_bnd(d_hyd%openbndsect_coll,ip1,-num_boundary_conditions) + call renum_bnd(d_hyd%openbndsect_coll,ip1,d_hyd%iglobal_bnd(-ip1)) end if end if else if (ip1 .ge. min_seg .and. ip1 .le. max_seg .and. ip2 .lt. 0) then !ip1 between min +1 and max @@ -509,11 +519,13 @@ subroutine merge_domains(hyd, domain_hyd_coll) d_hyd%iglobal_link(iq) = num_exchanges_u_dir if (abs(ip2) .le. d_hyd%nobndl) then num_boundary_conditions = num_boundary_conditions + 1 - d_hyd%iglobal_bnd(-ip2) = -num_boundary_conditions - call renum_bnd(d_hyd%openbndsect_coll,ip2,-num_boundary_conditions) + !d_hyd%iglobal_bnd(-ip2) = -num_boundary_conditions + !call renum_bnd(d_hyd%openbndsect_coll,ip2,-num_boundary_conditions) + call renum_bnd(d_hyd%openbndsect_coll,ip2,d_hyd%iglobal_bnd(-ip2)) end if end if end if + end do end do end do @@ -533,8 +545,9 @@ subroutine merge_domains(hyd, domain_hyd_coll) d_hyd%iglobal_link(iq) = num_exchanges_u_dir if (abs(ip1) .le. d_hyd%nobndl) then num_boundary_conditions = num_boundary_conditions + 1 - d_hyd%iglobal_bnd(-ip1) = -num_boundary_conditions - call renum_bnd(d_hyd%openbndsect_coll,ip1,-num_boundary_conditions) + !d_hyd%iglobal_bnd(-ip1) = -num_boundary_conditions + !call renum_bnd(d_hyd%openbndsect_coll,ip1,-num_boundary_conditions) + call renum_bnd(d_hyd%openbndsect_coll,ip1,d_hyd%iglobal_bnd(-ip1)) end if end if else if (ip1 .gt. 0 .and. ip2 .lt. 0) then @@ -545,8 +558,9 @@ subroutine merge_domains(hyd, domain_hyd_coll) d_hyd%iglobal_link(iq) = num_exchanges_u_dir if (abs(ip2) .le. d_hyd%nobndl) then num_boundary_conditions = num_boundary_conditions + 1 - d_hyd%iglobal_bnd(-ip2) = -num_boundary_conditions - call renum_bnd(d_hyd%openbndsect_coll,ip2,-num_boundary_conditions) + !d_hyd%iglobal_bnd(-ip2) = -num_boundary_conditions + !call renum_bnd(d_hyd%openbndsect_coll,ip2,-num_boundary_conditions) + call renum_bnd(d_hyd%openbndsect_coll,ip1,d_hyd%iglobal_bnd(-ip2)) end if end if end if @@ -603,7 +617,7 @@ subroutine merge_domains(hyd, domain_hyd_coll) endif enddo enddo - + ! For z-layer models, make sure the lower most boundary exist, so that the minimum value of the pointer ! equals the number of boundary conditions per layer times the number of layers if (minval(hyd%ipoint) /= -hyd%num_boundary_conditions) then @@ -1004,7 +1018,7 @@ subroutine merge_domains_old(hyd, domain_hyd_coll) end if end do end do - + ! For z-layer models, make sure the lower most boundary exist, so that the minimum value of the pointer ! equals the number of boundary conditions per layer times the number of layers if (minval(hyd%ipoint) /= -hyd%num_boundary_conditions) then @@ -1198,3 +1212,167 @@ subroutine merge_domains_old(hyd, domain_hyd_coll) end do return end + + subroutine sort_and_number_boundaries( hyd_pnts, num_boundaries ) + + use m_hydmod + + implicit none + + ! Use the properties of the boundary sections to sort + ! the individual boundary cells in such a way that a + ! numbering results that is independent of the partitioning. + ! Note: D-Flow FM does not provide such a numbering itself, + ! unlike with the internal cells. + + type(t_hydrodynamics), dimension(:), intent(inout) :: hyd_pnts + integer, intent(out) :: num_boundaries + + type :: t_cell_properties + character(len=40) :: name + real(kind=kind(1.0d0)), dimension(4) :: coordinates + integer :: domain_index + integer :: section_index + integer :: bound_index + integer :: original_index + integer :: sorted_index = -1 + end type t_cell_properties + + type(t_cell_properties), allocatable, dimension(:) :: cell_list + type(t_cell_properties), allocatable, dimension(:) :: new_cell_list + type(t_cell_properties) :: new_cell + type(t_openbnd_section) :: section + type(t_open_boundary_line) :: bound + integer :: b, c, d, e, s + + ! + ! Gather the information per domain into a single list + ! + allocate( cell_list(0) ) + + e = 0 + do d = 1,size(hyd_pnts) + do s = 1,hyd_pnts(d)%openbndsect_coll%current_size + section = hyd_pnts(d)%openbndsect_coll%openbndsect_pnts(s) + + do b = 1,section%openbndlin_coll%current_size + bound = section%openbndlin_coll%openbndlin_pnts(b) + e = e + 1 + + new_cell%name = section%name + new_cell%coordinates = [bound%x1, bound%y1, bound%x2, bound%y2] + new_cell%domain_index = d + new_cell%section_index = s + new_cell%bound_index = e + new_cell%original_index = bound%ibnd + + cell_list = [cell_list, new_cell] + enddo + enddo + enddo + + ! + ! Sort the list: + ! - First by section name + ! - Then by coordinates (one by one) + ! + do d = 1,size(hyd_pnts) + do s = 1,hyd_pnts(d)%openbndsect_coll%current_size + section = hyd_pnts(d)%openbndsect_coll%openbndsect_pnts(s) + + do b = 1,section%openbndlin_coll%current_size + bound = section%openbndlin_coll%openbndlin_pnts(b) + enddo + enddo + enddo + + call sort_lexicographically( cell_list ) + + call identify_duplicates( cell_list, new_cell_list ) + + ! + ! Fill the new numbering in in the domain list + ! + num_boundaries = size(new_cell_list) + + do c = 1,size(new_cell_list) + d = new_cell_list(c)%domain_index + s = new_cell_list(c)%section_index + b = abs( new_cell_list(c)%original_index ) + + !hyd_pnts(d)%openbndsect_coll%openbndsect_pnts(s)%openbndlin_coll%openbndlin_pnts(b)%ibnd_new = & + hyd_pnts(d)%iglobal_bnd(b) = -new_cell_list(c)%sorted_index + enddo + + contains + + + subroutine sort_lexicographically( cell_list ) + + ! Use a simple bubblesort to sort the boundary segments in a persistent way + + type(t_cell_properties), dimension(:), intent(inout) :: cell_list + + type(t_cell_properties) :: bstore + integer :: i, c + logical :: changed, swap + + changed = .true. + + do while (changed) + changed = .false. + + do i = 2,size(cell_list) + swap = .false. + if ( cell_list(i)%name < cell_list(i-1)%name ) then + swap = .true. + elseif ( cell_list(i)%name == cell_list(i-1)%name ) then + do c = 1,4 + if ( cell_list(i)%coordinates(c) < cell_list(i-1)%coordinates(c) ) then + swap = .true. + elseif ( cell_list(i)%coordinates(c) == cell_list(i-1)%coordinates(c) ) then + cycle + else + exit + endif + enddo + endif + if ( swap ) then + changed = .true. + bstore = cell_list(i) + cell_list(i) = cell_list(i-1) + cell_list(i-1) = bstore + endif + enddo + enddo + end subroutine sort_lexicographically + + + subroutine identify_duplicates( cell_list, new_cell_list ) + + ! Go through the sorted list and remove the duplicates + + type(t_cell_properties), dimension(:), intent(in) :: cell_list + type(t_cell_properties), dimension(:), intent(out), allocatable :: new_cell_list + + integer :: i, j, k + + new_cell_list = cell_list + j = 1 + k = 1 + new_cell_list(1)%sorted_index = k + + do i = 2,size(cell_list) + if ( new_cell_list(i)%name /= new_cell_list(j)%name .or. & + any( new_cell_list(i)%coordinates /= new_cell_list(j)%coordinates ) ) then + j = i ! Register this entry + k = k + 1 + endif + new_cell_list(i)%sorted_index = k + enddo + + end subroutine identify_duplicates + + end subroutine sort_and_number_boundaries + + end module m_merge_domains