Add initial support for gage-assisted diversions in channel routing

CMakeLists.txt

@@ -6,7 +6,7 @@ endif()
 
 
 # set project name and version numbers
-project (WRF_Hydro LANGUAGES Fortran)
+project (WRF_Hydro LANGUAGES Fortran C)
 set (WRF_Hydro_VERSION_MAJOR 5)
 set (WRF_Hydro_VERSION_MINOR 3)
 set (WRF_Hydro_VERSION_PATCH 0)
@@ -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")

src/CMakeLists.txt

@@ -8,6 +8,7 @@ add_subdirectory("Routing/Overland")
 add_subdirectory("Routing/Subsurface")
 #add_subdirectory("Routing/Groundwater")
 add_subdirectory("Routing/Reservoirs")
+add_subdirectory("Routing/Diversions")
 add_subdirectory("Data_Rec")
 add_subdirectory("Routing")
 add_subdirectory("HYDRO_drv")
@@ -103,7 +104,9 @@ 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 (WRF_HYDRO_NUDGING STREQUAL "1")
 
  # bash commands to copy namelists to the Run directory

src/Data_Rec/module_namelist.F

@@ -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=""
@@ -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, &
@@ -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

src/Data_Rec/module_namelist_inc.F

@@ -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=""

src/HYDRO_drv/module_HYDRO_drv.F

@@ -29,6 +29,8 @@ 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
@@ -1580,6 +1582,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

src/OrchestratorLayer/config.f90

@@ -123,6 +123,7 @@ module config_base
  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=""
@@ -169,7 +170,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
@@ -499,6 +500,7 @@ subroutine init_namelist_rt_field(did)
  logical :: compound_channel
  integer :: channel_loss_option = 0
  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=""
@@ -573,7 +575,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, route_lake_f, &
+ route_topo_f,route_chan_f,route_link_f, compound_channel, channel_loss_option, 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, &
@@ -801,6 +803,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

src/Routing/CMakeLists.txt

@@ -17,6 +17,8 @@ add_library(hydro_routing STATIC
  Noah_distr_routing_subsurface.F
 )
 
+add_dependencies(hydro_routing hydro_routing_diversions)
+
 target_link_libraries(hydro_routing PUBLIC hydro_mpp)
 target_link_libraries(hydro_routing PUBLIC hydro_utils)
 target_link_libraries(hydro_routing PUBLIC hydro_orchestrator)
@@ -31,4 +33,5 @@ target_link_libraries(hydro_routing PUBLIC hydro_routing_reservoirs_levelpool)
 target_link_libraries(hydro_routing PUBLIC hydro_routing_reservoirs_hybrid)
 target_link_libraries(hydro_routing PUBLIC hydro_data_rec)
 target_link_libraries(hydro_routing PUBLIC hydro_routing_reservoirs_rfc)
+target_link_libraries(hydro_routing PUBLIC hydro_routing_diversions)
 

src/Routing/Diversions/CMakeLists.txt

@@ -0,0 +1,10 @@
+add_library(hydro_routing_diversions STATIC
+ module_diversions.F
+ module_diversions_timeslice.F
+)
+
+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)

src/Routing/Diversions/module_diversions.F

@@ -0,0 +1,182 @@
+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, dst_link_out, qlink_src_in, diversion_quantity_out)
+ integer(kind=int64), intent(in) :: src_link_in
+ integer(kind=int64), intent(out) :: dst_link_out
+ real, intent(in) :: qlink_src_in
+ real, intent(out) :: diversion_quantity_out
+
+ integer :: i
+
+ dst_link_out = -99
+ diversion_quantity_out = 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_link_out = diversions(i)%id_dest
+ end if
+ end if
+ end do
+
+ 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
+
+ ! div_gage_flow = min(div_gage_flow, diversion%capacity) ! don't divert more than diversion can handle
+
+! #ifdef HYDRO_D
+! if (div_gage_flow /= 0) &
+! print free, "DEBUG: diverting", div_gage_flow, "of", qlink_src_curr, "from link id =", src_link
+! #endif
+! if (div_gage_flow <= qlink_src_curr) then
+! qlink_src_curr = qlink_src_curr - div_gage_flow
+! else
+! qlink_src_curr = 0
+! #ifdef HYDRO_D
+! print free, "DEBUG WARNING: diverted flow (", div_gage_flow, ") exceeds total flow, zeroing."
+! #endif
+! end if
+ end subroutine
+
+end module