Skip to content

Commit

Permalink
Add support for adding diversion to destination
Browse files Browse the repository at this point in the history
  • Loading branch information
@rcabell
rcabell committed Aug 1, 2024
1 parent bcb6ed7 commit ac60260
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 33 deletions.
60 changes: 28 additions & 32 deletions src/Routing/Diversions/module_diversions.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,14 +106,17 @@ subroutine init_diversions(diversions_file, timeslice_path)
end subroutine
subroutine calculate_diversion(src_link, dst_link, qlink_src_previous, qlink_src_current, qlink_dst)
integer(kind=int64), intent(in) :: src_link
integer(kind=int64), intent(in) :: dst_link
real, intent(inout) :: qlink_src_previous, qlink_src_current
real, intent(inout) :: qlink_dst
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

Expand All @@ -122,23 +125,25 @@ subroutine calculate_diversion(src_link, dst_link, qlink_src_previous, qlink_src
! call into sub-procedure to handle type=1, type=2, type=3, etc

do i = 1, ndivs
if (src_link == diversions(i)%id_src) then
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, diversions(i), qlink_src_previous, qlink_src_current)
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_prev, qlink_src_curr)
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(inout) :: qlink_src_prev, qlink_src_curr
real, intent(in) :: qlink_src
real, intent(out) :: div_gage_flow

real :: div_gage_flow, fraction
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).
Expand All @@ -155,32 +160,23 @@ subroutine gage_assisted_diversion(src_link, diversion, qlink_src_prev, qlink_sr
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_curr * fraction
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_prev) then
qlink_src_prev = qlink_src_prev - div_gage_flow
else
qlink_src_prev = 0
#ifdef HYDRO_D
print free, "DEBUG WARNING: diverted flow (", div_gage_flow, ") exceeds total flow, zeroing."
#endif
end if

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
! #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
36 changes: 35 additions & 1 deletion src/Routing/module_channel_routing.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1712,6 +1712,11 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, &
real, allocatable,dimension(:) :: tmpAssimilatedValue
character(len=256), allocatable,dimension(:) :: tmpAssimilatedSourceFile

! diversions
integer(kind=int64) :: div_dest
real :: div_qty
character(*), parameter :: free = '(*(g0,1x))'

#ifdef MPP_LAND
if(my_id .eq. io_id) then
#endif
Expand Down Expand Up @@ -1961,7 +1966,36 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, &
! QLateral(k), DTRT_CH, So(k), CHANLEN(k), &
! MannN(k), ChSSlp(k), Bw(k), Tw(k) )

call calculate_diversion(LINKID(k), -1_int64, Qup, Quc, tmpQLINK(k,2))
! HANDLE DIVERSIONS

call calculate_diversion(LINKID(k), div_dest, Quc, div_qty)
if (div_qty /= 0) then
! remove from upstream
#ifdef HYDRO_D
print free, "DEBUG: diverting", div_qty, "of", Quc, "from link id =", LINKID(k)
if (div_qty > Quc) &
print free, "DEBUG WARNING: diverted flow (", div_qty, ") exceeds total flow, zeroing."
#endif

! add to downstream if needed
if (div_dest /= -99) then
! find destination (it's not necessarily tmpQLINK(k,2)))
do kk = 1,NLINKSL
if (LINKID(kk) == div_dest) then
#ifdef HYDRO_D
print free, "Found diversion destination ID", LINKID(kk),"in local processor at index", kk, ", replacing", tmpQLINK(kk,2), "with", div_qty
#endif
tmpQLINK(kk,2) = div_qty
exit !do loop
end if
end do
Quc = max(0.0, Quc - div_qty)
Qup = max(0.0, Qup - div_qty)
end if
end if
#ifdef WRF_HYDRO_NUDGING
call nudge_apply_upstream_muskingumCunge( Qup, Quc, nudge(k), k )
Expand Down

0 comments on commit ac60260

Please sign in to comment.