Skip to content

Commit

Permalink
Add initial support for gage-assisted diversions in channel routing (#…
Browse files Browse the repository at this point in the history
…756)

* Add initial support for gage-assisted diversions (Type 3) in channel routing
---------

Co-authored-by: Soren Rasmussen <s.c.rasmussen@gmail.com>
  • Loading branch information
rcabell and scrasmussen authored Jan 7, 2025
1 parent 5da9e0c commit 2682cee
Show file tree
Hide file tree
Showing 18 changed files with 461 additions and 26 deletions.
16 changes: 8 additions & 8 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,13 @@ endif()


# set project name and version numbers
project(WRF_Hydro LANGUAGES Fortran)
set(WRF_Hydro_VERSION_MAJOR 5)
set(WRF_Hydro_VERSION_MINOR 3)
set(WRF_Hydro_VERSION_PATCH 0)
set(National_Water_Model_VERSION_MAJOR 3)
set(National_Water_Model_VERSION_MINOR 0)
set(National_Water_Model_VERSION_PATCH beta)
project (WRF_Hydro LANGUAGES Fortran C)
set (WRF_Hydro_VERSION_MAJOR 5)
set (WRF_Hydro_VERSION_MINOR 4)
set (WRF_Hydro_VERSION_PATCH 0)
set (National_Water_Model_VERSION_MAJOR 3)
set (National_Water_Model_VERSION_MINOR 1)
set (National_Water_Model_VERSION_PATCH beta)

# set cmake to work with MPI Fortran
find_package(MPI REQUIRED)
Expand Down Expand Up @@ -211,7 +211,7 @@ elseif (CMAKE_Fortran_COMPILER_ID MATCHES "Intel.*")
# set compile flags for ifort
message( "-- Using ifort")
set(CMAKE_Fortran_FLAGS "-fpp -w -ftz -align all -fno-alias -fp-model precise -FR -convert big_endian")
set(CMAKE_Fortran_FLAGS_RELEASE "-O2")
set(CMAKE_Fortran_FLAGS_RELEASE "-O2 -march=core-avx2")
set(CMAKE_Fortran_FLAGS_DEBUG "-O0 -g -traceback")
elseif ((CMAKE_Fortran_COMPILER_ID MATCHES "PGI.*") OR (CMAKE_Fortran_COMPILER_ID MATCHES "NVHPC.*"))
message("-- Using NVHPC / PGI")
Expand Down
5 changes: 4 additions & 1 deletion src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ add_subdirectory("Debug_Utilities")
add_subdirectory("Routing/Overland")
add_subdirectory("Routing/Subsurface")
add_subdirectory("Routing/Reservoirs")
add_subdirectory("Routing/Diversions")
add_subdirectory("Data_Rec")
add_subdirectory("Routing")
add_subdirectory("HYDRO_drv")
Expand Down Expand Up @@ -112,13 +113,15 @@ if (HYDRO_LSM MATCHES "NoahMP")

if (WRF_HYDRO_NUDGING STREQUAL "1")
target_link_libraries(wrfhydro.exe hydro_nudging)
target_link_libraries(wrfhydro.exe hydro_routing_diversions)
add_dependencies(wrfhydro.exe hydro_nudging)
add_dependencies(wrfhydro.exe hydro_routing_diversions)
endif()

# bash commands to copy namelists to the Run directory
set(BASH_CP_HRLDAS_NML "if [[ ! -f ${CMAKE_BINARY_DIR}/Run/namelist.hrldas ]]\; then cp ${PROJECT_SOURCE_DIR}/src/template/NoahMP/namelist.hrldas ${CMAKE_BINARY_DIR}/Run \; fi\;")
set(BASH_CP_HYDRO_NML "if [[ ! -f ${CMAKE_BINARY_DIR}/Run/hydro.namelist ]]\; then cp ${PROJECT_SOURCE_DIR}/src/template/HYDRO/hydro.namelist ${CMAKE_BINARY_DIR}/Run \; fi\;")

