speed up on streamorder domain decomposition

nmizukami · Jun 25, 2019 · e27caca · e27caca
1 parent 66c62b8
commit e27caca
Showing 1 changed file with 40 additions and 21 deletions.
diff --git a/route/build/src/domain_decomposition.f90 b/route/build/src/domain_decomposition.f90
@@ -273,7 +273,6 @@ subroutine omp_domain_decomposition_stro(nSeg, structNTOPO, river_basin_out, ier
 
          associate(upIndex => structNTOPO(ixOutlets(iOut))%var(ixNTOPO%allUpSegIndices)%dat)
 
-         ixUpSeg_tmp = 0
          nSegStreamOrder = 0
          do iSeg = 1, nUpSegs
            if (streamOrder(upIndex(iSeg))/=ix) cycle
@@ -324,30 +323,48 @@ subroutine omp_domain_decomposition_stro(nSeg, structNTOPO, river_basin_out, ier
 
          associate(upIndex => structNTOPO(ixTribOutlet(iTrib))%var(ixNTOPO%allUpSegIndices)%dat)
 
-         ixUpSeg_tmp = 0
-         nSegStreamOrder = 0
-         do iSeg = 1, nUpSegs
-           if (streamOrder(upIndex(iSeg))/=ix) cycle
-           nSegStreamOrder = nSegStreamOrder + 1
-           ixUpSeg_tmp(nSegStreamOrder) = upIndex(iSeg)
-         enddo
+         if (ix==1) then
 
-         if (allocated(subSegOrder)) deallocate(subSegOrder)
-         allocate(subSegOrder(nSegStreamOrder), stat=ierr)
-         if(ierr/=0)then; message=trim(message)//'problem allocating [subSegOrder]'; return; endif
+           if (allocated(subSegOrder)) deallocate(subSegOrder)
+           allocate(subSegOrder(nUpSegs), stat=ierr)
+           if(ierr/=0)then; message=trim(message)//'problem allocating [subSegOrder]'; return; endif
 
-         if (allocated(ixUpSeg)) deallocate(ixUpSeg)
-         allocate(ixUpSeg(nSegStreamOrder), stat=ierr)
-         if(ierr/=0)then; message=trim(message)//'problem allocating [ixUpSeg]'; return; endif
+           allocate(branch(iTrib+nOut)%segIndex(nUpSegs), stat=ierr)
+           if(ierr/=0)then; message=trim(message)//'problem allocating [branch(iTrib)%segIndex]'; return; endif
 
-         allocate(branch(iTrib+nOut)%segIndex(nSegStreamOrder), stat=ierr)
-         if(ierr/=0)then; message=trim(message)//'problem allocating [branch(iTrib)%segIndex]'; return; endif
+           call indexx(rankSegOrder(upIndex), subSegOrder)
 
-         ixUpSeg = ixUpSeg_tmp(1:nSegStreamOrder)
-         call indexx(rankSegOrder(ixUpSeg), subSegOrder)
+           branch(iTrib+nOut)%segIndex = upIndex(subSegOrder)
+           branch(iTrib+nOut)%nRch     = nUpSegs
 
-         branch(iTrib+nOut)%nRch     = nSegStreamOrder
-         branch(iTrib+nOut)%segIndex = ixUpSeg(subSegOrder)
+         else
+
+           nSegStreamOrder = 0
+           do iSeg = 1, nUpSegs
+             if (streamOrder(upIndex(iSeg))/=ix) cycle
+             nSegStreamOrder = nSegStreamOrder + 1
+             ixUpSeg_tmp(nSegStreamOrder) = upIndex(iSeg)
+           enddo
+
+           if (allocated(subSegOrder)) deallocate(subSegOrder)
+           allocate(subSegOrder(nSegStreamOrder), stat=ierr)
+           if(ierr/=0)then; message=trim(message)//'problem allocating [subSegOrder]'; return; endif
+
+           if (allocated(ixUpSeg)) deallocate(ixUpSeg)
+           allocate(ixUpSeg(nSegStreamOrder), stat=ierr)
+           if(ierr/=0)then; message=trim(message)//'problem allocating [ixUpSeg]'; return; endif
+
+           allocate(branch(iTrib+nOut)%segIndex(nSegStreamOrder), stat=ierr)
+           if(ierr/=0)then; message=trim(message)//'problem allocating [branch(iTrib)%segIndex]'; return; endif
+
+           ixUpSeg = ixUpSeg_tmp(1:nSegStreamOrder)
+
+           call indexx(rankSegOrder(ixUpSeg), subSegOrder)
+
+           branch(iTrib+nOut)%nRch     = nSegStreamOrder
+           branch(iTrib+nOut)%segIndex = ixUpSeg(subSegOrder)
+
+         endif
 
          nSegBranch(iTrib+nOut) = branch(iTrib+nOut)%nRch
 
@@ -448,6 +465,7 @@ subroutine classify_river_basin(nDivs,         & ! input:  number of divisions (
 
    ! populate domain(:)%hruIndex
    if (present(nContribHRU)) nContribHRU = 0 ! total number of HRUs that contribute to the reach
+
    do ix=1,nDomains
      associate (ixSeg => domains_out(ix)%segIndex)
 
@@ -470,8 +488,10 @@ subroutine classify_river_basin(nDivs,         & ! input:  number of divisions (
        ix2 = ix1+nHruLocal(iSeg)-1
        domains_out(ix)%hruIndex(ix1:ix2) = structNTOPO(ixSeg(iSeg))%var(ixNTOPO%hruContribIx)%dat(:)
      enddo
+
      deallocate(nHruLocal, stat=ierr)
      if(ierr/=0)then; message=trim(message)//'problem deallocating [nHruLocal]'; return; endif
+
      end associate
    enddo
 
@@ -820,5 +840,4 @@ subroutine basin_order(nSeg, structNTOPO_in, domains_omp, nDomain_omp, river_bas
 
  end subroutine basin_order
 
-
 end module domain_decomposition