From 30fd99c10da6b3b8e58cbacf14c1ef7a27b6a7aa Mon Sep 17 00:00:00 2001 From: Naoki Mizukami Date: Mon, 26 Jun 2023 05:58:31 -0600 Subject: [PATCH 1/7] New f90 for abstract routing routine for different methods --- route/build/Makefile | 1 + route/build/src/base_route.f90 | 60 ++++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+) create mode 100644 route/build/src/base_route.f90 diff --git a/route/build/Makefile b/route/build/Makefile index d9863b65..5e3d9995 100644 --- a/route/build/Makefile +++ b/route/build/Makefile @@ -153,6 +153,7 @@ DATATYPES = \ datetime_data.f90 \ dataTypes.f90 \ var_lookup.f90 \ + base_route.f90 \ csv_data.f90 \ globalData.f90 \ popMetadat.f90 \ diff --git a/route/build/src/base_route.f90 b/route/build/src/base_route.f90 new file mode 100644 index 00000000..6086625a --- /dev/null +++ b/route/build/src/base_route.f90 @@ -0,0 +1,60 @@ +MODULE base_route + + ! Abstract class for routing method + + implicit none + + private + public:: base_route_rch + public:: routeContainer + + ! --- routing method container + type :: routeContainer + class(base_route_rch), allocatable :: rch_route + end type + + + ! --- routing method container + type, abstract :: base_route_rch + CONTAINS + procedure(sub_route_rch), deferred :: route + end type + + ABSTRACT INTERFACE + + SUBROUTINE sub_route_rch(this, & ! + iEns, & ! input: ensemble index + segIndex, & ! input: reach indice + ixDesire, & ! input: index of verbose reach + T0, T1, & ! input: start and end of simulation time step since start [sec] + NETOPO_in, & ! input: reach topology data structure + RPARAM_in, & ! input: reach parameter data structure + RCHSTA_out, & ! inout: reach state data structure + RCHFLX_out, & ! inout: reach flux data structure + ierr, message) ! output: error control + + USE nrtype + USE dataTypes, ONLY: STRFLX ! fluxes in each reach + USE dataTypes, ONLY: STRSTA ! states in each reach + USE dataTypes, ONLY: RCHTOPO ! Network topology + USE dataTypes, ONLY: RCHPRP ! Reach parameter + + import base_route_rch + ! Arguments + class(base_route_rch) :: this + integer(i4b), intent(in) :: iEns ! ensemble member + integer(i4b), intent(in) :: segIndex ! ensemble member + integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output + real(dp), intent(in) :: T0, T1 ! start and end of the time step (seconds) + type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology + type(RCHPRP), intent(inout), allocatable :: RPARAM_in(:) ! River reach parameter + type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data + type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + + END SUBROUTINE sub_route_rch + + END INTERFACE + +END MODULE base_route From f0113dec74b622ef15b45d3e5f738b5327417710 Mon Sep 17 00:00:00 2001 From: Naoki Mizukami Date: Mon, 26 Jun 2023 06:11:35 -0600 Subject: [PATCH 2/7] added routing method container to hold routing method objects --- route/build/src/globalData.f90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/route/build/src/globalData.f90 b/route/build/src/globalData.f90 index 6123f530..805b9e01 100644 --- a/route/build/src/globalData.f90 +++ b/route/build/src/globalData.f90 @@ -38,6 +38,8 @@ MODULE globalData USE datetime_data, ONLY: datetime ! datetime data class + USE base_route, ONLY: routeContainer ! a container of instantiated routing methods + USE var_lookup, ONLY: nStructures ! number of variables in data structure (struct_info) USE var_lookup, ONLY: nDimensions ! number of variables in dimensions related to network topology USE var_lookup, ONLY: nStateDims ! number of variables in dimensions related to restart variables @@ -80,6 +82,7 @@ MODULE globalData integer(i4b), allocatable, public :: reachID(:) ! reach id in the whole river network ! ---------- routing methods ------------------------------------------------------------------------- + type(routeContainer), allocatable , public :: rch_routes(:) ! a collection of routing method objects integer(i4b) , public :: nRoutes ! number of active routing methods integer(i4b) , allocatable , public :: routeMethods(:) ! active routing method id logical(lgt) , public :: onRoute(0:nRouteMethods-1) ! logical to indicate active routing method(s) From 580305386fb2204ea8c40f32671914e1ba140c97 Mon Sep 17 00:00:00 2001 From: Naoki Mizukami Date: Mon, 26 Jun 2023 06:13:57 -0600 Subject: [PATCH 3/7] define specific each routing method --- route/build/src/accum_runoff.f90 | 175 +++++------------- route/build/src/dfw_route.f90 | 145 +++------------ route/build/src/irf_route.f90 | 152 ++++------------ route/build/src/kwe_route.f90 | 143 +++------------ route/build/src/kwt_route.f90 | 303 ++++++++++--------------------- route/build/src/mc_route.f90 | 140 +++----------- 6 files changed, 245 insertions(+), 813 deletions(-) diff --git a/route/build/src/accum_runoff.f90 b/route/build/src/accum_runoff.f90 index cf080dbd..ac3d53d2 100644 --- a/route/build/src/accum_runoff.f90 +++ b/route/build/src/accum_runoff.f90 @@ -1,147 +1,64 @@ -MODULE accum_runoff_module +MODULE accum_route_module ! Accumulate upstream flow instantaneously +! This is not used as routed runoff. +! This is used to get total instantaneous upstream runoff at each reach USE nrtype -! data type -USE dataTypes, ONLY: STRFLX ! fluxes in each reach -USE dataTypes, ONLY: RCHTOPO ! Network topology -USE dataTypes, ONLY: subbasin_omp ! mainstem+tributary data structures -! global data -USE public_var, ONLY: iulog ! i/o logical unit number -USE globalData, ONLY: idxSUM ! index of accumulation method -! subroutines: general -USE perf_mod, ONLY: t_startf,t_stopf ! timing start/stop -USE model_utils, ONLY: handle_err +USE dataTypes, ONLY: STRFLX ! fluxes in each reach +USE dataTypes, ONLY: STRSTA ! state in each reach +USE dataTypes, ONLY: RCHTOPO ! Network topology +USE dataTypes, ONLY: RCHPRP ! Reach parameter +USE public_var, ONLY: iulog ! i/o logical unit number +USE globalData, ONLY: idxSUM ! index of accumulation method +USE base_route, ONLY: base_route_rch implicit none private -public::accum_runoff +public::accum_route_rch -CONTAINS - - ! --------------------------------------------------------------------------------------- - ! Public subroutine main driver for basin routing - ! --------------------------------------------------------------------------------------- - 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 - ierr, message, & ! output: error controls - ixSubRch) ! optional input: subset of reach indices to be processed - ! ---------------------------------------------------------------------------------------- - ! Purpose: - ! - ! Accumulate all the upstream delayed runoff for each reach - ! This is not used as routed runoff. - ! This is used to get total instantaneous upstream runoff at each reach - ! - ! ---------------------------------------------------------------------------------------- - - implicit none - ! argument variables - integer(i4b), intent(in) :: iens ! runoff ensemble index - type(subbasin_omp), allocatable, intent(in) :: river_basin(:) ! river basin information (mainstem, tributary outlet etc.) - integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output - type(RCHTOPO), allocatable, intent(in) :: NETOPO_in(:) ! River Network topology - type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message ! error message - integer(i4b), optional, intent(in) :: ixSubRch(:) ! subset of reach indices to be processed - ! local variables - 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 - - ierr=0; message='accum_runoff/' - - if (size(NETOPO_in)/=size(RCHFLX_out(iens,:))) then - ierr=20; message=trim(message)//'sizes of NETOPO and RCHFLX mismatch'; return - endif - - nSeg = size(RCHFLX_out(iens,:)) - - allocate(doRoute(nSeg), stat=ierr) - if(ierr/=0)then; message=trim(message)//'unable to allocate space for [doRoute]'; return; endif - - ! if a subset of reaches is processed - if (present(ixSubRch))then - doRoute(:)=.false. - doRoute(ixSubRch) = .true. ! only subset of reaches are on - ! if all the reaches are processed - else - doRoute(:)=.true. ! every reach is on - endif - - nDom = size(river_basin) +type, extends(base_route_rch) :: accum_route_rch + CONTAINS + procedure, pass :: route => accum_inst_runoff +end type accum_route_rch - call t_startf('route/accum-runoff') - - 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) call handle_err(ierr, trim(message)//trim(cmessage)) - - end do - end do -!$OMP END PARALLEL DO - - end do - - call t_stopf('route/accum-runoff') - - END SUBROUTINE accum_runoff +CONTAINS ! ********************************************************************* ! subroutine: perform accumulate immediate upstream flow ! ********************************************************************* - SUBROUTINE accum_qupstream(iEns, & ! input: index of runoff ensemble to be processed - segIndex, & ! input: index of reach to be processed - ixDesire, & ! input: reachID to be checked by on-screen pringing - NETOPO_in, & ! input: reach topology data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr, message) ! output: error control + SUBROUTINE accum_inst_runoff(this, & + iEns, & ! input: index of runoff ensemble to be processed + segIndex, & ! input: index of reach to be processed + ixDesire, & ! input: reachID to be checked by on-screen pringing + T0,T1, & ! input: start and end of the time step + NETOPO_in, & ! input: reach topology data structure + RPARAM_in, & ! input: reach parameter data structure + RCHSTA_out, & ! inout: reach state data structure + RCHFLX_out, & ! inout: reach flux data structure + ierr, message) ! output: error control implicit none - ! argument variables - integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed - integer(i4b), intent(in) :: segIndex ! segment where routing is performed - integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output - type(RCHTOPO),intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(STRFLX), intent(inout) :: RCHFLX_out(:,:)! Reach fluxes (ensembles, space [reaches]) for decomposed domains - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message ! error message + ! Argument variables + class(accum_route_rch) :: this + integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed + integer(i4b), intent(in) :: segIndex ! segment where routing is performed + integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output + real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) + type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology + type(RCHPRP), intent(inout),allocatable :: RPARAM_in(:) ! River reach parameter + type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data + type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message ! Local variables - real(dp) :: q_upstream ! upstream Reach fluxes - integer(i4b) :: nUps ! number of upstream segment - integer(i4b) :: iUps ! upstream reach index - integer(i4b) :: iRch_ups ! index of upstream reach in NETOPO - character(len=strLen) :: fmt1,fmt2 ! format string + real(dp) :: q_upstream ! upstream Reach fluxes + integer(i4b) :: nUps ! number of upstream segment + integer(i4b) :: iUps ! upstream reach index + integer(i4b) :: iRch_ups ! index of upstream reach in NETOPO + character(len=strLen) :: fmt1,fmt2 ! format string - ierr=0; message='accum_qupstream/' + ierr=0; message='accum_inst_runoff/' ! identify number of upstream segments of the reach being processed nUps = size(NETOPO_in(segIndex)%UREACHI) @@ -176,6 +93,6 @@ SUBROUTINE accum_qupstream(iEns, & ! input: index of runoff ensemble to write(iulog,'(a,x,G15.4)') ' RCHFLX_out%ROUTE(idxSUM)%REACH_Q =', RCHFLX_out(iens,segIndex)%ROUTE(idxSUM)%REACH_Q endif - END SUBROUTINE accum_qupstream + END SUBROUTINE accum_inst_runoff -END MODULE accum_runoff_module +END MODULE accum_route_module diff --git a/route/build/src/dfw_route.f90 b/route/build/src/dfw_route.f90 index 170148f0..91e388d4 100644 --- a/route/build/src/dfw_route.f90 +++ b/route/build/src/dfw_route.f90 @@ -1,135 +1,38 @@ MODULE dfw_route_module +! diffusive wave routing + USE nrtype -! data types -USE dataTypes, ONLY: STRFLX ! fluxes in each reach -USE dataTypes, ONLY: STRSTA ! state in each reach -USE dataTypes, ONLY: RCHTOPO ! Network topology -USE dataTypes, ONLY: RCHPRP ! Reach parameter -USE dataTypes, ONLY: dwRCH ! dw specific state data structure -USE dataTypes, ONLY: subbasin_omp ! mainstem+tributary data strucuture -! global data -USE public_var, ONLY: iulog ! i/o logical unit number -USE public_var, ONLY: realMissing ! missing value for real number -USE public_var, ONLY: integerMissing ! missing value for integer number -USE public_var, ONLY: qmodOption ! qmod option (use 1==direct insertion) -USE globalData, ONLY: idxDW -! subroutines: general -USE water_balance, ONLY: comp_reach_wb ! compute water balance error -USE perf_mod, ONLY: t_startf,t_stopf ! timing start/stop -USE model_utils, ONLY: handle_err +USE dataTypes, ONLY: STRFLX ! fluxes in each reach +USE dataTypes, ONLY: STRSTA ! state in each reach +USE dataTypes, ONLY: RCHTOPO ! Network topology +USE dataTypes, ONLY: RCHPRP ! Reach parameter +USE dataTypes, ONLY: dwRCH ! dw specific state data structure +USE public_var, ONLY: iulog ! i/o logical unit number +USE public_var, ONLY: realMissing ! missing value for real number +USE public_var, ONLY: integerMissing ! missing value for integer number +USE public_var, ONLY: qmodOption ! qmod option (use 1==direct insertion) +USE globalData, ONLY: idxDW +USE water_balance, ONLY: comp_reach_wb ! compute water balance error +USE base_route, ONLY: base_route_rch implicit none private -public::dfw_route - -CONTAINS - - ! ********************************************************************* - ! subroutine: perform diffusive wave routing through the river network - ! ********************************************************************* - SUBROUTINE dfw_route(iens, & ! input: ensemble index - river_basin, & ! input: river basin information (mainstem, tributary outlet etc.) - T0,T1, & ! input: start and end of the time step - ixDesire, & ! input: reachID to be checked by on-screen pringing - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,message, & ! output: error control - ixSubRch) ! optional input: subset of reach indices to be processed - - USE public_var, ONLY: is_lake_sim ! logical whether or not lake should be simulated - - implicit none - ! Argument variables - integer(i4b), intent(in) :: iEns ! ensemble member - type(subbasin_omp), intent(in), allocatable :: river_basin(:) ! river basin information (mainstem, tributary outlet etc.) - real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) - integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output - type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter - type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data - type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message ! error message - integer(i4b), intent(in), optional :: ixSubRch(:) ! subset of reach indices to be processed - ! local variables - character(len=strLen) :: cmessage ! error message for downwind routine - logical(lgt), allocatable :: doRoute(:) ! logical to indicate which reaches are processed - integer(i4b) :: nOrder ! number of stream order - integer(i4b) :: nTrib ! number of tributary basins - integer(i4b) :: nSeg ! number of reaches in the network - integer(i4b) :: iSeg, jSeg ! loop indices - reach - integer(i4b) :: iTrib ! loop indices - branch - integer(i4b) :: ix ! loop indices stream order - - ierr=0; message='dfw_route/' - - ! number of reach check - if (size(NETOPO_in)/=size(RCHFLX_out(iens,:))) then - ierr=20; message=trim(message)//'sizes of NETOPO and RCHFLX mismatch'; return - endif +public::dfw_route_rch - nSeg = size(RCHFLX_out(iens,:)) +type, extends(base_route_rch) :: dfw_route_rch + CONTAINS + procedure, pass :: route => dfw_rch +end type dfw_route_rch - allocate(doRoute(nSeg), stat=ierr) - - if (present(ixSubRch))then - doRoute(:)=.false. - doRoute(ixSubRch) = .true. ! only subset of reaches are on - else - doRoute(:)=.true. ! every reach is on - endif - - nOrder = size(river_basin) - - call t_startf('route/dfw') - - do ix = 1, nOrder - - nTrib=size(river_basin(ix)%branch) - -!$OMP PARALLEL DO schedule(dynamic,1) & ! chunk size of 1 -!$OMP private(jSeg, iSeg) & ! private for a given thread -!$OMP private(ierr, cmessage) & ! private for a given thread -!$OMP shared(T0,T1) & ! 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(RPARAM_in) & ! data structure shared -!$OMP shared(RCHSTA_out) & ! data structure shared -!$OMP shared(RCHFLX_out) & ! data structure shared -!$OMP shared(ix, iEns, ixDesire) & ! indices shared -!$OMP firstprivate(nTrib) - trib:do iTrib = 1,nTrib - seg:do iSeg=1,river_basin(ix)%branch(iTrib)%nRch - jSeg = river_basin(ix)%branch(iTrib)%segIndex(iSeg) - if (.not. doRoute(jSeg)) cycle - call dfw_rch(iEns,jSeg, & ! input: array indices - ixDesire, & ! input: index of the desired reach - T0,T1, & ! input: start and end of the time step - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,cmessage) ! output: error control - if(ierr/=0) call handle_err(ierr, trim(message)//trim(cmessage)) - end do seg - end do trib -!$OMP END PARALLEL DO - - end do - - call t_stopf('route/dfw') - - END SUBROUTINE dfw_route +CONTAINS ! ********************************************************************* ! subroutine: perform diffusive wave routing for one segment ! ********************************************************************* - SUBROUTINE dfw_rch(iEns, segIndex, & ! input: index of runoff ensemble to be processed + SUBROUTINE dfw_rch(this, & + iEns, segIndex, & ! input: index of runoff ensemble to be processed ixDesire, & ! input: reachID to be checked by on-screen pringing T0,T1, & ! input: start and end of the time step NETOPO_in, & ! input: reach topology data structure @@ -140,12 +43,13 @@ SUBROUTINE dfw_rch(iEns, segIndex, & ! input: index of runoff ensemble to be pro implicit none ! Argument variables + class(dfw_route_rch) :: this integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed integer(i4b), intent(in) :: segIndex ! segment where routing is performed integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter + type(RCHPRP), intent(inout), allocatable :: RPARAM_in(:) ! River reach parameter type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains integer(i4b), intent(out) :: ierr ! error code @@ -550,5 +454,4 @@ SUBROUTINE TDMA(NX,MAT,b,T) END SUBROUTINE TDMA - END MODULE dfw_route_module diff --git a/route/build/src/irf_route.f90 b/route/build/src/irf_route.f90 index 1ecfa797..7cca3e3d 100644 --- a/route/build/src/irf_route.f90 +++ b/route/build/src/irf_route.f90 @@ -1,130 +1,42 @@ MODULE irf_route_module USE nrtype -! data type -USE dataTypes, ONLY: STRFLX ! fluxes in each reach -USE dataTypes, ONLY: RCHTOPO ! Network topology -USE dataTypes, ONLY: RCHPRP ! reach/lake property parameter -USE dataTypes, ONLY: irfRCH ! irf specific state data structure -USE dataTypes, ONLY: subbasin_omp ! mainstem+tributary data structures -! global parameters -USE public_var, ONLY: iulog ! i/o logical unit number -USE public_var, ONLY: realMissing ! missing value for real number -USE public_var, ONLY: integerMissing ! missing value for integer number -USE public_var, ONLY: dt -USE globalData, ONLY: idxIRF ! index of IRF method -! external subroutines -USE lake_route_module, ONLY: lake_route ! lake route module +USE dataTypes, ONLY: STRFLX ! fluxes in each reach +USE dataTypes, ONLY: STRSTA ! state in each reach +USE dataTypes, ONLY: RCHTOPO ! Network topology +USE dataTypes, ONLY: RCHPRP ! reach/lake property parameter +USE dataTypes, ONLY: irfRCH ! irf specific state data structure +USE public_var, ONLY: iulog ! i/o logical unit number +USE public_var, ONLY: realMissing ! missing value for real number +USE public_var, ONLY: integerMissing ! missing value for integer number +USE public_var, ONLY: dt +USE globalData, ONLY: idxIRF ! index of IRF method USE water_balance, ONLY: comp_reach_wb ! compute water balance error -USE perf_mod, ONLY: t_startf,t_stopf ! timing start/stop -USE model_utils, ONLY: handle_err +USE base_route, ONLY: base_route_rch implicit none private -public::irf_route +public::irf_route_rch -CONTAINS - - ! ********************************************************************* - ! subroutine: perform network UH routing - ! ********************************************************************* - SUBROUTINE irf_route(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 pringing - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr, message, & ! output: error control - ixSubRch) ! optional input: subset of reach indices to be processed +type, extends(base_route_rch) :: irf_route_rch + CONTAINS + procedure, pass :: route => irf_rch +end type irf_route_rch - USE public_var, ONLY: is_lake_sim ! logical whether or not lake should be simulated - - implicit none - ! argument variables - integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed - 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 ! Output - type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(RCHPRP), intent(inout), allocatable :: RPARAM_in(:) ! River Network parameters - TYPE(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message ! error message - integer(i4b), intent(in), optional :: ixSubRch(:) ! subset of reach indices to be processed - ! Local variables - character(len=strLen) :: cmessage ! error message from subroutine - logical(lgt), allocatable :: doRoute(:) ! logical to indicate which reaches are processed - integer(i4b) :: nOrder ! number of stream order - integer(i4b) :: nTrib ! number of tributary basins - integer(i4b) :: nSeg ! number of reaches in the network - integer(i4b) :: iSeg, jSeg ! loop indices - reach - integer(i4b) :: iTrib ! loop indices - branch - integer(i4b) :: ix ! loop indices stream order - - ierr=0; message='irf_route/' - - ! number of reach check - if (size(NETOPO_in)/=size(RCHFLX_out(iens,:))) then - ierr=20; message=trim(message)//'sizes of NETOPO and RCHFLX mismatch'; return - endif - - nSeg = size(RCHFLX_out(iens,:)) - - allocate(doRoute(nSeg), stat=ierr) - - if (present(ixSubRch))then - doRoute(:)=.false. - doRoute(ixSubRch) = .true. ! only subset of reaches are on - else - doRoute(:)=.true. ! every reach is on - endif - - nOrder = size(river_basin) - - call t_startf('route/irf') - - do ix = 1,nOrder - - nTrib=size(river_basin(ix)%branch) - - ! 1. Route tributary reaches (parallel) -!$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(RPARAM_in) & ! data structure shared -!$OMP shared(RCHFLX_out) & ! data structure shared -!$OMP shared(ix, iEns, ixDesire) & ! indices shared -!$OMP firstprivate(nTrib) - trib:do iTrib = 1,nTrib - seg:do iSeg=1,river_basin(ix)%branch(iTrib)%nRch - jSeg = river_basin(ix)%branch(iTrib)%segIndex(iSeg) - if (.not. doRoute(jSeg)) cycle - if ((NETOPO_in(jseg)%islake).and.(is_lake_sim)) then - call lake_route (iEns, jSeg, ixDesire, NETOPO_in, RPARAM_in, RCHFLX_out, ierr, cmessage) - else - call irf_rch(iEns, jSeg, ixDesire, NETOPO_in, RCHFLX_out, ierr, cmessage) - endif - if(ierr/=0) call handle_err(ierr, trim(message)//trim(cmessage)) - end do seg - end do trib -!$OMP END PARALLEL DO - - end do - - call t_stopf('route/irf') - - end subroutine irf_route +CONTAINS ! ********************************************************************* ! subroutine: perform one segment route UH routing ! ********************************************************************* - SUBROUTINE irf_rch(iEns, & ! input: index of runoff ensemble to be processed - segIndex, & ! input: index of runoff ensemble to be processed + SUBROUTINE irf_rch(this, & ! + iEns, & ! input: index of runoff ensemble + segIndex, & ! input: reach index ixDesire, & ! input: reachID to be checked by on-screen pringing + T0,T1, & ! input: start and end of the time step NETOPO_in, & ! input: reach topology data structure + RPARAM_in, & ! input: reach parameter data structure + RCHSTA_out, & ! inout: reach state data structure RCHFLX_out, & ! inout: reach flux data structure ierr, message) ! output: error control @@ -132,13 +44,17 @@ SUBROUTINE irf_rch(iEns, & ! input: index of runoff ensemble to be proce implicit none ! Argument variables - integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed - integer(i4b), intent(in) :: segIndex ! segment where routing is performed - integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output - type(RCHTOPO),intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(STRFLX), intent(inout) :: RCHFLX_out(:,:)! Reach fluxes (ensembles, space [reaches]) for decomposed domains - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message ! error message + class(irf_route_rch) :: this + integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed + integer(i4b), intent(in) :: segIndex ! segment where routing is performed + integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output + real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) + type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology + type(RCHPRP), intent(inout),allocatable :: RPARAM_in(:) ! River reach parameter + type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data + type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message ! Local variables logical(lgt) :: verbose ! check details of variables real(dp) :: q_upstream ! total discharge at top of the reach being processed diff --git a/route/build/src/kwe_route.f90 b/route/build/src/kwe_route.f90 index d38a2c79..6117d95b 100644 --- a/route/build/src/kwe_route.f90 +++ b/route/build/src/kwe_route.f90 @@ -3,138 +3,38 @@ MODULE kw_route_module ! kinematic wave routing USE nrtype -! data types -USE dataTypes, ONLY: STRFLX ! fluxes in each reach -USE dataTypes, ONLY: STRSTA ! state in each reach -USE dataTypes, ONLY: RCHTOPO ! Network topology -USE dataTypes, ONLY: RCHPRP ! Reach parameter -USE dataTypes, ONLY: kwRCH ! kw specific state data structure -USE dataTypes, ONLY: subbasin_omp ! mainstem+tributary data strucuture -! global data -USE public_var, ONLY: iulog ! i/o logical unit number -USE public_var, ONLY: realMissing ! missing value for real number -USE public_var, ONLY: integerMissing ! missing value for integer number -USE public_var, ONLY: qmodOption ! qmod option (use 1==direct insertion) -USE globalData, ONLY: idxKW ! loop index of kw routing -! subroutines: general -USE water_balance, ONLY: comp_reach_wb ! compute water balance error -USE perf_mod, ONLY: t_startf,t_stopf ! timing start/stop -USE model_utils, ONLY: handle_err +USE dataTypes, ONLY: STRFLX ! fluxes in each reach +USE dataTypes, ONLY: STRSTA ! state in each reach +USE dataTypes, ONLY: RCHTOPO ! Network topology +USE dataTypes, ONLY: RCHPRP ! Reach parameter +USE dataTypes, ONLY: kwRCH ! kw specific state data structure +USE public_var, ONLY: iulog ! i/o logical unit number +USE public_var, ONLY: realMissing ! missing value for real number +USE public_var, ONLY: integerMissing ! missing value for integer number +USE public_var, ONLY: qmodOption ! qmod option (use 1==direct insertion) +USE globalData, ONLY: idxKW ! loop index of kw routing +USE water_balance, ONLY: comp_reach_wb ! compute water balance error +USE base_route, ONLY: base_route_rch implicit none private -public::kw_route +public::kwe_route_rch real(dp), parameter :: critFactor=0.01 -CONTAINS - - ! ********************************************************************* - ! subroutine: route kinematic waves with Euler solution through the river network - ! ********************************************************************* - SUBROUTINE kw_route(iens, & ! input: ensemble index - river_basin, & ! input: river basin information (mainstem, tributary outlet etc.) - T0,T1, & ! input: start and end of the time step - ixDesire, & ! input: reachID to be checked by on-screen pringing - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,message, & ! output: error control - ixSubRch) ! optional input: subset of reach indices to be processed - - USE public_var, ONLY: is_lake_sim ! logical whether or not lake should be simulated - - implicit none - ! Argument variables - integer(i4b), intent(in) :: iEns ! ensemble member - type(subbasin_omp), intent(in), allocatable :: river_basin(:) ! river basin information (mainstem, tributary outlet etc.) - real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) - integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output - type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter - type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data - type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message ! error message - integer(i4b), intent(in), optional :: ixSubRch(:) ! subset of reach indices to be processed - ! local variables - character(len=strLen) :: cmessage ! error message for downwind routine - logical(lgt), allocatable :: doRoute(:) ! logical to indicate which reaches are processed - integer(i4b) :: nOrder ! number of stream order - integer(i4b) :: nTrib ! number of tributary basins - integer(i4b) :: nSeg ! number of reaches in the network - integer(i4b) :: iSeg, jSeg ! loop indices - reach - integer(i4b) :: iTrib ! loop indices - branch - integer(i4b) :: ix ! loop indices stream order - - ierr=0; message='kw_route/' - - ! number of reach check - if (size(NETOPO_in)/=size(RCHFLX_out(iens,:))) then - ierr=20; message=trim(message)//'sizes of NETOPO and RCHFLX mismatch'; return - endif - - nSeg = size(RCHFLX_out(iens,:)) - - allocate(doRoute(nSeg), stat=ierr) +type, extends(base_route_rch) :: kwe_route_rch + CONTAINS + procedure, pass :: route => kw_rch +end type kwe_route_rch - if (present(ixSubRch))then - doRoute(:)=.false. - doRoute(ixSubRch) = .true. ! only subset of reaches are on - else - doRoute(:)=.true. ! every reach is on - endif - - nOrder = size(river_basin) - - call t_startf('route/kw') - - ! route kinematic waves through the river network - do ix = 1, nOrder - - nTrib=size(river_basin(ix)%branch) - -!$OMP PARALLEL DO schedule(dynamic,1) & ! chunk size of 1 -!$OMP private(jSeg, iSeg) & ! private for a given thread -!$OMP private(ierr, cmessage) & ! private for a given thread -!$OMP shared(T0,T1) & ! 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(RPARAM_in) & ! data structure shared -!$OMP shared(RCHSTA_out) & ! data structure shared -!$OMP shared(RCHFLX_out) & ! data structure shared -!$OMP shared(ix, iEns, ixDesire) & ! indices shared -!$OMP firstprivate(nTrib) - trib:do iTrib = 1,nTrib - seg:do iSeg=1,river_basin(ix)%branch(iTrib)%nRch - jSeg = river_basin(ix)%branch(iTrib)%segIndex(iSeg) - if (.not. doRoute(jSeg)) cycle - call kw_rch(iEns,jSeg, & ! input: array indices - ixDesire, & ! input: index of the desired reach - T0,T1, & ! input: start and end of the time step - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,cmessage) ! output: error control - if(ierr/=0) call handle_err(ierr, trim(message)//trim(cmessage)) - end do seg - end do trib -!$OMP END PARALLEL DO - - end do - - call t_stopf('route/kw') - - END SUBROUTINE kw_route +CONTAINS ! ********************************************************************* ! subroutine: perform one segment route KW routing ! ********************************************************************* - SUBROUTINE kw_rch(iEns, segIndex, & ! input: index of runoff ensemble to be processed + SUBROUTINE kw_rch(this, & + iEns, segIndex, & ! input: index of runoff ensemble to be processed ixDesire, & ! input: reachID to be checked by on-screen pringing T0,T1, & ! input: start and end of the time step NETOPO_in, & ! input: reach topology data structure @@ -144,12 +44,13 @@ SUBROUTINE kw_rch(iEns, segIndex, & ! input: index of runoff ensemble to be proc ierr, message) ! output: error control implicit none ! Argument variables + class(kwe_route_rch) :: this integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed integer(i4b), intent(in) :: segIndex ! segment where routing is performed integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter + type(RCHPRP), intent(inout), allocatable :: RPARAM_in(:) ! River reach parameter type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains integer(i4b), intent(out) :: ierr ! error code diff --git a/route/build/src/kwt_route.f90 b/route/build/src/kwt_route.f90 index 52a09307..3ee0e0e3 100644 --- a/route/build/src/kwt_route.f90 +++ b/route/build/src/kwt_route.f90 @@ -1,150 +1,48 @@ MODULE kwt_route_module USE nrtype -! data types -USE dataTypes, ONLY: FPOINT ! particle -USE dataTypes, ONLY: STRFLX ! fluxes in each reach -USE dataTypes, ONLY: STRSTA ! states in each reach -USE dataTypes, ONLY: RCHTOPO ! Network topology -USE dataTypes, ONLY: RCHPRP ! Reach parameter -USE dataTypes, ONLY: kwtRCH ! kwt specific state data structure -USE dataTypes, ONLY: subbasin_omp ! mainstem+tributary data strucuture -! global data -USE public_var, ONLY: iulog ! i/o logical unit number -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: idxKWT ! index of KWT method -! utilities -USE nr_utils, ONLY: arth ! Num. Recipies utilities -! subroutines: general -USE perf_mod, ONLY: t_startf,t_stopf ! timing start/stop -USE model_utils, ONLY: handle_err +USE dataTypes, ONLY: FPOINT ! particle +USE dataTypes, ONLY: STRFLX ! fluxes in each reach +USE dataTypes, ONLY: STRSTA ! states in each reach +USE dataTypes, ONLY: RCHTOPO ! Network topology +USE dataTypes, ONLY: RCHPRP ! Reach parameter +USE dataTypes, ONLY: kwtRCH ! kwt specific state data structure +USE public_var, ONLY: iulog ! i/o logical unit number +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: idxKWT ! index of KWT method +USE nr_utils, ONLY: arth ! Num. Recipies utilities +USE base_route, ONLY: base_route_rch implicit none private -public::kwt_route +public::kwt_route_rch -CONTAINS - - ! ********************************************************************* - ! subroutine: route kinematic waves through the river network - ! ********************************************************************* - SUBROUTINE kwt_route(iens, & ! input: ensemble index - river_basin, & ! input: river basin information (mainstem, tributary outlet etc.) - T0,T1, & ! input: start and end of the time step - ixDesire, & ! input: reachID to be checked by on-screen pringing - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,message, & ! output: error control - ixSubRch) ! optional input: subset of reach indices to be processed +integer(i4b) :: LAKEFLAG=0 ! >0 if processing lakes +integer(i4b) :: RSTEP=0 ! retrospective time step offset (what is this for?? used for optional input) - implicit none - ! Argument variables - integer(i4b), intent(in) :: iEns ! ensemble member - type(subbasin_omp), intent(in), allocatable :: river_basin(:) ! river basin information (mainstem, tributary outlet etc.) - real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) - integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output - type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter - type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data - type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message ! error message - integer(i4b), intent(in), optional :: ixSubRch(:) ! subset of reach indices to be processed - ! local variables - character(len=strLen) :: cmessage ! error message for downwind routine - logical(lgt), allocatable :: doRoute(:) ! logical to indicate which reaches are processed - integer(i4b) :: LAKEFLAG=0 ! >0 if processing lakes - integer(i4b) :: nOrder ! number of stream order - integer(i4b) :: nTrib ! number of tributary basins - integer(i4b) :: nSeg ! number of reaches in the network - integer(i4b) :: iSeg, jSeg ! loop indices - reach - integer(i4b) :: iTrib ! loop indices - branch - integer(i4b) :: ix ! loop indices stream order - - ierr=0; message='kwt_route/' - - ! number of reach check - if (size(NETOPO_in)/=size(RCHFLX_out(iens,:))) then - ierr=20; message=trim(message)//'sizes of NETOPO and RCHFLX mismatch'; return - endif - - nSeg = size(RCHFLX_out(iens,:)) - - allocate(doRoute(nSeg), stat=ierr) - if(ierr/=0)then; message=trim(message)//'problem allocating space for [doRoute]'; return; endif - - if (present(ixSubRch))then - doRoute(:)=.false. - doRoute(ixSubRch) = .true. ! only subset of reaches are on - else - doRoute(:)=.true. ! every reach is on - endif - - nOrder = size(river_basin) - - call t_startf('route/kwt') - - ! route kinematic waves through the river network - do ix = 1, nOrder - - nTrib=size(river_basin(ix)%branch) - -!$OMP PARALLEL DO schedule(dynamic,1) & ! chunk size of 1 -!$OMP private(jSeg, iSeg) & ! private for a given thread -!$OMP private(ierr, cmessage) & ! private for a given thread -!$OMP shared(T0,T1) & ! private for a given thread -!$OMP shared(LAKEFLAG) & ! 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(RPARAM_in) & ! data structure shared -!$OMP shared(RCHSTA_out) & ! data structure shared -!$OMP shared(RCHFLX_out) & ! data structure shared -!$OMP shared(ix, iEns, ixDesire) & ! indices shared -!$OMP firstprivate(nTrib) - trib:do iTrib = 1,nTrib - seg:do iSeg=1,river_basin(ix)%branch(iTrib)%nRch - jSeg = river_basin(ix)%branch(iTrib)%segIndex(iSeg) - if (.not. doRoute(jSeg)) cycle - call qroute_rch(iEns,jSeg, & ! input: array indices - ixDesire, & ! input: index of verbose reach - T0,T1, & ! input: start and end of the time step - LAKEFLAG, & ! input: flag if lakes are to be processed - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,cmessage) ! output: error control - if(ierr/=0) call handle_err(ierr, trim(message)//trim(cmessage)) - end do seg - end do trib -!$OMP END PARALLEL DO - - end do ! basin loop - - call t_stopf('route/kwt') - - END SUBROUTINE kwt_route +type, extends(base_route_rch) :: kwt_route_rch + CONTAINS + procedure, pass :: route => kwt_rch +end type kwt_route_rch +CONTAINS ! ********************************************************************* ! subroutine: route kinematic waves at one segment ! ********************************************************************* - SUBROUTINE qroute_rch(IENS,JRCH, & ! input: array indices - ixDesire, & ! input: index of the reach for verbose output - T0,T1, & ! input: start and end of the time step - LAKEFLAG, & ! input: flag if lakes are to be processed - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,message, & ! output: error control - RSTEP) ! optional input: retrospective time step offset + SUBROUTINE kwt_rch(this, & + iEns, & ! input: ensemble index + segIndex, & ! input: index of reach to be processed + ixDesire, & ! input: index of the reach for verbose output + T0,T1, & ! input: start and end of the time step + NETOPO_in, & ! input: reach topology data structure + RPARAM_in, & ! input: reach parameter data structure + RCHSTA_out, & ! inout: reach state data structure + RCHFLX_out, & ! inout: reach flux data structure + ierr,message) ! output: error control USE public_var, ONLY: MAXQPAR ! maximum number of waves per reach USE public_var, ONLY: is_flux_wm ! logical whether or not water management (abstract/injection) is on ! ---------------------------------------------------------------------------------------- @@ -198,14 +96,13 @@ SUBROUTINE qroute_rch(IENS,JRCH, & ! input: array indices ! ---------------------------------------------------------------------------------------- implicit none ! Argument variables - integer(i4b), intent(in) :: IENS ! ensemble member - integer(i4b), intent(in) :: JRCH ! reach to process + class(kwt_route_rch) :: this + integer(i4b), intent(in) :: iEns ! ensemble member + integer(i4b), intent(in) :: segIndex ! reach to process integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) - integer(i4b), intent(in) :: LAKEFLAG ! >0 if processing lakes type(RCHTOPO),intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter - integer(i4b), intent(in), optional :: RSTEP ! retrospective time step offset + type(RCHPRP), intent(inout), allocatable :: RPARAM_in(:) ! River reach parameter type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains integer(i4b), intent(out) :: ierr ! error code @@ -237,24 +134,24 @@ SUBROUTINE qroute_rch(IENS,JRCH, & ! input: array indices character(len=strLen) :: fmt1,fmt2 ! format string character(len=strLen) :: CMESSAGE ! error message for downwind routine - ierr=0; message='qroute_rch/' + ierr=0; message='kwt_rch/' - if(JRCH==ixDesire) then + if(segIndex==ixDesire) then write(iulog,'(2a)') new_line('a'),'** Check kinematic wave tracking routing **' - write(iulog,"(a,x,I10,x,I10)") ' Reach index & ID =', JRCH, NETOPO_in(JRCH)%REACHID + write(iulog,"(a,x,I10,x,I10)") ' Reach index & ID =', segIndex, NETOPO_in(segIndex)%REACHID write(iulog,"(a,x,F20.7,1x,F20.7)") ' time step(T0,T1) =', T0, T1 - write(iulog,'(a,x,F15.7)') ' RPARAM_in%R_SLOPE =', RPARAM_in(JRCH)%R_SLOPE - write(iulog,'(a,x,F15.7)') ' RPARAM_in%R_MAN_N =', RPARAM_in(JRCH)%R_MAN_N - write(iulog,'(a,x,F15.7)') ' RPARAM_in%R_WIDTH =', RPARAM_in(JRCH)%R_WIDTH + write(iulog,'(a,x,F15.7)') ' RPARAM_in%R_SLOPE =', RPARAM_in(segIndex)%R_SLOPE + write(iulog,'(a,x,F15.7)') ' RPARAM_in%R_MAN_N =', RPARAM_in(segIndex)%R_MAN_N + write(iulog,'(a,x,F15.7)') ' RPARAM_in%R_WIDTH =', RPARAM_in(segIndex)%R_WIDTH end if ! ---------------------------------------------------------------------------------------- - ! (1) EXTRACT FLOW FROM UPSTREAM REACHES & APPEND TO THE NON-ROUTED FLOW PARTICLES IN JRCH + ! (1) EXTRACT FLOW FROM UPSTREAM REACHES & APPEND TO THE NON-ROUTED FLOW PARTICLES IN segIndex ! ---------------------------------------------------------------------------------------- - NUPS = count(NETOPO_in(JRCH)%goodBas) ! number of desired upstream reaches - !NUPS = size(NETOPO_in(JRCH)%UREACHI) ! number of upstream reaches + NUPS = count(NETOPO_in(segIndex)%goodBas) ! number of desired upstream reaches + !NUPS = size(NETOPO_in(segIndex)%UREACHI) ! number of upstream reaches if (NUPS.GT.0) then - call getusq_rch(IENS,JRCH,LAKEFLAG,T0,T1,ixDesire, & ! input + call getusq_rch(IENS,segIndex,LAKEFLAG,T0,T1,ixDesire, & ! input NETOPO_in,RPARAM_in,RCHFLX_out, & ! input RCHSTA_out, & ! inout Q_JRCH,TENTRY,T_EXIT,ierr,cmessage,& ! output @@ -265,29 +162,29 @@ SUBROUTINE qroute_rch(IENS,JRCH, & ! input: array indices ierr=20; message=trim(message)//'negative flow extracted from upstream reach'; return endif - if(JRCH==ixDesire)then + if(segIndex==ixDesire)then write(fmt1,'(A,I5,A)') '(A, 1X',size(Q_JRCH),'(1X,G15.4))' write(iulog,'(a)') ' * Wave discharge from upstream reaches (Q_JRCH) [m2/s]:' write(iulog,fmt1) ' Q_JRCH=',(Q_JRCH(IWV), IWV=0,size(Q_JRCH)-1) endif else ! set flow in headwater reaches to modelled streamflow from time delay histogram - RCHFLX_out(IENS,JRCH)%ROUTE(idxKWT)%REACH_Q = RCHFLX_out(IENS,JRCH)%BASIN_QR(1) - if (allocated(RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE)) then - deallocate(RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE,STAT=IERR) + RCHFLX_out(IENS,segIndex)%ROUTE(idxKWT)%REACH_Q = RCHFLX_out(IENS,segIndex)%BASIN_QR(1) + if (allocated(RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE)) then + deallocate(RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE,STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem deallocating space for RCHSTA_out'; return; endif endif - allocate(RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0:0),STAT=ierr) - if(ierr/=0)then; message=trim(message)//'problem allocating space for RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(1)'; return; endif - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0)%QF=-9999 - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0)%TI=-9999 - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0)%TR=-9999 - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0)%RF=.False. - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0)%QM=-9999 - - if(JRCH==ixDesire) then - write(iulog,'(a)') ' * Final discharge (RCHFLX_out(IENS,JRCH)%REACH_Q) [m3/s]:' - write(iulog,'(x,G15.4)') RCHFLX_out(IENS,JRCH)%ROUTE(idxKWT)%REACH_Q + allocate(RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0:0),STAT=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating space for RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(1)'; return; endif + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0)%QF=-9999 + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0)%TI=-9999 + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0)%TR=-9999 + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0)%RF=.False. + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0)%QM=-9999 + + if(segIndex==ixDesire) then + write(iulog,'(a)') ' * Final discharge (RCHFLX_out(IENS,segIndex)%REACH_Q) [m3/s]:' + write(iulog,'(x,G15.4)') RCHFLX_out(IENS,segIndex)%ROUTE(idxKWT)%REACH_Q end if return ! no upstream reaches (routing for sub-basins done using time-delay histogram) endif @@ -305,41 +202,37 @@ SUBROUTINE qroute_rch(IENS,JRCH, & ! input: array indices ! (x) Water use - take out (Qtake is negative) ! ---------------------------------------------------------------------------------------- ! set the retrospective offset - if (.not.present(RSTEP)) then - ROFFSET = 0 - else - ROFFSET = RSTEP - endif + ROFFSET = RSTEP ! set time boundaries T_START = T0 - (T1 - T0)*ROFFSET T_END = T1 - (T1 - T0)*ROFFSET - if (is_flux_wm .and. RCHFLX_out(iens,jrch)%REACH_WM_FLUX /= realMissing) then - call extract_from_rch(iens, jrch, & ! input: ensemble and reach indices - T_START, T_END, & ! input: time [sec] of current time step bounds - RPARAM_in, & ! input: river reach parameters - RCHFLX_out(iens,jrch)%REACH_WM_FLUX, & ! input: target Qtake (minus) - ixDesire, & ! input: - Q_JRCH, T_EXIT, TENTRY, & ! inout: discharge and exit time for particle + if (is_flux_wm .and. RCHFLX_out(iens,segIndex)%REACH_WM_FLUX /= realMissing) then + call extract_from_rch(iens, segIndex, & ! input: ensemble and reach indices + T_START, T_END, & ! input: time [sec] of current time step bounds + RPARAM_in, & ! input: river reach parameters + RCHFLX_out(iens,segIndex)%REACH_WM_FLUX, & ! input: target Qtake (minus) + ixDesire, & ! input: + Q_JRCH, T_EXIT, TENTRY, & ! inout: discharge and exit time for particle ierr,cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif end if ! ---------------------------------------------------------------------------------------- - ! (3) ROUTE FLOW WITHIN THE CURRENT [JRCH] RIVER SEGMENT + ! (3) ROUTE FLOW WITHIN THE CURRENT [segIndex] RIVER SEGMENT ! ---------------------------------------------------------------------------------------- allocate(FROUTE(0:NQ1),STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem allocating space for FROUTE'; return; endif FROUTE(0) = .true.; FROUTE(1:NQ1)=.false. ! init. routing flags - ! route flow through the current [JRCH] river segment (Q_JRCH in units of m2/s) - call kinwav_rch(JRCH,T_START,T_END,ixDesire, & ! input: location and time + ! route flow through the current [segIndex] river segment (Q_JRCH in units of m2/s) + call kinwav_rch(segIndex,T_START,T_END,ixDesire, & ! input: location and time NETOPO_in, RPARAM_in, & ! input: river data structure Q_JRCH(1:NQ1),TENTRY(1:NQ1),T_EXIT(1:NQ1),FROUTE(1:NQ1), & ! inout: kwt states NQ2,ierr,cmessage) ! output: if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - if(JRCH == ixDesire)then + if(segIndex == ixDesire)then write(fmt1,'(A,I5,A)') '(A,1X',NQ1+1,'(1X,G15.4))' write(fmt2,'(A,I5,A)') '(A,1X',NQ1+1,'(1X,L))' write(iulog,'(a)') ' * After routed: wave discharge (Q_JRCH) [m2/s], isExit(FROUTE), entry time (TENTRY) [s], and exit time (T_EXIT) [s]:' @@ -361,13 +254,13 @@ SUBROUTINE qroute_rch(IENS,JRCH, & ! input: array indices if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! m2/s --> m3/s + instantaneous runoff from basin - RCHFLX_out(IENS,JRCH)%ROUTE(idxKWT)%REACH_Q = QNEW(1)*RPARAM_in(JRCH)%R_WIDTH + RCHFLX_out(IENS,JRCH)%BASIN_QR(1) + RCHFLX_out(IENS,segIndex)%ROUTE(idxKWT)%REACH_Q = QNEW(1)*RPARAM_in(segIndex)%R_WIDTH + RCHFLX_out(IENS,segIndex)%BASIN_QR(1) - if(JRCH == ixDesire)then + if(segIndex == ixDesire)then write(iulog,'(a)') ' * Time-ave. wave discharge that exit (QNEW(1)) [m2/s], local-area discharge (RCHFLX_out%BASIN_QR(1)) [m3/s] and Final discharge (RCHFLX_out%REACH_Q) [m3/s]:' write(iulog,"(A,1x,G15.4)") ' QNEW(1) =', QNEW(1) - write(iulog,"(A,1x,G15.4)") ' RCHFLX_out%BASIN_QR(1) =', RCHFLX_out(IENS,JRCH)%BASIN_QR(1) - write(iulog,"(A,1x,G15.4)") ' RCHFLX_out%REACH_Q =', RCHFLX_out(IENS,JRCH)%ROUTE(idxKWT)%REACH_Q + write(iulog,"(A,1x,G15.4)") ' RCHFLX_out%BASIN_QR(1) =', RCHFLX_out(IENS,segIndex)%BASIN_QR(1) + write(iulog,"(A,1x,G15.4)") ' RCHFLX_out%REACH_Q =', RCHFLX_out(IENS,segIndex)%ROUTE(idxKWT)%REACH_Q endif ! ---------------------------------------------------------------------------------------- @@ -381,24 +274,24 @@ SUBROUTINE qroute_rch(IENS,JRCH, & ! input: array indices TIMEI = TENTRY(NR) + & ! (dT/dT) (dT) ( (TENTRY(NR+1)-TENTRY(NR)) / (T_EXIT(NR+1)-T_EXIT(NR)) ) * (T_END-T_EXIT(NR)) ! allocate space for the routed data (+1 to allocate space for the interpolated point) - if (.not.allocated(RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE)) then + if (.not.allocated(RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE)) then ierr=20; message=trim(message)//'RCHSTA_out is not associated'; return else - deallocate(RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE, STAT=ierr) - if(ierr/=0)then; message=trim(message)//'problem deallocating space for RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE'; return; endif - allocate(RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0:NQ2+1),STAT=ierr) ! NQ2 is number of points for kinematic routing - if(ierr/=0)then; message=trim(message)//'problem allocating space for RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0:NQ2+1)'; return; endif + deallocate(RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE, STAT=ierr) + if(ierr/=0)then; message=trim(message)//'problem deallocating space for RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE'; return; endif + allocate(RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0:NQ2+1),STAT=ierr) ! NQ2 is number of points for kinematic routing + if(ierr/=0)then; message=trim(message)//'problem allocating space for RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0:NQ2+1)'; return; endif endif ! insert the interpolated point (TI is irrelevant, as the point is "routed") - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(NR+1)%QF=Q_END; RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(NR+1)%TI=TIMEI - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(NR+1)%TR=T_END; RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(NR+1)%RF=.true. + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(NR+1)%QF=Q_END; RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(NR+1)%TI=TIMEI + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(NR+1)%TR=T_END; RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(NR+1)%RF=.true. ! add the output from kinwave... - skip NR+1 - ! (when JRCH becomes IR routed points will be stripped out & the structures updated again) - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0:NR)%QF=Q_JRCH(0:NR); RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(NR+2:NQ2+1)%QF=Q_JRCH(NR+1:NQ2) - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0:NR)%TI=TENTRY(0:NR); RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(NR+2:NQ2+1)%TI=TENTRY(NR+1:NQ2) - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0:NR)%TR=T_EXIT(0:NR); RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(NR+2:NQ2+1)%TR=T_EXIT(NR+1:NQ2) - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0:NR)%RF=FROUTE(0:NR); RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(NR+2:NQ2+1)%RF=FROUTE(NR+1:NQ2) - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0:NQ2+1)%QM=-9999 + ! (when segIndex becomes IR routed points will be stripped out & the structures updated again) + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0:NR)%QF=Q_JRCH(0:NR); RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(NR+2:NQ2+1)%QF=Q_JRCH(NR+1:NQ2) + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0:NR)%TI=TENTRY(0:NR); RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(NR+2:NQ2+1)%TI=TENTRY(NR+1:NQ2) + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0:NR)%TR=T_EXIT(0:NR); RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(NR+2:NQ2+1)%TR=T_EXIT(NR+1:NQ2) + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0:NR)%RF=FROUTE(0:NR); RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(NR+2:NQ2+1)%RF=FROUTE(NR+1:NQ2) + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0:NQ2+1)%QM=-9999 ! implement water use !if (NUSER.GT.0 .and. UCFFLAG.GE.1) then @@ -412,8 +305,8 @@ SUBROUTINE qroute_rch(IENS,JRCH, & ! input: array indices ! *** ! remove flow particles from the most downstream reach ! if the last reach or lake inlet (and lakes are enabled), remove routed elements from memory - if ((NETOPO_in(JRCH)%DREACHK<=0 ).OR. & ! if the last reach (down reach ID:DREACHK is negative), then there is no downstream reach - (LAKEFLAG.EQ.1.AND.NETOPO_in(JRCH)%LAKINLT)) then ! if lake inlet + if ((NETOPO_in(segIndex)%DREACHK<=0 ).OR. & ! if the last reach (down reach ID:DREACHK is negative), then there is no downstream reach + (LAKEFLAG.EQ.1.AND.NETOPO_in(segIndex)%LAKINLT)) then ! if lake inlet ! copy data to a temporary wave if (allocated(NEW_WAVE)) then deallocate(NEW_WAVE,STAT=IERR) @@ -421,19 +314,19 @@ SUBROUTINE qroute_rch(IENS,JRCH, & ! input: array indices endif allocate(NEW_WAVE(0:NN),STAT=IERR) ! NN = number non-routed (the zero element is the last routed point) if(ierr/=0)then; message=trim(message)//'problem allocating space for NEW_WAVE'; return; endif - NEW_WAVE(0:NN) = RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(NR+1:NQ2+1) ! +1 because of the interpolated point + NEW_WAVE(0:NN) = RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(NR+1:NQ2+1) ! +1 because of the interpolated point ! re-size wave structure - if (allocated(RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE)) then - deallocate(RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE,STAT=IERR) + if (allocated(RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE)) then + deallocate(RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE,STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem deallocating space for RCHSTA_out'; return; endif endif - allocate(RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0:NN),STAT=IERR) ! again, the zero element for the last routed point + allocate(RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0:NN),STAT=IERR) ! again, the zero element for the last routed point if(ierr/=0)then; message=trim(message)//'problem allocating space for RCHSTA_out'; return; endif ! copy data back to the wave structure and deallocate space for the temporary wave - RCHSTA_out(IENS,JRCH)%LKW_ROUTE%KWAVE(0:NN) = NEW_WAVE(0:NN) - endif ! (if JRCH is the last reach) + RCHSTA_out(IENS,segIndex)%LKW_ROUTE%KWAVE(0:NN) = NEW_WAVE(0:NN) + endif ! (if segIndex is the last reach) - END SUBROUTINE qroute_rch + END SUBROUTINE kwt_rch ! ********************************************************************* ! subroutine: wave discharge mod to extract/infect water from the JRCH reach diff --git a/route/build/src/mc_route.f90 b/route/build/src/mc_route.f90 index 5023bf48..251c5c76 100644 --- a/route/build/src/mc_route.f90 +++ b/route/build/src/mc_route.f90 @@ -3,135 +3,36 @@ MODULE mc_route_module ! muskingum-cunge routing USE nrtype -! data types -USE dataTypes, ONLY: STRFLX ! fluxes in each reach -USE dataTypes, ONLY: STRSTA ! state in each reach -USE dataTypes, ONLY: RCHTOPO ! Network topology -USE dataTypes, ONLY: RCHPRP ! Reach parameter -USE dataTypes, ONLY: mcRCH ! MC specific state data structure -USE dataTypes, ONLY: subbasin_omp ! mainstem+tributary data strucuture -! global data -USE public_var, ONLY: iulog ! i/o logical unit number -USE public_var, ONLY: realMissing ! missing value for real number -USE public_var, ONLY: integerMissing ! missing value for integer number -USE public_var, ONLY: qmodOption ! qmod option (use 1==direct insertion) -USE globalData, ONLY: idxMC ! index of IRF method -! subroutines: general +USE dataTypes, ONLY: STRFLX ! fluxes in each reach +USE dataTypes, ONLY: STRSTA ! state in each reach +USE dataTypes, ONLY: RCHTOPO ! Network topology +USE dataTypes, ONLY: RCHPRP ! Reach parameter +USE dataTypes, ONLY: mcRCH ! MC specific state data structure +USE public_var, ONLY: iulog ! i/o logical unit number +USE public_var, ONLY: realMissing ! missing value for real number +USE public_var, ONLY: integerMissing ! missing value for integer number +USE public_var, ONLY: qmodOption ! qmod option (use 1==direct insertion) +USE globalData, ONLY: idxMC ! index of IRF method USE water_balance, ONLY: comp_reach_wb ! compute water balance error -USE perf_mod, ONLY: t_startf,t_stopf ! timing start/stop -USE model_utils, ONLY: handle_err +USE base_route, ONLY: base_route_rch implicit none private -public::mc_route +public::mc_route_rch -CONTAINS - - ! ********************************************************************* - ! subroutine: perform muskingum-cunge routing through the river network - ! ********************************************************************* - SUBROUTINE mc_route(iens, & ! input: ensemble index - river_basin, & ! input: river basin information (mainstem, tributary outlet etc.) - T0,T1, & ! input: start and end of the time step - ixDesire, & ! input: reachID to be checked by on-screen pringing - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,message, & ! output: error control - ixSubRch) ! optional input: subset of reach indices to be processed - - USE public_var, ONLY: is_lake_sim ! logical whether or not lake should be simulated - - implicit none - ! Argument variables - integer(i4b), intent(in) :: iEns ! ensemble member - type(subbasin_omp), intent(in), allocatable :: river_basin(:) ! river basin information (mainstem, tributary outlet etc.) - real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) - integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output - type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter - type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data - type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message ! error message - integer(i4b), intent(in), optional :: ixSubRch(:) ! subset of reach indices to be processed - ! local variables - character(len=strLen) :: cmessage ! error message for downwind routine - logical(lgt), allocatable :: doRoute(:) ! logical to indicate which reaches are processed - integer(i4b) :: nOrder ! number of stream order - integer(i4b) :: nTrib ! number of tributary basins - integer(i4b) :: nSeg ! number of reaches in the network - integer(i4b) :: iSeg, jSeg ! loop indices - reach - integer(i4b) :: iTrib ! loop indices - branch - integer(i4b) :: ix ! loop indices stream order - - ierr=0; message='mc_route/' - - ! number of reach check - if (size(NETOPO_in)/=size(RCHFLX_out(iens,:))) then - ierr=20; message=trim(message)//'sizes of NETOPO and RCHFLX mismatch'; return - endif - - nSeg = size(NETOPO_in) - - allocate(doRoute(nSeg), stat=ierr) +type, extends(base_route_rch) :: mc_route_rch + CONTAINS + procedure, pass :: route => mc_rch +end type mc_route_rch - if (present(ixSubRch))then - doRoute(:)=.false. - doRoute(ixSubRch) = .true. ! only subset of reaches are on - else - doRoute(:)=.true. ! every reach is on - endif - - nOrder = size(river_basin) - - call t_startf('route/mc') - - do ix = 1, nOrder - - nTrib=size(river_basin(ix)%branch) - -!$OMP PARALLEL DO schedule(dynamic,1) & ! chunk size of 1 -!$OMP private(jSeg, iSeg) & ! private for a given thread -!$OMP private(ierr, cmessage) & ! private for a given thread -!$OMP shared(T0,T1) & ! 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(RPARAM_in) & ! data structure shared -!$OMP shared(RCHSTA_out) & ! data structure shared -!$OMP shared(RCHFLX_out) & ! data structure shared -!$OMP shared(ix, iEns, ixDesire) & ! indices shared -!$OMP firstprivate(nTrib) - trib:do iTrib = 1,nTrib - seg:do iSeg=1,river_basin(ix)%branch(iTrib)%nRch - jSeg = river_basin(ix)%branch(iTrib)%segIndex(iSeg) - if (.not. doRoute(jSeg)) cycle - call mc_rch(iEns,jSeg, & ! input: array indices - ixDesire, & ! input: index of the desired reach - T0,T1, & ! input: start and end of the time step - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,cmessage) ! output: error control - if(ierr/=0) call handle_err(ierr, trim(message)//trim(cmessage)) - end do seg - end do trib -!$OMP END PARALLEL DO - - end do - - call t_stopf('route/mc') - - END SUBROUTINE mc_route +CONTAINS ! ********************************************************************* ! subroutine: perform muskingum-cunge routing for one segment ! ********************************************************************* - SUBROUTINE mc_rch(iEns, segIndex, & ! input: index of runoff ensemble to be processed + SUBROUTINE mc_rch(this, & + iEns, segIndex, & ! input: index of runoff ensemble to be processed ixDesire, & ! input: reachID to be checked by on-screen pringing T0,T1, & ! input: start and end of the time step NETOPO_in, & ! input: reach topology data structure @@ -142,12 +43,13 @@ SUBROUTINE mc_rch(iEns, segIndex, & ! input: index of runoff ensemble to be proc implicit none ! Argument variables + class(mc_route_rch) :: this integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed integer(i4b), intent(in) :: segIndex ! segment where routing is performed integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter + type(RCHPRP), intent(inout), allocatable :: RPARAM_in(:) ! River reach parameter type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains integer(i4b), intent(out) :: ierr ! error code From 901438119dc6888ed6d2431ae9a7549bbc552257 Mon Sep 17 00:00:00 2001 From: Naoki Mizukami Date: Mon, 26 Jun 2023 06:15:39 -0600 Subject: [PATCH 4/7] routine to instantiate active routing method objects --- route/build/src/init_model_data.f90 | 81 ++++++++++++++++++++++++----- 1 file changed, 69 insertions(+), 12 deletions(-) diff --git a/route/build/src/init_model_data.f90 b/route/build/src/init_model_data.f90 index cfc006b1..e2f6d14e 100644 --- a/route/build/src/init_model_data.f90 +++ b/route/build/src/init_model_data.f90 @@ -90,18 +90,15 @@ SUBROUTINE init_model(cfile_name, ierr, message) ! used to read control files and namelist and broadcast to all processors - ! shared data used - USE public_var, ONLY: ancil_dir - USE public_var, ONLY: input_dir - USE public_var, ONLY: param_nml - USE public_var, ONLY: gageMetaFile - USE public_var, ONLY: outputAtGage - ! subroutines: populate metadata - USE popMetadat_module, ONLY: popMetadat ! populate metadata - ! subroutines: model control - USE read_control_module, ONLY: read_control ! read the control file - USE read_param_module, ONLY: read_param ! read the routing parameters - USE process_gage_meta, ONLY: read_gage_meta ! process gauge metadata + USE public_var, ONLY: ancil_dir + USE public_var, ONLY: input_dir + USE public_var, ONLY: param_nml + USE public_var, ONLY: gageMetaFile + USE public_var, ONLY: outputAtGage + USE popMetadat_module, ONLY: popMetadat ! populate metadata + USE read_control_module, ONLY: read_control ! read the control file + USE read_param_module, ONLY: read_param ! read the routing parameters + USE process_gage_meta, ONLY: read_gage_meta ! process gauge metadata implicit none ! Argument variables @@ -129,6 +126,9 @@ SUBROUTINE init_model(cfile_name, ierr, message) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif end if + call init_route_method(ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + END SUBROUTINE init_model @@ -608,4 +608,61 @@ SUBROUTINE init_ntopo(nHRU_out, nRch_out, END SUBROUTINE init_ntopo + ! ********************************************************************* + ! public subroutine: initialize routing method object + ! ********************************************************************* + SUBROUTINE init_route_method(ierr, message) + ! + ! DESCRIPTION: + ! Instantiate a collection of routing method objects + ! + USE globalData, ONLY: rch_routes ! routing methods instantiated + USE globalData, ONLY: routeMethods ! Active routing method + USE public_var, ONLY: accumRunoff ! routing method ID + USE public_var, ONLY: impulseResponseFunc ! routing method ID + USE public_var, ONLY: kinematicWaveTracking ! routing method ID + USE public_var, ONLY: kinematicWave ! routing method ID + USE public_var, ONLY: muskingumCunge ! routing method ID + USE public_var, ONLY: diffusiveWave ! routing method ID + USE accum_route_module, ONLY: accum_route_rch ! routing routine: accumulation instantaneous runoff + USE irf_route_module, ONLY: irf_route_rch ! routing routine: Impulse response function + USE kwt_route_module, ONLY: kwt_route_rch ! routing routine: Lagrangian kinematic + USE kw_route_module, ONLY: kwe_route_rch ! routing routine: kinematic + USE mc_route_module, ONLY: mc_route_rch ! routing routine: muskingum + USE dfw_route_module, ONLY: dfw_route_rch ! routing routine: diffusive + + implicit none + ! Argument variables: + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + ! Local variables: + character(len=strLen) :: cmessage ! error message from subroutine + integer(i4b) :: ix + + ierr=0; message='init_route_method/' + + allocate(rch_routes(size(routeMethods)), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [rch_routes]'; return; endif + + do ix=1, size(routeMethods) + select case (routeMethods(ix)) + case (accumRunoff) + allocate(accum_route_rch :: rch_routes(ix)%rch_route) + case (impulseResponseFunc) + allocate(irf_route_rch :: rch_routes(ix)%rch_route) + case (kinematicWaveTracking) + allocate(kwt_route_rch :: rch_routes(ix)%rch_route) + case (kinematicWave) + allocate(kwe_route_rch :: rch_routes(ix)%rch_route) + case (muskingumCunge) + allocate(mc_route_rch :: rch_routes(ix)%rch_route) + case (diffusiveWave) + allocate(dfw_route_rch :: rch_routes(ix)%rch_route) + case default + ierr=20; message=trim(message)//'no valid routing method'; return + end select + end do + + END SUBROUTINE init_route_method + END MODULE init_model_data From bf00cfe89f23c91acf0e473651a6e04664d7334b Mon Sep 17 00:00:00 2001 From: Naoki Mizukami Date: Mon, 26 Jun 2023 06:16:34 -0600 Subject: [PATCH 5/7] use new routing objects to route through the network --- route/build/src/main_route.f90 | 244 ++++++++++++++++++--------------- 1 file changed, 136 insertions(+), 108 deletions(-) diff --git a/route/build/src/main_route.f90 b/route/build/src/main_route.f90 index 992d8e12..c9a43bdf 100644 --- a/route/build/src/main_route.f90 +++ b/route/build/src/main_route.f90 @@ -1,26 +1,16 @@ MODULE main_route_module -USE nrtype ! variable types, etc. -! data structures -USE dataTypes, ONLY: STRSTA ! state in each reach -USE dataTypes, ONLY: STRFLX ! fluxes in each reach -USE dataTypes, ONLY: RCHTOPO ! Network topology -USE dataTypes, ONLY: RCHPRP ! Reach parameter -USE dataTypes, ONLY: runoff ! runoff data type -USE dataTypes, ONLY: subbasin_omp ! mainstem+tributary data structures -USE globalData, ONLY: routeMethods -USE globalData, ONLY: nRoutes -! mapping HRU runoff to reach -USE process_remap_module, ONLY: basin2reach -! subroutines: basin routing -USE basinUH_module, ONLY: IRF_route_basin ! perform UH convolution for basin routing -! subroutines: river routing -USE accum_runoff_module, ONLY: accum_runoff ! upstream flow accumulation -USE irf_route_module, ONLY: irf_route ! unit hydrograph (impulse response function) routing method -USE kwt_route_module, ONLY: kwt_route ! lagrangian kinematic wave routing method -USE kw_route_module, ONLY: kw_route ! kinematic wave routing method -USE mc_route_module, ONLY: mc_route ! muskingum-cunge routing method -USE dfw_route_module, ONLY: dfw_route ! diffusive wave routing method +USE nrtype ! variable types, etc. +USE dataTypes, ONLY: STRSTA ! state in each reach +USE dataTypes, ONLY: STRFLX ! fluxes in each reach +USE dataTypes, ONLY: RCHTOPO ! Network topology +USE dataTypes, ONLY: RCHPRP ! Reach parameter +USE dataTypes, ONLY: runoff ! runoff data type +USE dataTypes, ONLY: subbasin_omp ! mainstem+tributary data structures +USE globalData, ONLY: routeMethods ! Active routing method IDs +USE public_var, ONLY: is_lake_sim ! lake simulation flag +USE process_remap_module, ONLY: basin2reach ! remap HRU variable to reach +USE basinUH_module, ONLY: IRF_route_basin ! perform UH convolution for basin routing implicit none @@ -53,16 +43,9 @@ SUBROUTINE main_route(iens, & ! input: ensemble index ! 2. Process a list of reach indices (in terms of NETOPO_in etc.) given by ixRchProcessed ! 3. basinRunoff_in is given in the order of NETOPO_in(:)%HRUIX. - USE public_var, ONLY: doesBasinRoute - USE public_var, ONLY: accumRunoff - USE public_var, ONLY: impulseResponseFunc - USE public_var, ONLY: kinematicWaveTracking - USE public_var, ONLY: kinematicWave - USE public_var, ONLY: muskingumCunge - USE public_var, ONLY: diffusiveWave - USE globalData, ONLY: onRoute ! logical to indicate which routing method(s) is on USE globalData, ONLY: TSEC ! beginning/ending of simulation time step [sec] - USE public_var, ONLY: is_lake_sim ! logical whether or not lake should be simulated + USE globalData, ONLY: rch_routes ! + USE public_var, ONLY: doesBasinRoute USE public_var, ONLY: is_flux_wm ! logical whether or not fluxes should be passed USE public_var, ONLY: is_vol_wm ! logical whether or not target volume should be passed @@ -85,18 +68,15 @@ SUBROUTINE main_route(iens, & ! input: ensemble index character(len=strLen), intent(out) :: message ! error message ! local variables character(len=strLen) :: cmessage ! error message of downwind routine - real(dp) :: T0,T1 ! beginning/ending of simulation time step [sec] real(dp), allocatable :: reachRunoff_local(:) ! reach runoff (m/s) real(dp), allocatable :: reachEvapo_local(:) ! reach evaporation (m/s) real(dp), allocatable :: reachPrecip_local(:) ! reach precipitation (m/s) integer(i4b) :: nSeg ! number of reach to be processed integer(i4b) :: iSeg ! index of reach + integer(i4b) :: ix ! loop index ierr=0; message = "main_routing/" - ! define the start and end of the time step - T0=TSEC(1); T1=TSEC(2) - ! number of reaches to be processed nSeg = size(ixRchProcessed) @@ -191,86 +171,134 @@ SUBROUTINE main_route(iens, & ! input: ensemble index endif ! 3. subroutine: river reach routing - ! perform upstream flow accumulation - if (onRoute(accumRunoff)) then - call accum_runoff(iens, & ! input: ensemble index - river_basin, & ! input: river basin data type - ixDesire, & ! input: index of verbose reach - NETOPO_in, & ! input: reach topology data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr, cmessage, & ! output: error controls - ixRchProcessed) ! optional input: indices of reach to be routed + do ix=1,size(routeMethods) + call route_network(rch_routes(ix)%rch_route, & + iens, & ! input: ensemble index + river_basin, & ! input: river basin data type + TSEC(1), TSEC(2), & ! input: start and end of the time step since simulation start [sec] + ixDesire, & ! input: index of verbose reach + NETOPO_in, & ! input: reach topology data structure + RPARAM_in, & ! input: reach parameter + RCHSTA_out, & ! inout: reach state data structure + RCHFLX_out, & ! inout: reach flux data structure + ierr, cmessage, & ! output: error controls + ixRchProcessed) ! optional input: indices of reach to be routed if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - endif + end do - if (onRoute(kinematicWaveTracking)) then - call kwt_route(iens, & ! input: ensemble index - river_basin, & ! input: river basin data type - T0,T1, & ! input: start and end of the time step - ixDesire, & ! input: index of verbose reach - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,cmessage, & ! output: error control - ixRchProcessed) ! optional input: indices of reach to be routed - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - end if + END SUBROUTINE main_route - if (onRoute(impulseResponseFunc)) then - call irf_route(iens, & ! input: ensemble index - river_basin, & ! input: river basin data type - ixDesire, & ! input: index of verbose reach - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,cmessage, & ! output: error control - ixRchProcessed) ! optional input: indices of reach to be routed - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - end if + ! ********************************************************************* + ! private subroutine: route throughout river network + ! ********************************************************************* + SUBROUTINE route_network(rch_route, & ! input: + iens, & ! input: ensemble index + river_basin, & ! input: river basin information (mainstem, tributary outlet etc.) + T0,T1, & ! input: start and end of the time step since simulation start [sec] + ixDesire, & ! input: reachID to be checked by on-screen pringing + NETOPO_in, & ! input: reach topology data structure + RPARAM_in, & ! input: reach parameter data structure + RCHSTA_out, & ! inout: reach state data structure + RCHFLX_out, & ! inout: reach flux data structure + ierr,message, & ! output: error control + ixSubRch) ! optional input: subset of reach indices to be processed - if (onRoute(kinematicWave)) then - call kw_route(iens, & ! input: ensemble index - river_basin, & ! input: river basin data type - T0,T1, & ! input: start and end of the time step - ixDesire, & ! input: index of verbose reach - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,cmessage, & ! output: error control - ixRchProcessed) ! optional input: indices of reach to be routed - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - end if + USE perf_mod, ONLY: t_startf,t_stopf ! timing start/stop + USE lake_route_module, ONLY: lake_route ! lake route module + USE base_route, ONLY: base_route_rch ! + USE model_utils, ONLY: handle_err - if (onRoute(muskingumCunge)) then - call mc_route(iens, & ! input: ensemble index - river_basin, & ! input: river basin data type - T0,T1, & ! input: start and end of the time step - ixDesire, & ! input: index of verbose reach - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,cmessage, & ! output: error control - ixRchProcessed) ! optional input: indices of reach to be routed - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - end if + implicit none + ! Argument variables + class(base_route_rch), intent(in), allocatable :: rch_route + integer(i4b), intent(in) :: iEns ! ensemble member + type(subbasin_omp), intent(in), allocatable :: river_basin(:) ! river basin information (mainstem, tributary outlet etc.) + real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) + integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output + type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology + type(RCHPRP), intent(inout), allocatable :: RPARAM_in(:) ! River reach parameter + type(STRSTA), intent(inout) :: RCHSTA_out(:,:) ! reach state data + type(STRFLX), intent(inout) :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + integer(i4b), intent(in), optional :: ixSubRch(:) ! subset of reach indices to be processed + ! local variables + character(len=strLen) :: cmessage ! error message for downwind routine + logical(lgt), allocatable :: doRoute(:) ! logical to indicate which reaches are processed + integer(i4b) :: nOrder ! number of stream order + integer(i4b) :: nTrib ! number of tributary basins + integer(i4b) :: nSeg ! number of reaches in the network + integer(i4b) :: iSeg, jSeg ! loop indices - reach + integer(i4b) :: iTrib ! loop indices - branch + integer(i4b) :: ix ! loop indices stream order - if (onRoute(diffusiveWave)) then - call dfw_route(iens, & ! input: ensemble index - river_basin, & ! input: river basin data type - T0,T1, & ! input: start and end of the time step - ixDesire, & ! input: index of verbose reach - NETOPO_in, & ! input: reach topology data structure - RPARAM_in, & ! input: reach parameter data structure - RCHSTA_out, & ! inout: reach state data structure - RCHFLX_out, & ! inout: reach flux data structure - ierr,cmessage, & ! output: error control - ixRchProcessed) ! optional input: indices of reach to be routed - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - end if + ierr=0; message='route_network/' - END SUBROUTINE main_route + nSeg = size(RCHFLX_out(iens,:)) + + ! number of reach check + if (size(NETOPO_in)/=nSeg) then + ierr=20; message=trim(message)//'sizes of NETOPO and RCHFLX mismatch'; return + endif + + allocate(doRoute(nSeg), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating space for [doRoute]'; return; endif + + if (present(ixSubRch))then + doRoute(:)=.false. + doRoute(ixSubRch) = .true. ! only subset of reaches are on + else + doRoute(:)=.true. ! every reach is on + endif + + nOrder = size(river_basin) + + call t_startf('route_network') ! timing start + + ! routing through river network + do ix = 1, nOrder + nTrib=size(river_basin(ix)%branch) +!$OMP PARALLEL DO schedule(dynamic,1) & ! chunk size of 1 +!$OMP private(jSeg, iSeg) & ! private for a given thread +!$OMP private(ierr, cmessage) & ! private for a given thread +!$OMP shared(T0,T1) & ! 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(RPARAM_in) & ! data structure shared +!$OMP shared(RCHSTA_out) & ! 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 + if ((NETOPO_in(jseg)%islake).and.(is_lake_sim)) then + call lake_route(iEns, jSeg, & ! input: ensemble and reach indices + ixDesire, & ! input: index of verbose reach + NETOPO_in, & ! input: reach topology data structure + RPARAM_in, & ! input: reach parameter data structure + RCHFLX_out, & ! inout: reach flux data structure + ierr,cmessage) ! output: error control + else + call rch_route%route(iEns,jSeg, & ! input: array indices + ixDesire, & ! input: index of verbose reach + T0,T1, & ! input: start and end of the time step + NETOPO_in, & ! input: reach topology data structure + RPARAM_in, & ! input: reach parameter data structure + RCHSTA_out, & ! inout: reach state data structure + RCHFLX_out, & ! inout: reach flux data structure + ierr,cmessage) ! output: error control + end if + if(ierr/=0) call handle_err(ierr, trim(message)//trim(cmessage)) + end do ! reach index + end do ! tributary +!$OMP END PARALLEL DO + end do ! basin loop + + call t_stopf('route_network') + + END SUBROUTINE route_network END MODULE main_route_module From 48252d40f723bf2597ff0ebf12c4b8ca1fddcfff Mon Sep 17 00:00:00 2001 From: Naoki Mizukami Date: Sat, 1 Jul 2023 15:13:58 -0600 Subject: [PATCH 6/7] added comments as suggested by Erik --- route/build/src/accum_runoff.f90 | 14 +++++++------- route/build/src/base_route.f90 | 20 ++++++++++++++------ route/build/src/dfw_route.f90 | 4 ++-- route/build/src/init_model_data.f90 | 4 ++-- route/build/src/irf_route.f90 | 6 +++--- route/build/src/kwe_route.f90 | 4 ++-- route/build/src/kwt_route.f90 | 4 ++-- route/build/src/mc_route.f90 | 4 ++-- 8 files changed, 34 insertions(+), 26 deletions(-) diff --git a/route/build/src/accum_runoff.f90 b/route/build/src/accum_runoff.f90 index ac3d53d2..fe3a2e90 100644 --- a/route/build/src/accum_runoff.f90 +++ b/route/build/src/accum_runoff.f90 @@ -1,4 +1,4 @@ -MODULE accum_route_module +MODULE accum_runoff_module ! Accumulate upstream flow instantaneously ! This is not used as routed runoff. @@ -16,19 +16,19 @@ MODULE accum_route_module implicit none private -public::accum_route_rch +public::accum_runoff_rch -type, extends(base_route_rch) :: accum_route_rch +type, extends(base_route_rch) :: accum_runoff_rch CONTAINS procedure, pass :: route => accum_inst_runoff -end type accum_route_rch +end type accum_runoff_rch CONTAINS ! ********************************************************************* ! subroutine: perform accumulate immediate upstream flow ! ********************************************************************* - SUBROUTINE accum_inst_runoff(this, & + SUBROUTINE accum_inst_runoff(this, & ! "accum_runoff_rchr" object to bound this procedure iEns, & ! input: index of runoff ensemble to be processed segIndex, & ! input: index of reach to be processed ixDesire, & ! input: reachID to be checked by on-screen pringing @@ -40,7 +40,7 @@ SUBROUTINE accum_inst_runoff(this, & ierr, message) ! output: error control implicit none ! Argument variables - class(accum_route_rch) :: this + class(accum_runoff_rch) :: this ! "accum_runoff_rchr" object to bound this procedure integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed integer(i4b), intent(in) :: segIndex ! segment where routing is performed integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output @@ -95,4 +95,4 @@ SUBROUTINE accum_inst_runoff(this, & END SUBROUTINE accum_inst_runoff -END MODULE accum_route_module +END MODULE accum_runoff_module diff --git a/route/build/src/base_route.f90 b/route/build/src/base_route.f90 index 6086625a..6c3f40d2 100644 --- a/route/build/src/base_route.f90 +++ b/route/build/src/base_route.f90 @@ -1,20 +1,23 @@ MODULE base_route - ! Abstract class for routing method + ! Description: Definition of base (or template) reach routing method class. + ! this abstract class needs to be extended to specific routing method types for instantiation. implicit none private - public:: base_route_rch - public:: routeContainer + public:: base_route_rch ! base (abstract) reach routing method class (to be extended to specific) + public:: routeContainer ! a holder of instantiated reach routing method object ! --- routing method container + ! This container (holder) include instantiated reach routing method + ! type :: routeContainer class(base_route_rch), allocatable :: rch_route end type - ! --- routing method container + ! --- base (abstract or template) reach routing method type, abstract :: base_route_rch CONTAINS procedure(sub_route_rch), deferred :: route @@ -22,7 +25,7 @@ MODULE base_route ABSTRACT INTERFACE - SUBROUTINE sub_route_rch(this, & ! + SUBROUTINE sub_route_rch(this, & ! object to bound the procedure iEns, & ! input: ensemble index segIndex, & ! input: reach indice ixDesire, & ! input: index of verbose reach @@ -33,6 +36,11 @@ SUBROUTINE sub_route_rch(this, & ! RCHFLX_out, & ! inout: reach flux data structure ierr, message) ! output: error control + ! Description: perform a routing at a given reach (segIndex) and time step + ! reach parameters (RPARAM), river network topology (NETOPO) to get upstream location, + ! state (RCHSTA) and flux (RCHFLX) are required for a set of input + ! ixDesire is index of reach where more information is writting in log along the computation + USE nrtype USE dataTypes, ONLY: STRFLX ! fluxes in each reach USE dataTypes, ONLY: STRSTA ! states in each reach @@ -41,7 +49,7 @@ SUBROUTINE sub_route_rch(this, & ! import base_route_rch ! Arguments - class(base_route_rch) :: this + class(base_route_rch) :: this ! object to bound the procedure integer(i4b), intent(in) :: iEns ! ensemble member integer(i4b), intent(in) :: segIndex ! ensemble member integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output diff --git a/route/build/src/dfw_route.f90 b/route/build/src/dfw_route.f90 index 91e388d4..2531ac8b 100644 --- a/route/build/src/dfw_route.f90 +++ b/route/build/src/dfw_route.f90 @@ -12,9 +12,9 @@ MODULE dfw_route_module USE public_var, ONLY: realMissing ! missing value for real number USE public_var, ONLY: integerMissing ! missing value for integer number USE public_var, ONLY: qmodOption ! qmod option (use 1==direct insertion) -USE globalData, ONLY: idxDW +USE globalData, ONLY: idxDW ! routing method index for diffusive wave USE water_balance, ONLY: comp_reach_wb ! compute water balance error -USE base_route, ONLY: base_route_rch +USE base_route, ONLY: base_route_rch ! base (abstract) reach routing method class implicit none diff --git a/route/build/src/init_model_data.f90 b/route/build/src/init_model_data.f90 index e2f6d14e..311c3f9f 100644 --- a/route/build/src/init_model_data.f90 +++ b/route/build/src/init_model_data.f90 @@ -624,7 +624,7 @@ SUBROUTINE init_route_method(ierr, message) USE public_var, ONLY: kinematicWave ! routing method ID USE public_var, ONLY: muskingumCunge ! routing method ID USE public_var, ONLY: diffusiveWave ! routing method ID - USE accum_route_module, ONLY: accum_route_rch ! routing routine: accumulation instantaneous runoff + USE accum_runoff_module,ONLY: accum_runoff_rch ! routing routine: accumulation instantaneous runoff USE irf_route_module, ONLY: irf_route_rch ! routing routine: Impulse response function USE kwt_route_module, ONLY: kwt_route_rch ! routing routine: Lagrangian kinematic USE kw_route_module, ONLY: kwe_route_rch ! routing routine: kinematic @@ -647,7 +647,7 @@ SUBROUTINE init_route_method(ierr, message) do ix=1, size(routeMethods) select case (routeMethods(ix)) case (accumRunoff) - allocate(accum_route_rch :: rch_routes(ix)%rch_route) + allocate(accum_runoff_rch :: rch_routes(ix)%rch_route) case (impulseResponseFunc) allocate(irf_route_rch :: rch_routes(ix)%rch_route) case (kinematicWaveTracking) diff --git a/route/build/src/irf_route.f90 b/route/build/src/irf_route.f90 index 7cca3e3d..32a73564 100644 --- a/route/build/src/irf_route.f90 +++ b/route/build/src/irf_route.f90 @@ -9,10 +9,10 @@ MODULE irf_route_module USE public_var, ONLY: iulog ! i/o logical unit number USE public_var, ONLY: realMissing ! missing value for real number USE public_var, ONLY: integerMissing ! missing value for integer number -USE public_var, ONLY: dt -USE globalData, ONLY: idxIRF ! index of IRF method +USE public_var, ONLY: dt ! simulation time step [sec] +USE globalData, ONLY: idxIRF ! routing method index for IRF method USE water_balance, ONLY: comp_reach_wb ! compute water balance error -USE base_route, ONLY: base_route_rch +USE base_route, ONLY: base_route_rch ! base (abstract) reach routing method class implicit none diff --git a/route/build/src/kwe_route.f90 b/route/build/src/kwe_route.f90 index 6117d95b..23d6e981 100644 --- a/route/build/src/kwe_route.f90 +++ b/route/build/src/kwe_route.f90 @@ -12,9 +12,9 @@ MODULE kw_route_module USE public_var, ONLY: realMissing ! missing value for real number USE public_var, ONLY: integerMissing ! missing value for integer number USE public_var, ONLY: qmodOption ! qmod option (use 1==direct insertion) -USE globalData, ONLY: idxKW ! loop index of kw routing +USE globalData, ONLY: idxKW ! routing method index for kinematic wwave USE water_balance, ONLY: comp_reach_wb ! compute water balance error -USE base_route, ONLY: base_route_rch +USE base_route, ONLY: base_route_rch ! base (abstract) reach routing method class implicit none diff --git a/route/build/src/kwt_route.f90 b/route/build/src/kwt_route.f90 index 3ee0e0e3..c9eb4c15 100644 --- a/route/build/src/kwt_route.f90 +++ b/route/build/src/kwt_route.f90 @@ -11,9 +11,9 @@ 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: idxKWT ! index of KWT method +USE globalData, ONLY: idxKWT ! routing method index for lagrangian kinematic wave method USE nr_utils, ONLY: arth ! Num. Recipies utilities -USE base_route, ONLY: base_route_rch +USE base_route, ONLY: base_route_rch ! base (abstract) reach routing method class implicit none diff --git a/route/build/src/mc_route.f90 b/route/build/src/mc_route.f90 index 251c5c76..e4730c37 100644 --- a/route/build/src/mc_route.f90 +++ b/route/build/src/mc_route.f90 @@ -12,9 +12,9 @@ MODULE mc_route_module USE public_var, ONLY: realMissing ! missing value for real number USE public_var, ONLY: integerMissing ! missing value for integer number USE public_var, ONLY: qmodOption ! qmod option (use 1==direct insertion) -USE globalData, ONLY: idxMC ! index of IRF method +USE globalData, ONLY: idxMC ! routing method index for muskingum method USE water_balance, ONLY: comp_reach_wb ! compute water balance error -USE base_route, ONLY: base_route_rch +USE base_route, ONLY: base_route_rch ! base (abstract) reach routing method class implicit none From 8c37e9c1aaf0d7ccdeaddb5d142216b8ff0b7963 Mon Sep 17 00:00:00 2001 From: Naoki Mizukami Date: Sat, 1 Jul 2023 15:40:23 -0600 Subject: [PATCH 7/7] more comments --- route/build/src/accum_runoff.f90 | 4 ++-- route/build/src/base_route.f90 | 7 ++++--- route/build/src/dfw_route.f90 | 4 ++-- route/build/src/irf_route.f90 | 4 ++-- route/build/src/kwe_route.f90 | 4 ++-- route/build/src/kwt_route.f90 | 4 ++-- route/build/src/mc_route.f90 | 4 ++-- 7 files changed, 16 insertions(+), 15 deletions(-) diff --git a/route/build/src/accum_runoff.f90 b/route/build/src/accum_runoff.f90 index fe3a2e90..2310f976 100644 --- a/route/build/src/accum_runoff.f90 +++ b/route/build/src/accum_runoff.f90 @@ -10,8 +10,8 @@ MODULE accum_runoff_module USE dataTypes, ONLY: RCHTOPO ! Network topology USE dataTypes, ONLY: RCHPRP ! Reach parameter USE public_var, ONLY: iulog ! i/o logical unit number -USE globalData, ONLY: idxSUM ! index of accumulation method -USE base_route, ONLY: base_route_rch +USE globalData, ONLY: idxSUM ! routing method index for runoff accumulation method +USE base_route, ONLY: base_route_rch ! base (abstract) reach routing method class implicit none diff --git a/route/build/src/base_route.f90 b/route/build/src/base_route.f90 index 6c3f40d2..aab4092e 100644 --- a/route/build/src/base_route.f90 +++ b/route/build/src/base_route.f90 @@ -1,7 +1,8 @@ MODULE base_route ! Description: Definition of base (or template) reach routing method class. - ! this abstract class needs to be extended to specific routing method types for instantiation. + ! this abstract class needs to be extended to specific routing method types for + ! implementation and instantiation. implicit none @@ -11,7 +12,6 @@ MODULE base_route ! --- routing method container ! This container (holder) include instantiated reach routing method - ! type :: routeContainer class(base_route_rch), allocatable :: rch_route end type @@ -36,7 +36,8 @@ SUBROUTINE sub_route_rch(this, & ! object to bound the procedure RCHFLX_out, & ! inout: reach flux data structure ierr, message) ! output: error control - ! Description: perform a routing at a given reach (segIndex) and time step + ! Description: template interfade for reach routing subroutine + ! to perform a routing (after instantiated) at a given reasch (segIndex) and time step ! reach parameters (RPARAM), river network topology (NETOPO) to get upstream location, ! state (RCHSTA) and flux (RCHFLX) are required for a set of input ! ixDesire is index of reach where more information is writting in log along the computation diff --git a/route/build/src/dfw_route.f90 b/route/build/src/dfw_route.f90 index 2531ac8b..a7134257 100644 --- a/route/build/src/dfw_route.f90 +++ b/route/build/src/dfw_route.f90 @@ -31,7 +31,7 @@ MODULE dfw_route_module ! ********************************************************************* ! subroutine: perform diffusive wave routing for one segment ! ********************************************************************* - SUBROUTINE dfw_rch(this, & + SUBROUTINE dfw_rch(this, & ! dfw_route_rch object to bound this procedure iEns, segIndex, & ! input: index of runoff ensemble to be processed ixDesire, & ! input: reachID to be checked by on-screen pringing T0,T1, & ! input: start and end of the time step @@ -43,7 +43,7 @@ SUBROUTINE dfw_rch(this, & implicit none ! Argument variables - class(dfw_route_rch) :: this + class(dfw_route_rch) :: this ! dfw_route_rch object to bound this procedure integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed integer(i4b), intent(in) :: segIndex ! segment where routing is performed integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output diff --git a/route/build/src/irf_route.f90 b/route/build/src/irf_route.f90 index 32a73564..d649bfd3 100644 --- a/route/build/src/irf_route.f90 +++ b/route/build/src/irf_route.f90 @@ -29,7 +29,7 @@ MODULE irf_route_module ! ********************************************************************* ! subroutine: perform one segment route UH routing ! ********************************************************************* - SUBROUTINE irf_rch(this, & ! + SUBROUTINE irf_rch(this, & ! irf_route_rch object to bound this procedure iEns, & ! input: index of runoff ensemble segIndex, & ! input: reach index ixDesire, & ! input: reachID to be checked by on-screen pringing @@ -44,7 +44,7 @@ SUBROUTINE irf_rch(this, & ! implicit none ! Argument variables - class(irf_route_rch) :: this + class(irf_route_rch) :: this ! irf_route_rch object to bound this procedure integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed integer(i4b), intent(in) :: segIndex ! segment where routing is performed integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output diff --git a/route/build/src/kwe_route.f90 b/route/build/src/kwe_route.f90 index 23d6e981..3c2d6b02 100644 --- a/route/build/src/kwe_route.f90 +++ b/route/build/src/kwe_route.f90 @@ -33,7 +33,7 @@ MODULE kw_route_module ! ********************************************************************* ! subroutine: perform one segment route KW routing ! ********************************************************************* - SUBROUTINE kw_rch(this, & + SUBROUTINE kw_rch(this, & ! kwe_route_rch object to bound this procedure iEns, segIndex, & ! input: index of runoff ensemble to be processed ixDesire, & ! input: reachID to be checked by on-screen pringing T0,T1, & ! input: start and end of the time step @@ -44,7 +44,7 @@ SUBROUTINE kw_rch(this, & ierr, message) ! output: error control implicit none ! Argument variables - class(kwe_route_rch) :: this + class(kwe_route_rch) :: this ! kwe_route_rch object to bound this procedure integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed integer(i4b), intent(in) :: segIndex ! segment where routing is performed integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output diff --git a/route/build/src/kwt_route.f90 b/route/build/src/kwt_route.f90 index c9eb4c15..e9edd77b 100644 --- a/route/build/src/kwt_route.f90 +++ b/route/build/src/kwt_route.f90 @@ -33,7 +33,7 @@ MODULE kwt_route_module ! ********************************************************************* ! subroutine: route kinematic waves at one segment ! ********************************************************************* - SUBROUTINE kwt_rch(this, & + SUBROUTINE kwt_rch(this, & ! kwt_route_rch object to bound this procedure iEns, & ! input: ensemble index segIndex, & ! input: index of reach to be processed ixDesire, & ! input: index of the reach for verbose output @@ -96,7 +96,7 @@ SUBROUTINE kwt_rch(this, & ! ---------------------------------------------------------------------------------------- implicit none ! Argument variables - class(kwt_route_rch) :: this + class(kwt_route_rch) :: this ! kwt_route_rch object to bound this procedure integer(i4b), intent(in) :: iEns ! ensemble member integer(i4b), intent(in) :: segIndex ! reach to process integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output diff --git a/route/build/src/mc_route.f90 b/route/build/src/mc_route.f90 index e4730c37..42ff3d19 100644 --- a/route/build/src/mc_route.f90 +++ b/route/build/src/mc_route.f90 @@ -31,7 +31,7 @@ MODULE mc_route_module ! ********************************************************************* ! subroutine: perform muskingum-cunge routing for one segment ! ********************************************************************* - SUBROUTINE mc_rch(this, & + SUBROUTINE mc_rch(this, & ! mc_route_rch object to bound this procedure iEns, segIndex, & ! input: index of runoff ensemble to be processed ixDesire, & ! input: reachID to be checked by on-screen pringing T0,T1, & ! input: start and end of the time step @@ -43,7 +43,7 @@ SUBROUTINE mc_rch(this, & implicit none ! Argument variables - class(mc_route_rch) :: this + class(mc_route_rch) :: this ! mc_route_rch object to bound this procedure integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed integer(i4b), intent(in) :: segIndex ! segment where routing is performed integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output