diff --git a/route/build/src/dataTypes.f90 b/route/build/src/dataTypes.f90 index d074a0b3..a303e26f 100644 --- a/route/build/src/dataTypes.f90 +++ b/route/build/src/dataTypes.f90 @@ -170,6 +170,8 @@ MODULE dataTypes integer(i4b) , allocatable :: hru_id(:) ! id of HM_HRUs or RN_HRUs at which runoff is stored (size: nSpace(1)) integer(i4b) , allocatable :: hru_ix(:) ! Index of RN_HRUs associated with river network (used only if HM_HRUs = RN_HRUs) real(dp) , allocatable :: basinRunoff(:)! remapped river network catchment runoff (size: number of nHRU) + real(dp) , allocatable :: basinEvapo(:) ! remapped river network catchment runoff (size: number of nHRU) + real(dp) , allocatable :: basinPrecip(:)! remapped river network catchment runoff (size: number of nHRU) end type runoff ! ---------- reach parameters ---------------------------------------------------------------------------- diff --git a/route/build/src/globalData.f90 b/route/build/src/globalData.f90 index 345b68a9..e377f9fa 100644 --- a/route/build/src/globalData.f90 +++ b/route/build/src/globalData.f90 @@ -180,6 +180,10 @@ module globalData type(runoff) , public :: runoff_data ! HRU runoff data structure for one time step for LSM HRUs and River network HRUs real(dp) , allocatable , public :: basinRunoff_trib(:) ! HRU runoff array (m/s) for tributaries real(dp) , allocatable , public :: basinRunoff_main(:) ! HRU runoff array (m/s) for mainstem + real(dp) , allocatable , public :: basinEvapo_trib(:) ! HRU evaporation array (m/s) for tributaries + real(dp) , allocatable , public :: basinEvapo_main(:) ! HRU evaporation array (m/s) for mainstem + real(dp) , allocatable , public :: basinPrecip_trib(:) ! HRU precipitation array (m/s) for tributaries + real(dp) , allocatable , public :: basinPrecip_main(:) ! HRU precipitation array (m/s) for mainstem ! domain data ! MPI diff --git a/route/build/src/main_route.f90 b/route/build/src/main_route.f90 index e6933c4a..29ca0fce 100644 --- a/route/build/src/main_route.f90 +++ b/route/build/src/main_route.f90 @@ -38,6 +38,8 @@ SUBROUTINE main_route(& ! input iens, & ! ensemble index basinRunoff_in, & ! basin (i.e.,HRU) runoff (m/s) + basinEvapo_in, & ! basin (i.e.,HRU) evaporation (m/s) + basinPrecip_in, & ! basin (i.e.,HRU) precipitation (m/s) ixRchProcessed, & ! indices of reach to be routed river_basin, & ! OMP basin decomposition NETOPO_in, & ! reach topology data structure @@ -64,12 +66,15 @@ SUBROUTINE main_route(& USE public_var, ONLY: kinematicWaveEuler USE public_var, ONLY: impulseResponseFunc 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 implicit none ! input integer(i4b), intent(in) :: iens ! ensemble member real(dp), allocatable, intent(in) :: basinRunoff_in(:) ! basin (i.e.,HRU) runoff (m/s) + real(dp), allocatable, intent(in) :: basinEvapo_in(:) ! basin (i.e.,HRU) evaporation (m/s) + real(dp), allocatable, intent(in) :: basinPrecip_in(:) ! basin (i.e.,HRU) precipitation (m/s) integer(i4b), allocatable, intent(in) :: ixRchProcessed(:) ! indices of reach to be routed type(subbasin_omp), allocatable, intent(in) :: river_basin(:) ! OMP basin decomposition type(RCHTOPO), allocatable, intent(in) :: NETOPO_in(:) ! River Network topology @@ -85,12 +90,15 @@ SUBROUTINE main_route(& 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 ! initialize errors ierr=0; message = "main_routing/" + ! define the start and end of the time step T0=TSEC(0); T1=TSEC(1) diff --git a/route/build/src/mpi_process.f90 b/route/build/src/mpi_process.f90 index d47458b4..43487597 100644 --- a/route/build/src/mpi_process.f90 +++ b/route/build/src/mpi_process.f90 @@ -705,6 +705,10 @@ subroutine mpi_route(pid, & ! input: proc id USE globalData, ONLY: nHRU_mainstem ! number of mainstem HRUs USE globalData, ONLY: basinRunoff_main ! mainstem only HRU runoff USE globalData, ONLY: basinRunoff_trib ! tributary only HRU runoff + USE globalData, ONLY: basinEvapo_main ! mainstem only HRU Evaporation + USE globalData, ONLY: basinEvapo_trib ! tributary only HRU Evaporation + USE globalData, ONLY: basinPrecip_main ! mainstem only HRU Precipitation + USE globalData, ONLY: basinPrecip_trib ! tributary only HRU Precipitation USE globalData, ONLY: river_basin_trib ! tributary OMP domain data structure USE globalData, ONLY: river_basin_main ! mainstem OMP domain data structure USE globalData, ONLY: nRch_mainstem ! number of mainstem reaches @@ -798,6 +802,8 @@ subroutine mpi_route(pid, & ! input: proc id ! Perform routing call main_route(iens, & ! input: ensemble index basinRunoff_trib, & ! input: basin (i.e.,HRU) runoff (m/s) + basinEvapo_trib, & ! input: basin (i.e. HRU) Evapo (m/s) + basinPrecip_trib, & ! input: basin (i.e. HRU) Precip (m/s) ixRchProcessed, & ! input: indices of reach to be routed river_basin_trib, & ! input: OMP basin decomposition NETOPO_trib, & ! input: reach topology data structure @@ -874,6 +880,8 @@ subroutine mpi_route(pid, & ! input: proc id call main_route(iens, & ! input: ensemble index basinRunoff_main, & ! input: basin (i.e.,HRU) runoff (m/s) + basinEvapo_main, & ! input: basin (i.e. HRU) Evapo (m/s) + basinPrecip_main, & ! input: basin (i.e. HRU) Precip (m/s) ixRchProcessed, & ! input: indices of reach to be routed river_basin_main, & ! input: OMP basin decomposition NETOPO_main, & ! input: reach topology data structure @@ -933,6 +941,12 @@ SUBROUTINE scatter_runoff(nNodes, comm, & ! mpi variables: number nodes, communi USE globalData, ONLY: hru_per_proc ! number of hrus assigned to each proc (i.e., node) USE globalData, ONLY: basinRunoff_main ! HRU runoff holder for mainstem USE globalData, ONLY: basinRunoff_trib ! HRU runoff holder for tributary + USE globalData, ONLY: basinEvapo_main ! HRU evaporation holder for mainstem + USE globalData, ONLY: basinEvapo_trib ! HRU evaporation holder for tributary + USE globalData, ONLY: basinPrecip_main ! HRU precipitation holder for mainstem + USE globalData, ONLY: basinPrecip_trib ! HRU precipitation holder for tributary + USE public_var, ONLY: is_lake_sim ! logical whether or not lake should be simulated + ! input variables integer(i4b), intent(in) :: nNodes ! number of processes (MPI) @@ -941,7 +955,10 @@ SUBROUTINE scatter_runoff(nNodes, comm, & ! mpi variables: number nodes, communi integer(i4b), intent(out) :: ierr ! error code character(len=strLen), intent(out) :: message ! error message ! local variables + integer(i4b) :: iHru,jHru ! loop indices real(dp) :: basinRunoff_local(nHRU) ! temporal basin runoff (m/s) for whole domain + real(dp) :: basinEvapo_local(nHRU) ! temporal basin runoff (m/s) for whole domain + real(dp) :: basinPrecip_local(nHRU) ! temporal basin runoff (m/s) for whole domain character(len=strLen) :: cmessage ! error message from a subroutine ierr=0; message='scatter_runoff/' @@ -955,9 +972,27 @@ SUBROUTINE scatter_runoff(nNodes, comm, & ! mpi variables: number nodes, communi end if basinRunoff_main(:) = runoff_data%basinRunoff(:) + if (is_lake_sim) then + + ! if only single proc is used, all evaporation is stored in mainstem evaporation array + if (.not. allocated(basinEvapo_main)) then + allocate(basinEvapo_main(nHRU), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating array for [basinEvapo_main]'; return; endif + end if + basinEvapo_main(:) = runoff_data%basinEvapo(:) + + ! if only single proc is used, all precipitation is stored in mainstem precipitation array + if (.not. allocated(basinPrecip_main)) then + allocate(basinPrecip_main(nHRU), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating array for [basinPrecip_main]'; return; endif + end if + basinPrecip_main(:) = runoff_data%basinPrecip(:) + + end if + else - ! sort the basin runoff in terms of nodes/domains + ! sort the basin runoff, precipitation and evaporation in terms of nodes/domains if (masterproc) then ! this is a root process if (.not. allocated(basinRunoff_main)) then @@ -965,11 +1000,35 @@ SUBROUTINE scatter_runoff(nNodes, comm, & ! mpi variables: number nodes, communi if(ierr/=0)then; message=trim(message)//'problem allocating array for [basinRunoff_main]'; return; endif endif - ! runoff at hru in mainstem - basinRunoff_local(1:nHRU) = runoff_data%basinRunoff(1:nHRU) + if (is_lake_sim) then + + if (.not. allocated(basinEvapo_main)) then + allocate(basinEvapo_main(nHRU_mainstem), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating array for [basinEvapo_main]'; return; endif + + endif + + if (.not. allocated(basinPrecip_main)) then + allocate(basinPrecip_main(nHRU_mainstem), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating array for [basinPrecip_main]'; return; endif + endif + + end if + + ! runoff at hru in mainstem and tributaries + basinRunoff_local(1:nHRU) = runoff_data%basinRunoff(1:nHRU) basinRunoff_main(1:nHRU_mainstem) = basinRunoff_local(1:nHRU_mainstem) + + ! evaporation and precipitation at main channel and tributaries + if (is_lake_sim) then + basinEvapo_local (1:nHRU) = runoff_data%basinEvapo(1:nHRU) + basinPrecip_local(1:nHRU) = runoff_data%basinPrecip(1:nHRU) + basinEvapo_main (1:nHRU_mainstem) = basinEvapo_local (1:nHRU_mainstem) + basinPrecip_main(1:nHRU_mainstem) = basinPrecip_local(1:nHRU_mainstem) + end if + end if call shr_mpi_barrier(comm, message) @@ -981,6 +1040,20 @@ SUBROUTINE scatter_runoff(nNodes, comm, & ! mpi variables: number nodes, communi ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + if (is_lake_sim) then + call shr_mpi_scatterV(basinEvapo_local(nHRU_mainstem+1:nHRU), & + hru_per_proc(0:nNodes-1), & + basinEvapo_trib, & + ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + call shr_mpi_scatterV(basinPrecip_local(nHRU_mainstem+1:nHRU), & + hru_per_proc(0:nNodes-1), & + basinPrecip_trib, & + ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end if + end if END SUBROUTINE scatter_runoff diff --git a/route/build/src/public_var.f90 b/route/build/src/public_var.f90 index 5d126342..9d75a388 100644 --- a/route/build/src/public_var.f90 +++ b/route/build/src/public_var.f90 @@ -94,6 +94,8 @@ module public_var integer(i4b) ,public :: routOpt = integerMissing ! routing scheme options 0-> both, 1->IRF, 2->KWT, otherwise error integer(i4b) ,public :: doesBasinRoute = 1 ! basin routing options 0-> no, 1->IRF, otherwise error integer(i4b) ,public :: doesAccumRunoff = 1 ! option to delayed runoff accumulation over all the upstream reaches + logical(lgt),public :: is_lake_sim = .false. ! logical whether or not lakes are activated in simulation + logical(lgt),public :: is_wm_sim = .false. ! logical whether or not water management componenets, abstraction, injections and target volums are provided and will be used in the simulation ! RIVER NETWORK TOPOLOGY character(len=strLen),public :: fname_ntopOld = '' ! old filename containing stream network topology information logical(lgt) ,public :: ntopAugmentMode = .false. ! option for river network augmentation mode. terminate the program after writing augmented ntopo. @@ -104,6 +106,8 @@ module public_var ! RUNOFF FILE character(len=strLen),public :: fname_qsim = '' ! simulated runoff netCDF name character(len=strLen),public :: vname_qsim = '' ! variable name for simulated runoff + character(len=strLen),public :: vname_evapo = '' ! variable name for actual evapoartion + character(len=strLen),public :: vname_precip = '' ! variable name for precipitation character(len=strLen),public :: vname_time = '' ! variable name for time character(len=strLen),public :: vname_hruid = '' ! variable name for runoff hru id character(len=strLen),public :: dname_time = '' ! dimension name for time diff --git a/route/build/src/read_control.f90 b/route/build/src/read_control.f90 index c2b9f3b7..9e47803d 100644 --- a/route/build/src/read_control.f90 +++ b/route/build/src/read_control.f90 @@ -117,6 +117,8 @@ subroutine read_control(ctl_fname, err, message) case(''); read(cData,*,iostat=io_error) doesBasinRoute ! basin routing options 0-> no, 1->IRF, otherwise error case(''); read(cData,*,iostat=io_error) doesAccumRunoff ! option to delayed runoff accumulation over all the upstream reaches. 0->no, 1->yes case('' ); read(cData,*,iostat=io_error) idSegOut ! desired outlet reach id (if -9999 --> route over the entire network) + case(''); read(cData,*,iostat=io_error) is_lake_sim ! logical whether or not lakes are simulated + case(''); read(cData,*,iostat=io_error) is_wm_sim ! logical whether or not water balance components, abstraction, injection and target volumes should be used in simulation ! RIVER NETWORK TOPOLOGY case(''); fname_ntopOld = trim(cData) ! name of file containing stream network topology information case(''); read(cData,*,iostat=io_error) ntopAugmentMode ! option for river network augmentation mode. terminate the program after writing augmented ntopo. @@ -126,6 +128,8 @@ subroutine read_control(ctl_fname, err, message) ! RUNOFF FILE case(''); fname_qsim = trim(cData) ! name of text file containing nc file names and their order for runoff case(''); vname_qsim = trim(cData) ! name of runoff variable + case(''); vname_evapo = trim(cData) ! name of actual evapoartion variable + case(''); vname_precip = trim(cData) ! name of precipitation variable case(''); vname_time = trim(cData) ! name of time variable in the runoff file case(''); vname_hruid = trim(cData) ! name of the HRU id case(''); dname_time = trim(cData) ! name of time variable in the runoff file diff --git a/route/build/src/standalone/get_basin_runoff.f90 b/route/build/src/standalone/get_basin_runoff.f90 index 81c5ba6b..a164bb5e 100644 --- a/route/build/src/standalone/get_basin_runoff.f90 +++ b/route/build/src/standalone/get_basin_runoff.f90 @@ -21,7 +21,11 @@ subroutine get_hru_runoff(ierr, message) ! output: error control USE public_var, only:input_dir ! directory containing input data USE public_var, only:fname_qsim ! simulated runoff netCDF name USE public_var, only:vname_qsim ! varibale runoff in netCDF file + USE public_var, only:vname_evapo ! varibale actual evaporation in netCDF file + USE public_var, only:vname_precip ! varibale precipitation in netCDF file USE public_var, only:is_remap ! logical whether or not runnoff needs to be mapped to river network HRU + USE public_var, only:is_lake_sim ! logical whether or not lake should be simulated + USE public_var, only:is_wm_sim ! logical whether or not water management components should be read, abstraction, injection and target volume USE globalData, only:iTime_local ! iTime index for the given netcdf file USE globalData, only:nHRU ! number of routing sub-basin USE globalData, only:runoff_data ! data structure to hru runoff data @@ -67,6 +71,60 @@ subroutine get_hru_runoff(ierr, message) ! output: error control if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif end if + if (is_lake_sim) then ! if is_lake_sim if true then read actual evaporation and preciptation + + ! get the actual evaporation - runoff_data%sim(:) or %sim2D(:,:) + call read_runoff_data(trim(input_dir)//trim(fname_qsim), & ! input: filename + trim(vname_evapo), & ! input: varname + iTime_local, & ! input: time index + runoff_data, & ! inout: runoff data structure + ierr, cmessage) ! output: error control + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + ! initialize runoff_data%basinEvapo + if ( allocated(runoff_data%basinEvapo) ) then + deallocate(runoff_data%basinEvapo, stat=ierr) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end if + allocate(runoff_data%basinEvapo(nHRU), stat=ierr) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + ! Get river network HRU runoff into runoff_data data structure + if (is_remap) then ! remap LSM simulated flux to the HRUs in the river network + call remap_runoff(runoff_data, remap_data, runoff_data%basinEvapo, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + else ! runoff is already remapped to river network HRUs + call sort_runoff(runoff_data, runoff_data%basinEvapo, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end if + + + ! get the precepitation - runoff_data%sim(:) or %sim2D(:,:) + call read_runoff_data(trim(input_dir)//trim(fname_qsim), & ! input: filename + trim(vname_precip), & ! input: varname + iTime_local, & ! input: time index + runoff_data, & ! inout: runoff data structure + ierr, cmessage) ! output: error control + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + ! initialize runoff_data%basinPrecip + if ( allocated(runoff_data%basinPrecip) ) then + deallocate(runoff_data%basinPrecip, stat=ierr) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end if + allocate(runoff_data%basinPrecip(nHRU), stat=ierr) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + ! Get river network HRU runoff into runoff_data data structure + if (is_remap) then ! remap LSM simulated flux to the HRUs in the river network + call remap_runoff(runoff_data, remap_data, runoff_data%basinPrecip, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + else ! runoff is already remapped to river network HRUs + call sort_runoff(runoff_data, runoff_data%basinPrecip, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end if + end if + end subroutine get_hru_runoff -end module get_runoff \ No newline at end of file +end module get_runoff