add_custom_command(TARGET wrfhydro.exe POST_BUILD
COMMAND mkdir -p ${CMAKE_BINARY_DIR}/Run
COMMAND cp ${PROJECT_SOURCE_DIR}/tests/ctests/run_dir_makefile.mk ${CMAKE_BINARY_DIR}/Run/Makefile
Expand Down
6 changes: 4 additions & 2 deletions src/Data_Rec/module_namelist.F90
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ subroutine read_rt_nlst(nlst)
character(len=256) :: route_lake_f=""
character(len=256) :: route_direction_f=""
character(len=256) :: route_order_f=""
character(len=256) :: diversions_file=""
logical :: reservoir_persistence_usgs
logical :: reservoir_persistence_usace
character(len=256) :: reservoir_parameter_file=""
Expand Down Expand Up @@ -105,8 +106,8 @@ subroutine read_rt_nlst(nlst)
RT_OPTION, CHANRTSWCRT, channel_option, &
SUBRTSWCRT,OVRTSWCRT,AGGFACTRT, dtrt_ter,dtrt_ch,dxrt,&
GwSpinCycles, GwPreCycles, GwSpinUp, GwPreDiag, GwPreDiagInterval, gwIhShift, &
GWBASESWCRT, gwChanCondSw, gwChanCondConstIn, gwChanCondConstOut , &
route_topo_f,route_chan_f,route_link_f, compound_channel, route_lake_f, &
GWBASESWCRT, gwChanCondSw, gwChanCondConstIn, gwChanCondConstOut, &
route_topo_f,route_chan_f,route_link_f, compound_channel, route_lake_f, diversions_file, &
reservoir_persistence_usgs, reservoir_persistence_usace, reservoir_parameter_file, reservoir_usgs_timeslice_path, &
reservoir_usace_timeslice_path, reservoir_observation_lookback_hours, reservoir_observation_update_time_interval_seconds, &
reservoir_rfc_forecasts, reservoir_rfc_forecasts_time_series_path, reservoir_rfc_forecasts_lookback_hours, &
Expand Down Expand Up @@ -248,6 +249,7 @@ subroutine read_rt_nlst(nlst)
nlst%DEEPGWSPIN = DEEPGWSPIN
nlst%SOLVEG_INITSWC = SOLVEG_INITSWC
nlst%reservoir_obs_dir = "testDirectory"
nlst%diversions_file = diversions_file
nlst%reservoir_persistence_usgs = reservoir_persistence_usgs
nlst%reservoir_persistence_usace = reservoir_persistence_usace
nlst%reservoir_parameter_file = reservoir_parameter_file
Expand Down
1 change: 1 addition & 0 deletions src/Data_Rec/module_namelist_inc.F90
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ module module_namelist_inc
character(len=256) :: route_chan_f=""
character(len=256) :: route_link_f=""
character(len=256) :: route_lake_f=""
character(len=256) :: diversions_file=""
logical :: reservoir_persistence_usgs
logical :: reservoir_persistence_usace
character(len=256) :: reservoir_parameter_file=""
Expand Down
6 changes: 6 additions & 0 deletions src/HYDRO_drv/module_HYDRO_drv.F90
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ module module_HYDRO_drv
#endif
use module_hydro_stop, only: HYDRO_stop
use module_UDMAP, only: get_basn_area_nhd
use module_channel_diversions, only: init_diversions
use netcdf

implicit none
Expand Down Expand Up @@ -1583,6 +1584,11 @@ subroutine HYDRO_ini(ntime, did,ix0,jx0, vegtyp,soltyp)
if(nlst(did)%CHANRTSWCRT .ne. 0) call init_stream_nudging
#endif

!#ifdef WRF_HYDRO_DIVERSIONS
! TODO: should this check to make sure we have nudging on too? [RC]
call init_diversions(nlst(did)%diversions_file, nlst(did)%timeSlicePath)
!#endif


! if (trim(nlst_rt(did)%restart_file) == "") then
! output at the initial time
Expand Down
8 changes: 6 additions & 2 deletions src/OrchestratorLayer/config.F90
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,7 @@ module config_base
character(len=256) :: route_link_f=""
character(len=256) :: route_lake_f=""
integer :: lake_option
character(len=256) :: diversions_file=""
logical :: reservoir_persistence_usgs
logical :: reservoir_persistence_usace
character(len=256) :: reservoir_parameter_file=""
Expand Down Expand Up @@ -170,7 +171,7 @@ module config_base
integer :: maxAgePairsBiasPersist
logical :: invDistTimeWeightBias
logical :: noConstInterfBias
character(len=256) :: timeSlicePath
character(len=256) :: timeSlicePath = "./nudgingTimeSliceObs/"
integer :: nLastObs
integer :: bucket_loss
integer :: imperv_adj
Expand Down Expand Up @@ -554,6 +555,7 @@ subroutine init_namelist_rt_field(did)
integer :: channel_loss_option = 0
character(len=256) :: route_lake_f=""
integer :: lake_option !0: lakes off 1: level pool 2: passthrough, 3: reservoir da
character(len=256) :: diversions_file=""
logical :: reservoir_persistence_usgs
logical :: reservoir_persistence_usace
character(len=256) :: reservoir_parameter_file=""
Expand Down Expand Up @@ -628,7 +630,7 @@ subroutine init_namelist_rt_field(did)
SUBRTSWCRT,OVRTSWCRT,AGGFACTRT, dtrt_ter,dtrt_ch,dxrt,&
GwSpinCycles, GwPreCycles, GwSpinUp, GwPreDiag, GwPreDiagInterval, gwIhShift, &
GWBASESWCRT, gwChanCondSw, gwChanCondConstIn, gwChanCondConstOut , &
route_topo_f,route_chan_f,route_link_f, compound_channel, channel_loss_option, lake_option, route_lake_f, &
route_topo_f,route_chan_f,route_link_f, compound_channel, channel_loss_option, lake_option, route_lake_f, diversions_file, &
route_direction_f,route_order_f,gwbasmskfil, &
geo_finegrid_flnm, gwstrmfil,GW_RESTART,RSTRT_SWC,TERADJ_SOLAR, sys_cpl, &
order_to_write , rst_typ, rst_bi_in, rst_bi_out, gwsoilcpl, &
Expand Down Expand Up @@ -876,6 +878,8 @@ subroutine init_namelist_rt_field(did)
nlst(did)%route_link_f = route_link_f
nlst(did)%route_lake_f = route_lake_f

nlst(did)%diversions_file = diversions_file

nlst(did)%reservoir_persistence_usgs = reservoir_persistence_usgs
nlst(did)%reservoir_persistence_usace = reservoir_persistence_usace
nlst(did)%reservoir_parameter_file = reservoir_parameter_file
Expand Down
1 change: 1 addition & 0 deletions src/Routing/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -29,4 +29,5 @@ target_link_libraries(hydro_routing
hydro_routing_reservoirs_hybrid
hydro_data_rec
hydro_routing_reservoirs_rfc
hydro_routing_diversions
)
10 changes: 10 additions & 0 deletions src/Routing/Diversions/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
add_library(hydro_routing_diversions STATIC
module_diversions.F90
module_diversions_timeslice.F90
)

add_dependencies(hydro_routing_diversions hydro_orchestrator)
add_dependencies(hydro_routing_diversions fortglob)

target_link_libraries(hydro_routing_diversions PUBLIC hydro_orchestrator)
target_link_libraries(hydro_routing_diversions PUBLIC fortglob)
177 changes: 177 additions & 0 deletions src/Routing/Diversions/module_diversions.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,177 @@
module module_channel_diversions
use netcdf
use iso_fortran_env, only: int8, int16, int64
use ieee_arithmetic, only: ieee_is_nan

use module_diversions_timeslice, only: get_flow_for_gage, init_timeslices
use module_hydro_stop, only: hydro_stop

implicit none

type diversion_t
character(len=128) :: name

character(len=16) :: da_src, da_dest
integer(kind=int8) :: type_div, type_src, type_dest
integer(kind=int64) :: id_src, id_dest, src_index, dest_index
real :: capacity, fraction
integer(kind=int16) :: lookback

real :: persisted_flow_src, persisted_flow_dest
end type

logical :: diversions_active = .false.
integer :: ndivs = 0

type(diversion_t), allocatable :: diversions(:)

character(*), parameter :: free = '(*(g0,1x))'

contains
subroutine init_diversions(diversions_file, timeslice_path)
character(*), intent(in) :: diversions_file
character(*), intent(in) :: timeslice_path

integer :: g, i, ierr = 0
character(len=20) :: istr
character(len=256) :: char_tmp

integer :: ncid, dimid
integer :: name_vid, type_div_vid, type_src_vid, type_dest_vid, da_src_vid, da_dest_vid
integer :: id_src_vid, id_dest_vid, capacity_vid, fraction_vid, lookback_vid

if (len_trim(diversions_file) > 0) then
print *, "Loading diversions data from " // trim(diversions_file)
ierr = nf90_open(trim(diversions_file), NF90_NOWRITE, ncid)
if (ierr /= 0) call hydro_stop("Could not open diversions file: " // trim(diversions_file))
ierr = nf90_inq_dimid(ncid, "diversion", dimid)
if (ierr /= 0) call hydro_stop("Error reading diversions file: " // trim(diversions_file))
ierr = nf90_inquire_dimension(ncid, dimid, len=ndivs)
if (ierr /= 0) call hydro_stop("Error reading diversions file: " // trim(diversions_file))

write (istr, *) ndivs
print *, "Diversions file has " // trim(adjustl(istr)) // " diversions"

! get fields
ierr = 0
ierr = ierr + nf90_inq_varid(ncid, "Diversion_Name", name_vid)
ierr = ierr + nf90_inq_varid(ncid, "DivType", type_div_vid)
ierr = ierr + nf90_inq_varid(ncid, "FromType", type_src_vid)
ierr = ierr + nf90_inq_varid(ncid, "ToType", type_dest_vid)
ierr = ierr + nf90_inq_varid(ncid, "DA_Src", da_src_vid)
ierr = ierr + nf90_inq_varid(ncid, "DA_Dest", da_dest_vid)
ierr = ierr + nf90_inq_varid(ncid, "DivFrom", id_src_vid)
ierr = ierr + nf90_inq_varid(ncid, "DivTo", id_dest_vid)
ierr = ierr + nf90_inq_varid(ncid, "DivCap", capacity_vid)
ierr = ierr + nf90_inq_varid(ncid, "DivFrac", fraction_vid)
ierr = ierr + nf90_inq_varid(ncid, "Lookback", lookback_vid)

if (ierr /= 0) call hydro_stop("Error occurred accessing diversion file variables")

if (ndivs > 0) then
diversions_active = .true.

! Read the timeslice data
ierr = init_timeslices(timeslice_path)
if (ierr /= 0) call hydro_stop("No timeslice files available when initializing diversions")

allocate(diversions(ndivs))
do i = 1, ndivs
associate (div => diversions(i))
div = diversion_t('', '', '', -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1)

ierr = 0
!ierr = ierr + nf_get_var1(ncid, name_vid, i, div%name) !! can't read string with Fortran :-(

ierr = ierr + nf90_get_var(ncid, da_src_vid, div%da_src, start=(/i/), count=(/15/))
div%persisted_flow_src = get_flow_for_gage(div%da_src)
ierr = ierr + nf90_get_var(ncid, da_dest_vid, div%da_dest, start=(/i/), count=(/15/))
div%persisted_flow_dest = get_flow_for_gage(div%da_dest)

ierr = ierr + nf90_get_var(ncid, type_div_vid, div%type_div, start=(/i/))
ierr = ierr + nf90_get_var(ncid, type_src_vid, div%type_src, start=(/i/))
ierr = ierr + nf90_get_var(ncid, type_dest_vid, div%type_dest, start=(/i/))
ierr = ierr + nf90_get_var(ncid, id_src_vid, div%id_src, start=(/i/))
ierr = ierr + nf90_get_var(ncid, id_dest_vid, div%id_dest, start=(/i/))
ierr = ierr + nf90_get_var(ncid, capacity_vid, div%capacity, start=(/i/))
ierr = ierr + nf90_get_var(ncid, fraction_vid, div%fraction, start=(/i/))
ierr = ierr + nf90_get_var(ncid, lookback_vid, div%lookback, start=(/i/))

if (ierr /= 0) call hydro_stop("Error occurred reading diversion variables from diversion file")
end associate
end do

end if
end if

end subroutine

subroutine calculate_diversion(src_link_in, qlink_src_in, diversion_quantity_out, diversion_quantity_in)
integer(kind=int64), intent(in) :: src_link_in
! integer(kind=int64), intent(out) :: dst_out
real, intent(in) :: qlink_src_in
real, intent(out) :: diversion_quantity_out, diversion_quantity_in

integer :: i

diversion_quantity_out = 0
diversion_quantity_in = 0

! bail if we're inactive
if (.not. diversions_active) return

! link to gage
! look to see what type of diversion it is
! call into sub-procedure to handle type=1, type=2, type=3, etc

do i = 1, ndivs
if (src_link_in == diversions(i)%id_src) then
if (diversions(i)%type_div /= 3) then
print free, "!!! UNSUPPORTED DIVERSION TYPE (", diversions(i)%type_div, "), skipping"
else
call gage_assisted_diversion(src_link_in, diversions(i), qlink_src_in, diversion_quantity_out)
! dst_out = diversions(i)%id_dest
end if
end if

if (src_link_in == diversions(i)%id_dest) then
if (diversions(i)%type_div /= 3) then
print free, "!!! UNSUPPORTED DIVERSION TYPE (", diversions(i)%type_div, "), skipping"
else
diversion_quantity_in = diversions(i)%persisted_flow_dest
end if
end if
end do

! subtract dst_out from source gage

end subroutine

subroutine gage_assisted_diversion(src_link, diversion, qlink_src, div_gage_flow)
integer(kind=int64), intent(in) :: src_link
type(diversion_t), intent(in) :: diversion
real, intent(in) :: qlink_src
real, intent(out) :: div_gage_flow

real :: fraction

! This is the so-called "Type 3" diversion. We take the observed flow from div_gage,
! and subtract it from the upstream qlink_src, if it's a valid flow (not-NaN).
!
! If it's not a valid flow, we try to use the Fraction property of the diversion,
! and if -that's- not available, we just leave the flow untouched.

div_gage_flow = diversion%persisted_flow_dest
if (ieee_is_nan(div_gage_flow)) then
fraction = diversion%fraction
if (fraction == -1) then
print free, "WARNING: No fractional diversion value specified for diversion at gage '" // trim(adjustl(diversion%da_dest)) // "', skipping"
fraction = 0
else
print free, "INFO: No gage discharge available for diversion '" // trim(adjustl(diversion%da_dest)) // "', using fixed fractional diversion of", fraction
end if
div_gage_flow = qlink_src * fraction
end if
end subroutine

end module
Loading

0 comments on commit 2682cee

Please sign in to comment.