From 2e1148174262d65c98bcad60ed57180bbfa66154 Mon Sep 17 00:00:00 2001 From: Naoki Mizukami Date: Thu, 20 Jun 2019 14:50:19 -0600 Subject: [PATCH] put back the routine for upstream threshold based decomposition for OMP. Implemented flow accumulation omp parallel routine. Removed timing for openMP within OMP parallel do loop --- route/build/src/accum_runoff.f90 | 75 +++++--- route/build/src/domain_decomposition.f90 | 217 +++++++++++++++++++++-- route/build/src/irf_route.f90 | 37 +--- route/build/src/kwt_route.f90 | 38 +--- route/build/src/main_route.f90 | 1 + route/build/src/mpi_process.f90 | 11 +- 6 files changed, 267 insertions(+), 112 deletions(-) diff --git a/route/build/src/accum_runoff.f90 b/route/build/src/accum_runoff.f90 index 8c4fd7d6..41ab03ed 100644 --- a/route/build/src/accum_runoff.f90 +++ b/route/build/src/accum_runoff.f90 @@ -18,8 +18,8 @@ MODULE accum_runoff_module ! --------------------------------------------------------------------------------------- ! Public subroutine main driver for basin routing ! --------------------------------------------------------------------------------------- - SUBROUTINE accum_runoff(& - iEns, & ! input: index of runoff ensemble to be processed + SUBROUTINE accum_runoff(iEns, & ! input: index of runoff ensemble to be processed + river_basin, & ! input: river basin information (mainstem, tributary outlet etc.) ixDesire, & ! input: ReachID to be checked by on-screen printing NETOPO_in, & ! input: reach topology data structure RCHFLX_out, & ! inout: reach flux data structure @@ -33,28 +33,32 @@ SUBROUTINE accum_runoff(& ! ! ---------------------------------------------------------------------------------------- - USE nr_utility_module, ONLY : arth + USE dataTypes, ONLY: subbasin_omp ! mainstem+tributary data structures implicit none ! input - integer(i4b), intent(in) :: iens ! runoff ensemble index - integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output - type(RCHTOPO),intent(in), allocatable :: NETOPO_in(:) ! River Network topology + integer(i4b), intent(in) :: iens ! runoff ensemble index + type(subbasin_omp), intent(in), allocatable :: river_basin(:) ! river basin information (mainstem, tributary outlet etc.) + integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output + type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology ! inout - TYPE(STRFLX), intent(inout), allocatable :: RCHFLX_out(:,:)! Reach fluxes (ensembles, space [reaches]) for decomposed domains + TYPE(STRFLX), intent(inout), allocatable :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains ! output - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message ! error message + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message ! input (optional) - integer(i4b), intent(in), optional :: ixSubRch(:) ! subset of reach indices to be processed + integer(i4b), intent(in), optional :: ixSubRch(:) ! subset of reach indices to be processed ! local variables - integer(i4b) :: nSeg ! number of segments in the network - integer(i4b) :: iSeg, jSeg ! reach segment indices - logical(lgt), allocatable :: doRoute(:) ! logical to indicate which reaches are processed - character(len=strLen) :: cmessage ! error message from subroutines - integer*8 :: cr ! rate - integer*8 :: startTime,endTime ! date/time for the start and end of the initialization - real(dp) :: elapsedTime ! elapsed time for the process + integer(i4b) :: nSeg ! number of segments in the network + integer(i4b) :: nTrib ! number of tributaries + integer(i4b) :: nDom ! number of domains defined by e.g., stream order, tributary/mainstem + integer(i4b) :: iSeg, jSeg ! reach segment indices + integer(i4b) :: iTrib, ix ! loop indices + logical(lgt), allocatable :: doRoute(:) ! logical to indicate which reaches are processed + character(len=strLen) :: cmessage ! error message from subroutines + integer*8 :: cr ! rate + integer*8 :: startTime,endTime ! date/time for the start and end of the initialization + real(dp) :: elapsedTime ! elapsed time for the process ierr=0; message='accum_runoff/' call system_clock(count_rate=cr) @@ -78,17 +82,36 @@ SUBROUTINE accum_runoff(& doRoute(:)=.true. ! every reach is on endif - call system_clock(startTime) - - ! compute the sum of all upstream runoff at each point in the river network - do iSeg=1,nSeg - - jSeg = NETOPO_in(iSeg)%RHORDER + nDom = size(river_basin) - if (.not. doRoute(jSeg)) cycle + call system_clock(startTime) - call accum_qupstream(iens, jSeg, ixDesire, NETOPO_in, RCHFLX_out, ierr, cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + do ix = 1,nDom + ! 1. Route tributary reaches (parallel) + ! compute the sum of all upstream runoff at each point in the river network + nTrib=size(river_basin(ix)%branch) + +!$OMP PARALLEL DO schedule(dynamic,1) & +!$OMP private(jSeg, iSeg) & ! private for a given thread +!$OMP private(ierr, cmessage) & ! private for a given thread +!$OMP shared(river_basin) & ! data structure shared +!$OMP shared(doRoute) & ! data array shared +!$OMP shared(NETOPO_in) & ! data structure shared +!$OMP shared(RCHFLX_out) & ! data structure shared +!$OMP shared(ix, iEns, ixDesire) & ! indices shared +!$OMP firstprivate(nTrib) + do iTrib = 1,nTrib + do iSeg=1,river_basin(ix)%branch(iTrib)%nRch + jSeg = river_basin(ix)%branch(iTrib)%segIndex(iSeg) + + if (.not. doRoute(jSeg)) cycle + + call accum_qupstream(iens, jSeg, ixDesire, NETOPO_in, RCHFLX_out, ierr, cmessage) + !if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + end do + end do +!$OMP END PARALLEL DO end do ! looping through stream segments diff --git a/route/build/src/domain_decomposition.f90 b/route/build/src/domain_decomposition.f90 index b8690230..9033ba9f 100644 --- a/route/build/src/domain_decomposition.f90 +++ b/route/build/src/domain_decomposition.f90 @@ -18,10 +18,15 @@ MODULE domain_decomposition implicit none +! common parameters within this module +integer(i4b), parameter :: tributary=1 +integer(i4b), parameter :: mainstem=2 + private public :: mpi_domain_decomposition public :: omp_domain_decomposition +public :: omp_domain_decomposition_stro contains @@ -126,15 +131,49 @@ end subroutine print_screen end subroutine mpi_domain_decomposition - ! *************************************************************** - ! public subroutine: OMP domain decomposition + ! public subroutine: OMP domain decomposition - method 1 ! *************************************************************** subroutine omp_domain_decomposition(nSeg, structNTOPO, river_basin_out, ierr, message) + USE globalData, only: nThreads ! number of threads + + implicit none + ! Input variables + integer(i4b), intent(in) :: nSeg ! number of stream segments + type(var_ilength), allocatable, intent(in) :: structNTOPO(:) ! network topology + ! Output variables + type(subbasin_omp),allocatable, intent(out) :: river_basin_out(:) + integer(i4b), intent(out) :: ierr + character(len=strLen), intent(out) :: message ! error message + ! Local variables + type(subbasin_mpi) :: domains_omp(maxDomain)! domain decomposition data structure (maximum domain is set to maxDomain) + integer(i4b) :: nDomain_omp + character(len=strLen) :: cmessage ! error message from subroutine + + ierr=0; message='omp_domain_decomposition/' + + call classify_river_basin(nThreads, & ! input: number of threads + nSeg, & ! input: number of reaches in the entire river network + structNTOPO, & ! input: river network data structure + domains_omp, & ! output: domain data structure + nDomain_omp, & ! output: number of domains + ierr, cmessage) ! output: error handling + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + call basin_order(nSeg, structNTOPO, domains_omp, nDomain_omp, river_basin_out, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + end subroutine omp_domain_decomposition + + ! *************************************************************** + ! public subroutine: OMP domain decomposition - method2 + ! *************************************************************** + subroutine omp_domain_decomposition_stro(nSeg, structNTOPO, river_basin_out, ierr, message) + ! External modules USE pfafstetter_module, only: lgc_tributary_outlet - + USE dataTypes, only: reach ! reach data structure implicit none ! Input variables integer(i4b), intent(in) :: nSeg ! number of stream segments @@ -145,6 +184,9 @@ subroutine omp_domain_decomposition(nSeg, structNTOPO, river_basin_out, ierr, me character(len=strLen), intent(out) :: message ! error message ! Local variables character(len=strLen) :: cmessage ! error message from subroutine + type(reach), allocatable :: branch(:) + integer(i4b), allocatable :: nSegBranch(:) ! number of reaches in a branch of a stream order + integer(i4b), allocatable :: rankBranch(:) ! ranked branches with a stream order integer(i4b) :: segIndex(nSeg) ! reach indices for reaches within a domain integer(i4b) :: downIndex(nSeg) ! downstream reach indices for reaches within a domain integer(i4b) :: streamOrder(nSeg) ! stream order for reaches within a domain @@ -160,12 +202,12 @@ subroutine omp_domain_decomposition(nSeg, structNTOPO, river_basin_out, ierr, me integer(i4b), allocatable :: subSegOrder(:) ! reach order in a subset domain logical(lgt), allocatable :: isTribOutlet(:) integer(i4b) :: iSeg ! reach loop indices - integer(i4b) :: iTrib,iOut,ix ! loop indices + integer(i4b) :: iTrib,iOut,iBrn,ix ! loop indices integer(i4b) :: nTrib,nOut ! number of tributaries, and basin outlets, respectively integer(i4b) :: nSegStreamOrder ! number of reachs of the same stream order integer(i4b) :: nUpSegs ! numpber of upstream segments - ierr=0; message='omp_domain_decomposition/' + ierr=0; message='omp_domain_decomposition_stro/' do iSeg = 1,nSeg segIndex(iSeg) = structNTOPO(iSeg)%var(ixNTOPO%segIndex)%dat(1) @@ -211,6 +253,12 @@ subroutine omp_domain_decomposition(nSeg, structNTOPO, river_basin_out, ierr, me allocate(river_basin_out(ix)%branch(nTrib+nOut), stat=ierr) if(ierr/=0)then; message=trim(message)//'problem allocating [river_basin_out(ix)%branch]'; return; endif + if (allocated(branch)) deallocate(branch) + if (allocated(rankBranch)) deallocate(rankBranch) + if (allocated(nSegBranch)) deallocate(nSegBranch) + allocate(branch(nTrib+nOut),nSegBranch(nTrib+nOut), rankBranch(nTrib+nOut), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating [branch]'; return; endif + outlet:if (nOut>0) then if (allocated(ixOutlets)) deallocate(ixOutlets) allocate(ixOutlets(nOut), stat=ierr) @@ -241,15 +289,18 @@ subroutine omp_domain_decomposition(nSeg, structNTOPO, river_basin_out, ierr, me allocate(ixUpSeg(nSegStreamOrder), stat=ierr) if(ierr/=0)then; message=trim(message)//'problem allocating [ixUpSeg]'; return; endif - allocate(river_basin_out(ix)%branch(iOut)%segIndex(nSegStreamOrder), stat=ierr) - if(ierr/=0)then; message=trim(message)//'problem allocating [river_basin_out(ix)%branch(iTrib)%segIndex]'; return; endif + allocate(branch(iOut)%segIndex(nSegStreamOrder), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating [branch(iOut)%segIndex]'; return; endif ixUpSeg = ixUpSeg_tmp(1:nSegStreamOrder) call indexx(rankSegOrder(ixUpSeg), subSegOrder) - river_basin_out(ix)%branch(iOut)%nRch = nSegStreamOrder - river_basin_out(ix)%branch(iOut)%segIndex = ixUpSeg(subSegOrder) + branch(iOut)%nRch = nSegStreamOrder + branch(iOut)%segIndex = ixUpSeg(subSegOrder) + + nSegBranch(iOut) = branch(iOut)%nRch + end associate end do @@ -289,22 +340,32 @@ subroutine omp_domain_decomposition(nSeg, structNTOPO, river_basin_out, ierr, me allocate(ixUpSeg(nSegStreamOrder), stat=ierr) if(ierr/=0)then; message=trim(message)//'problem allocating [ixUpSeg]'; return; endif - allocate(river_basin_out(ix)%branch(iTrib+nOut)%segIndex(nSegStreamOrder), stat=ierr) - if(ierr/=0)then; message=trim(message)//'problem allocating [river_basin_out(ix)%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 ixUpSeg = ixUpSeg_tmp(1:nSegStreamOrder) call indexx(rankSegOrder(ixUpSeg), subSegOrder) - river_basin_out(ix)%branch(iTrib+nOut)%nRch = nSegStreamOrder - river_basin_out(ix)%branch(iTrib+nOut)%segIndex = ixUpSeg(subSegOrder) + branch(iTrib+nOut)%nRch = nSegStreamOrder + branch(iTrib+nOut)%segIndex = ixUpSeg(subSegOrder) + + nSegBranch(iTrib+nOut) = branch(iTrib+nOut)%nRch + end associate end do endif trib + ! re-order branch based on number of reaches + call indexx(nSegBranch, rankBranch) + + do iBrn = nTrib+nOut, 1, -1 + river_basin_out(ix)%branch(iBrn) = branch(rankBranch(iBrn)) + end do + end do sorder - end subroutine omp_domain_decomposition + end subroutine omp_domain_decomposition_stro ! *************************************************************** ! private subroutine: Domain decomposition @@ -459,8 +520,6 @@ subroutine decomposeDomain(structNTOPO, & ! input: integer(i4b), allocatable :: ixSubset(:) ! subset indices based on logical array from global index array logical(lgt), allocatable :: isTribOutlet(:) logical(lgt), allocatable :: isMainstem2d(:,:) - integer(i4b), parameter :: mainstem=2 - integer(i4b), parameter :: tributary=1 ierr=0; message='decomposeDomain/' @@ -570,8 +629,6 @@ subroutine assign_node(nNodes, ierr, message) integer(i4b) :: nTrib ! integer(i4b) :: nEven ! integer(i4b) :: nSmallTrib ! number of small tributaries to be processed in root node - integer(i4b), parameter :: tributary=1 - integer(i4b), parameter :: mainstem=2 logical(lgt),allocatable :: isAssigned(:) ierr=0; message='assign_node/' @@ -640,4 +697,128 @@ subroutine assign_node(nNodes, ierr, message) end subroutine assign_node +! *************************************************************** + ! private subroutine: Assign decomposed domain into procs + ! *************************************************************** + subroutine basin_order(nSeg, structNTOPO_in, domains_omp, nDomain_omp, river_basin_out, ierr, message) + + implicit none + ! Input variables + integer(i4b), intent(in) :: nSeg + type(var_ilength), allocatable, intent(in) :: structNTOPO_in(:) ! network topology + type(subbasin_mpi), intent(in) :: domains_omp(:) ! domain decomposition data structure (maximum domain is set to maxDomain) + integer(i4b) :: nDomain_omp + ! Output variables + type(subbasin_omp),allocatable, intent(out) :: river_basin_out(:)! + integer(i4b), intent(out) :: ierr + character(len=strLen), intent(out) :: message ! error message + ! Local variables + character(len=strLen) :: cmessage ! error message from subroutine + integer(i4b) :: segOrder(nSeg) ! reach order + integer(i4b) :: rankSegOrder(nSeg) ! ranked reach order + integer(i4b) :: nTrib ! number of tributary reaches + integer(i4b) :: nMain ! number of mainstem reaches + integer(i4b) :: ix,ixx, iSeg ! loop indices + integer(i4b),allocatable :: nSubSeg(:) ! + integer(i4b),allocatable :: subSegOrder(:) ! reach order in a subset domain + integer(i4b),allocatable :: rankDomain(:) ! ranked domain based on size + logical(lgt),allocatable :: isAssigned(:) + + ierr=0; message='basin_order/' + + do iSeg=1,nSeg + segOrder(iSeg) = structNTOPO_in(iSeg)%var(ixNTOPO%rchOrder)%dat(1) + enddo + + ! sorting reach processing order + call indexx(segOrder,rankSegOrder) + + allocate(river_basin_out(2), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating [river_basin_out]'; return; endif + + allocate(nSubSeg(nDomain_omp),rankDomain(nDomain_omp),isAssigned(nDomain_omp),stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating [nSubSeg,rankDomain,isAssigned]'; return; endif + + ! rank domains based on number of reaches i.e., nSubSeg - rankDomain + ! count tributaries + + nTrib = 0; nMain = 0 ! initialize number of tributaries and mainstems + do ix = 1,nDomain_omp + nSubSeg(ix) = size(domains_omp(ix)%segIndex) + if (domains_omp(ix)%basinType==tributary) nTrib=nTrib+1 + if (domains_omp(ix)%basinType==mainstem) nMain=nMain+1 + end do + call indexx(nSubSeg, rankDomain) + + if (nTrib/=0) then + allocate(river_basin_out(tributary)%branch(nTrib), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating [river_basin_out%branch]'; return; endif + endif + + if (nMain/=0) then + allocate(river_basin_out(mainstem)%branch(nMain), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating [river_basin_out%branch]'; return; endif + endif + + ! put reaches in tributaries and mainstem in the processing order within domain + nTrib=0; nMain=0 + isAssigned = .false. + domain:do ix = nDomain_omp,1,-1 ! Going through domain from the largest size + + ixx = rankDomain(ix) + + assigned:if (.not. isAssigned(ixx)) then + + associate(ixSegs => domains_omp(ixx)%segIndex) + + ! Compute reach order for only small basin + if (allocated(subSegOrder)) deallocate(subSegOrder) + allocate(subSegOrder(size(ixSegs)), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating subSegOrder'; return; endif + + call indexx(rankSegOrder(ixSegs), subSegOrder) + + domain:if (domains_omp(ixx)%basinType==mainstem) then ! if domain is mainstem + + nMain = nMain + 1 + + allocate(river_basin_out(mainstem)%branch(nMain)%segIndex(size(ixSegs)), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating river_basin_out(mainstem)%branch(1)%segIndex'; return; endif + + river_basin_out(mainstem)%branch(nMain)%segIndex(:) = ixSegs(subSegOrder) + river_basin_out(mainstem)%branch(nMain)%nRch = size(ixSegs) + + isAssigned(ixx) = .true. + + elseif (domains_omp(ixx)%basinType==tributary) then ! if domain is tributary + + nTrib = nTrib + 1 + + allocate(river_basin_out(tributary)%branch(nTrib)%segIndex(size(ixSegs)), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating river_basin_out(tributary)%branch(ix)%segIndex'; return; endif + + river_basin_out(tributary)%branch(nTrib)%segIndex(:) = ixSegs(subSegOrder) + river_basin_out(tributary)%branch(nTrib)%nRch = size(ixSegs) + + isAssigned(ixx) = .true. + + endif domain + + end associate + + endif assigned + + end do domain + + ! check + do ix = 1,nDomain_omp + if (.not.isAssigned(ix)) then + write(cmessage, "(A,I1,A)") 'Domain ', ix, 'is not assigned to any nodes' + ierr = 10; message=trim(message)//trim(cmessage); return + endif + enddo + + end subroutine basin_order + + end module domain_decomposition diff --git a/route/build/src/irf_route.f90 b/route/build/src/irf_route.f90 index c38ef9bb..32a43cef 100644 --- a/route/build/src/irf_route.f90 +++ b/route/build/src/irf_route.f90 @@ -58,11 +58,6 @@ subroutine irf_route(iEns, & ! input: index of runoff ensemble to be p integer(i4b) :: iTrib ! loop indices - branch integer(i4b) :: ix ! loop indices stream order ! variables needed for timing - integer(i4b) :: omp_get_thread_num - integer(i4b), allocatable :: ixThread(:) ! thread id - integer*8, allocatable :: openMPend(:) ! time for the start of the parallelization section - integer*8, allocatable :: timeTribStart(:) ! time Tributaries start - real(dp), allocatable :: timeTrib(:) ! time spent on each Tributary integer*8 :: cr ! rate integer*8 :: startTime,endTime ! date/time for the start and end of the initialization real(dp) :: elapsedTime ! elapsed time for the process @@ -92,17 +87,12 @@ subroutine irf_route(iEns, & ! input: index of runoff ensemble to be p nOrder = size(river_basin) + call system_clock(startTime) + do ix = 1,nOrder nTrib=size(river_basin(ix)%branch) - allocate(ixThread(nTrib), openMPend(nTrib), timeTrib(nTrib), timeTribStart(nTrib), stat=ierr) - if(ierr/=0)then; message=trim(message)//trim(cmessage)//': unable to allocate space for Trib timing'; return; endif - timeTrib(:) = realMissing - ixThread(:) = integerMissing - - call system_clock(startTime) - ! 1. Route tributary reaches (parallel) !$OMP PARALLEL default(none) & !$OMP private(jSeg, iSeg) & ! private for a given thread @@ -112,40 +102,25 @@ subroutine irf_route(iEns, & ! input: index of runoff ensemble to be p !$OMP shared(NETOPO_in) & ! data structure shared !$OMP shared(RCHFLX_out) & ! data structure shared !$OMP shared(ix, iEns, ixDesire) & ! indices shared -!$OMP shared(openMPend, nThreads) & ! timing variables shared -!$OMP shared(timeTribStart) & ! timing variables shared -!$OMP shared(timeTrib) & ! timing variables shared -!$OMP shared(ixThread) & ! thread id array shared !$OMP firstprivate(nTrib) !$OMP DO schedule(dynamic,1) do iTrib = 1,nTrib -!$ ixThread(iTrib) = omp_get_thread_num() - call system_clock(timeTribStart(iTrib)) do iSeg=1,river_basin(ix)%branch(iTrib)%nRch jSeg = river_basin(ix)%branch(iTrib)%segIndex(iSeg) if (.not. doRoute(jSeg)) cycle call segment_irf(iEns, jSeg, ixDesire, NETOPO_IN, RCHFLX_out, ierr, cmessage) ! if(ierr/=0)then; ixmessage(iTrib)=trim(message)//trim(cmessage); exit; endif end do - call system_clock(openMPend(iTrib)) - timeTrib(iTrib) = real(openMPend(iTrib)-timeTribStart(iTrib), kind(dp)) end do !$OMP END DO !$OMP END PARALLEL - call system_clock(endTime) - elapsedTime = real(endTime-startTime, kind(dp))/real(cr) -! write(*,"(A,1PG15.7,A)") ' elapsed-time [routing/irf/tributary] = ', elapsedTime, ' s' - -! write(*,'(a)') 'iTrib nSeg ixThread nThreads StartTime EndTime' -! do iTrib=1,nTrib -! write(*,'(4(i5,1x),2(I20,1x))') iTrib, river_basin(1)%tributary(iTrib)%nRch, ixThread(iTrib), nThreads, timeTribStart(iTrib), openMPend(iTrib) -! enddo - deallocate(ixThread, openMPend, timeTrib, timeTribStart, stat=ierr) - if(ierr/=0)then; message=trim(message)//trim(cmessage)//': unable to deallocate space for Trib timing'; return; endif + end do - end do + call system_clock(endTime) + elapsedTime = real(endTime-startTime, kind(dp))/real(cr) + write(*,"(A,1PG15.7,A)") ' elapsed-time [routing/irf] = ', elapsedTime, ' s' end subroutine irf_route diff --git a/route/build/src/kwt_route.f90 b/route/build/src/kwt_route.f90 index 29b0ad7e..fb666546 100644 --- a/route/build/src/kwt_route.f90 +++ b/route/build/src/kwt_route.f90 @@ -12,7 +12,6 @@ module kwt_route_module USE public_var, only : verySmall ! a very small value USE public_var, only : realMissing ! missing value for real number USE public_var, only : integerMissing ! missing value for integer number -USE globalData, only : nThreads ! number of threads used for openMP ! utilities use nr_utility_module, only : arth ! Num. Recipies utilities USE time_utils_module, only : elapsedSec ! calculate the elapsed time @@ -68,11 +67,6 @@ SUBROUTINE kwt_route(iens, & ! input: ensemble index integer(i4b) :: iTrib ! loop indices - branch integer(i4b) :: ix ! loop indices stream order ! variables needed for timing - integer(i4b) :: omp_get_thread_num - integer(i4b), allocatable :: ixThread(:) ! thread id - integer*8, allocatable :: openMPend(:) ! time for the start of the parallelization section - integer*8, allocatable :: timeTribStart(:) ! time Tributaries start - real(dp), allocatable :: timeTrib(:) ! time spent on each Tributary integer*8 :: cr ! rate integer*8 :: startTime,endTime ! date/time for the start and end of the initialization real(dp) :: elapsedTime ! elapsed time for the process @@ -99,17 +93,12 @@ SUBROUTINE kwt_route(iens, & ! input: ensemble index nOrder = size(river_basin) + call system_clock(startTime) + do ix = 1, nOrder nTrib=size(river_basin(ix)%branch) - allocate(ixThread(nTrib), openMPend(nTrib), timeTrib(nTrib), timeTribStart(nTrib), stat=ierr) - if(ierr/=0)then; message=trim(message)//trim(cmessage)//': unable to allocate space for Trib timing'; return; endif - timeTrib(:) = realMissing - ixThread(:) = integerMissing - - call system_clock(startTime) - ! 1. Route tributary reaches (parallel) !$OMP parallel default(none) & !$OMP private(jSeg, iSeg) & ! private for a given thread @@ -123,16 +112,10 @@ SUBROUTINE kwt_route(iens, & ! input: ensemble index !$OMP shared(KROUTE_out) & ! data structure shared !$OMP shared(RCHFLX_out) & ! data structure shared !$OMP shared(ix, iEns, ixDesire) & ! indices shared -!$OMP shared(openMPend, nThreads) & ! timing variables shared -!$OMP shared(timeTribStart) & ! timing variables shared -!$OMP shared(timeTrib) & ! timing variables shared -!$OMP shared(ixThread) & ! thread id array shared !$OMP firstprivate(nTrib) !$OMP DO schedule(dynamic, 1) ! chunk size of 1 do iTrib = 1,nTrib -!$ ixThread(iTrib) = omp_get_thread_num() - call system_clock(timeTribStart(iTrib)) do iSeg=1,river_basin(ix)%branch(iTrib)%nRch jSeg = river_basin(ix)%branch(iTrib)%segIndex(iSeg) if (.not. doRoute(jSeg)) cycle @@ -148,25 +131,16 @@ SUBROUTINE kwt_route(iens, & ! input: ensemble index ierr,cmessage) ! output: error control !if (ierr/=0)then; message=trim(message)//trim(cmessage); return; endif end do ! (looping through stream segments) - call system_clock(openMPend(iTrib)) - timeTrib(iTrib) = real(openMPend(iTrib)-timeTribStart(iTrib), kind(dp)) end do ! (looping through stream segments) !$OMP END DO !$OMP END PARALLEL - call system_clock(endTime) - elapsedTime = real(endTime-startTime, kind(dp))/real(cr) -! write(*,"(A,1PG15.7,A)") ' elapsed-time [routing/kwt/tributary] = ', elapsedTime, ' s' - -! write(*,'(a)') 'iTrib nSeg ixThread nThreads StartTime EndTime' -! do iTrib=1,nTrib -! write(*,'(4(i5,1x),2(I20,1x))') iTrib, river_basin(1)%tributary(iTrib)%nRch, ixThread(iTrib), nThreads, timeTribStart(iTrib), openMPend(iTrib) -! enddo - deallocate(ixThread, openMPend, timeTrib, timeTribStart, stat=ierr) - if(ierr/=0)then; message=trim(message)//trim(cmessage)//': unable to deallocate space for Trib timing'; return; endif - end do + call system_clock(endTime) + elapsedTime = real(endTime-startTime, kind(dp))/real(cr) + write(*,"(A,1PG15.7,A)") ' elapsed-time [routing/kwt] = ', elapsedTime, ' s' + END SUBROUTINE kwt_route diff --git a/route/build/src/main_route.f90 b/route/build/src/main_route.f90 index 410d948c..e33e95ab 100644 --- a/route/build/src/main_route.f90 +++ b/route/build/src/main_route.f90 @@ -138,6 +138,7 @@ subroutine main_route(& ! perform upstream flow accumulation if (doesAccumRunoff == 1) then call accum_runoff(iens, & ! input: ensemble index + river_basin, & ! input: river basin data type ixPrint, & ! input: index of verbose reach NETOPO_in, & ! input: reach topology data structure RCHFLX_out, & ! inout: reach flux data structure diff --git a/route/build/src/mpi_process.f90 b/route/build/src/mpi_process.f90 index d3513361..d6e374c3 100644 --- a/route/build/src/mpi_process.f90 +++ b/route/build/src/mpi_process.f90 @@ -32,9 +32,9 @@ MODULE mpi_routine USE mpi_mod, ONLY: shr_mpi_allgather implicit none - -integer(i4b),parameter :: scatter=1 ! communication identifier -integer(i4b),parameter :: gather=2 ! communication identifier +! common parameters within this module +integer(i4b),parameter :: scatter=1 +integer(i4b),parameter :: gather=2 private @@ -77,8 +77,9 @@ subroutine comm_ntopo_data(pid, & ! input: proc id USE globalData, ONLY: local_ix_comm ! local reach index at tributary reach outlets to mainstem (size = sum of tributary outlets within entire network) USE alloc_data, ONLY: alloc_struct USE process_ntopo, ONLY: augment_ntopo ! compute all the additional network topology (only compute option = on) - USE process_ntopo, ONLY: put_data_struct ! - USE domain_decomposition,ONLY: omp_domain_decomposition ! domain decomposition for omp + USE process_ntopo, ONLY: put_data_struct ! + USE domain_decomposition,ONLY: omp_domain_decomposition & ! domain decomposition for omp + => omp_domain_decomposition_stro implicit none ! Input variables integer(i4b), intent(in) :: pid ! process id (MPI)