From ac602601d483a51b4b92de544890ab9a96e44ce0 Mon Sep 17 00:00:00 2001 From: Ryan Cabell Date: Thu, 1 Aug 2024 15:11:07 -0600 Subject: [PATCH] Add support for adding diversion to destination --- src/Routing/Diversions/module_diversions.F | 60 ++++++++++------------ src/Routing/module_channel_routing.F | 36 ++++++++++++- 2 files changed, 63 insertions(+), 33 deletions(-) diff --git a/src/Routing/Diversions/module_diversions.F b/src/Routing/Diversions/module_diversions.F index 3fec5dfcf..35c2f7998 100644 --- a/src/Routing/Diversions/module_diversions.F +++ b/src/Routing/Diversions/module_diversions.F @@ -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 @@ -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). @@ -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 diff --git a/src/Routing/module_channel_routing.F b/src/Routing/module_channel_routing.F index 9e0d4d773..e30b522de 100644 --- a/src/Routing/module_channel_routing.F +++ b/src/Routing/module_channel_routing.F @@ -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 @@ -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 )